Apply stm to GTEx data

Introduction

Gtex Data from Yang Li. Data consist of per-base read counts from “GeuvadisLCL”, “GTEXAdipose”, “WholeBlood”.

Gene: RPS13. It has 6 exons and 10 splicing variants.

load('~/RPS13fit.RData')
library(stm)

Full Data

Fit stm to the full dataset with \(K=3\).

lf = poisson2multinom(t(fit_stm_bmsm_RPS13$qf$Ef),fit_stm_bmsm_RPS13$ql$El)
for(i in 1:3){
  barplot(lf$L[,i],main=paste0('Estimate of loading ',i),xlab='',ylab = '',ylim = c(0,max(lf$L)),col=c(rep(2,88),rep(3,226),rep(4,927)))
}

Past versions of unnamed-chunk-2-1.png

Past versions of unnamed-chunk-2-2.png

Past versions of unnamed-chunk-2-3.png

for(i in 1:3){
  barplot(lf$FF[,i],main=paste0('Estimate of factor ',i),xlab='',ylab = '',ylim = c(0,max(lf$FF)))
}

Past versions of unnamed-chunk-2-4.png

Past versions of unnamed-chunk-2-5.png

Past versions of unnamed-chunk-2-6.png

Remove WholeBlood

Fit stm to the dataset removing “WholeBlood” tissue with \(K=3\).

lf = poisson2multinom(t(fit_stm_bmsm_RPS13_2$qf$Ef),fit_stm_bmsm_RPS13_2$ql$El)
plot(lf$FF[,1],ylim=range(lf$FF),type='l',xlab = 'base',ylab='Intensity',main='Factors')
lines(lf$FF[,2],col=2)
lines(lf$FF[,3],col=4)

Past versions of unnamed-chunk-3-1.png

for(i in 1:3){
  barplot(lf$L[,i],main=paste0('Estimated loading ',i),ylab = '',ylim = c(0,max(lf$L)),col=c(rep(2,88),rep(3,226)))
}

Past versions of unnamed-chunk-3-2.png

Past versions of unnamed-chunk-3-3.png

Past versions of unnamed-chunk-3-4.png

Fit Xing’s stm method to the dataset removing “WholeBlood” tissue with \(K=3\).

plot(xing_rps13$phi[1,],type='l',ylim = range(xing_rps13$phi),xlab = 'base',ylab='Intensity',main='Factors')
lines(xing_rps13$phi[2,],type='l',col=2)
lines(xing_rps13$phi[3,],type='l',col=4)

Past versions of unnamed-chunk-4-1.png

Fit stm to the dataset removing “WholeBlood” tissue with \(K=10\).

lf = poisson2multinom(t(fit_stm_bmsm_RPS13_3$qf$Ef),fit_stm_bmsm_RPS13_3$ql$El)


for(i in 1:10){
  barplot(lf$L[,i],main=paste0('Estimate of loading ',i),xlab='',ylab = '',ylim = c(0,max(lf$L)),col=c(rep(2,88),rep(3,226)))
}

Past versions of unnamed-chunk-5-1.png

Past versions of unnamed-chunk-5-2.png

Past versions of unnamed-chunk-5-3.png

Past versions of unnamed-chunk-5-4.png

Past versions of unnamed-chunk-5-5.png

Past versions of unnamed-chunk-5-6.png

Past versions of unnamed-chunk-5-7.png

Past versions of unnamed-chunk-5-8.png

Past versions of unnamed-chunk-5-9.png

Past versions of unnamed-chunk-5-10.png

for(i in 1:10){
  barplot(lf$FF[,i],main=paste0('Estimate of factor ',i),xlab='',ylab = '',ylim = c(0,max(lf$FF)))
}

Past versions of unnamed-chunk-5-11.png

Past versions of unnamed-chunk-5-12.png

Past versions of unnamed-chunk-5-13.png

Past versions of unnamed-chunk-5-14.png

Past versions of unnamed-chunk-5-15.png

Past versions of unnamed-chunk-5-16.png

Past versions of unnamed-chunk-5-17.png

Past versions of unnamed-chunk-5-18.png

Past versions of unnamed-chunk-5-19.png

Past versions of unnamed-chunk-5-20.png

GTEx adipose tissue

Fit stm to the dataset using only adipose tissue with \(K=5\).

lf = poisson2multinom(t(fit_stm_bmsm_RPS13_4$qf$Ef),fit_stm_bmsm_RPS13_4$ql$El)
plot(lf$FF[,1],ylim=range(lf$FF),type='l',xlab = 'base',ylab='Intensity',main='Factors')
lines(lf$FF[,2],col=2)
lines(lf$FF[,3],col=4)
lines(lf$FF[,4],col=3)
lines(lf$FF[,5],col=6)

Past versions of unnamed-chunk-6-1.png

for(i in 1:5){
  barplot(lf$L[,i],main=paste0('Estimate of loading ',i),xlab='',ylab = '',ylim = c(0,max(lf$L)))
}

Past versions of unnamed-chunk-7-1.png

Past versions of unnamed-chunk-7-2.png

Past versions of unnamed-chunk-7-3.png

Past versions of unnamed-chunk-7-4.png

Past versions of unnamed-chunk-7-5.png

for(i in 1:5){
  barplot(lf$FF[,i],main=paste0('Estimate of factor ',i),xlab='',ylab = '',ylim = c(0,max(lf$FF)))
}

Past versions of unnamed-chunk-7-6.png

Past versions of unnamed-chunk-7-7.png

Past versions of unnamed-chunk-7-8.png

Past versions of unnamed-chunk-7-9.png

Past versions of unnamed-chunk-7-10.png

Fit stm with smashgen using only adipose tissue with \(K=3\).

lf = poisson2multinom(t(fit_stm_sg_RPS13_4$qf$Ef),fit_stm_sg_RPS13_4$ql$El)
plot(lf$FF[,1],ylim=range(lf$FF),type='l',xlab = 'base',ylab='Intensity',main='Factors')
lines(lf$FF[,2],col=2)
lines(lf$FF[,3],col=4)

Past versions of unnamed-chunk-8-1.png

for(i in 1:3){
  barplot(lf$L[,i],main=paste0('Estimate of loading ',i),xlab='',ylab = '',ylim = c(0,max(lf$L)))
}

Past versions of unnamed-chunk-9-1.png

Past versions of unnamed-chunk-9-2.png

Past versions of unnamed-chunk-9-3.png

for(i in 1:3){
  barplot(lf$FF[,i],main=paste0('Estimate of factor ',i),xlab='',ylab = '',ylim = c(0,max(lf$FF)))
}

Past versions of unnamed-chunk-9-4.png

Past versions of unnamed-chunk-9-5.png

Past versions of unnamed-chunk-9-6.png

Estimated Nugget effect:

fit_stm_sg_RPS13_4$nugget$nugget_f

[1] 0.03478134 0.12772478 0.04643110

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] stm_1.0.0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.2        compiler_3.5.1    later_0.7.5      
 [4] git2r_0.26.1      workflowr_1.5.0   bitops_1.0-6     
 [7] iterators_1.0.10  tools_3.5.1       digest_0.6.18    
[10] evaluate_0.12     lattice_0.20-38   Matrix_1.2-15    
[13] foreach_1.4.4     yaml_2.2.0        parallel_3.5.1   
[16] smashr_1.2-7      stringr_1.3.1     knitr_1.20       
[19] fs_1.3.1          caTools_1.17.1.1  rprojroot_1.3-2  
[22] grid_3.5.1        data.table_1.12.0 glue_1.3.0       
[25] R6_2.3.0          rmarkdown_1.10    mixsqp_0.2-2     
[28] ashr_2.2-39       magrittr_1.5      whisker_0.3-2    
[31] backports_1.1.2   promises_1.0.1    codetools_0.2-15 
[34] htmltools_0.3.6   MASS_7.3-51.1     httpuv_1.4.5     
[37] wavethresh_4.6.8  stringi_1.2.4     doParallel_1.0.14
[40] pscl_1.5.2        truncnorm_1.0-8   NNLM_0.4.2       
[43] SQUAREM_2017.10-1 logitnorm_0.8.37