Introduction

This code uses the variancePartition package to quantify the percent of variance in the data explained by some biological/technical factors.

Load packages and normalized data

library(variancePartition)

Loading required package: ggplot2

Loading required package: limma

Loading required package: foreach

Loading required package: scales

Loading required package: Biobase

Loading required package: BiocGenerics

Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following object is masked from 'package:limma':

    plotMA

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    anyDuplicated, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter,
    Find, get, grep, grepl, intersect, is.unsorted, lapply, Map,
    mapply, match, mget, order, paste, pmax, pmax.int, pmin,
    pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
    setdiff, sort, table, tapply, union, unique, unsplit, which,
    which.max, which.min

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Attaching package: 'variancePartition'

The following object is masked from 'package:limma':

    classifyTestsF

library(RUVSeq)

Loading required package: EDASeq

Loading required package: ShortRead

Loading required package: BiocParallel

Loading required package: Biostrings

Loading required package: S4Vectors

Loading required package: stats4

Attaching package: 'S4Vectors'

The following object is masked from 'package:base':

    expand.grid

Loading required package: IRanges

Loading required package: XVector

Attaching package: 'Biostrings'

The following object is masked from 'package:base':

    strsplit

Loading required package: Rsamtools

Loading required package: GenomeInfoDb

Loading required package: GenomicRanges

Loading required package: GenomicAlignments

Loading required package: SummarizedExperiment

Loading required package: DelayedArray

Loading required package: matrixStats

Attaching package: 'matrixStats'

The following objects are masked from 'package:Biobase':

    anyMissing, rowMedians

Attaching package: 'DelayedArray'

The following objects are masked from 'package:matrixStats':

    colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

The following object is masked from 'package:Biostrings':

    type

The following objects are masked from 'package:base':

    aperm, apply, rowsum

Loading required package: edgeR

library(edgeR)
library(dplyr)

Attaching package: 'dplyr'

The following object is masked from 'package:ShortRead':

    id

The following objects are masked from 'package:GenomicAlignments':

    first, last

The following object is masked from 'package:matrixStats':

    count

The following objects are masked from 'package:GenomicRanges':

    intersect, setdiff, union

The following object is masked from 'package:GenomeInfoDb':

    intersect

The following objects are masked from 'package:Biostrings':

    collapse, intersect, setdiff, setequal, union

The following object is masked from 'package:XVector':

    slice

The following objects are masked from 'package:IRanges':

    collapse, desc, intersect, setdiff, slice, union

The following objects are masked from 'package:S4Vectors':

    first, intersect, rename, setdiff, setequal, union

The following object is masked from 'package:Biobase':

    combine

The following objects are masked from 'package:BiocGenerics':

    combine, intersect, setdiff, union

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

#raw(filtered) counts
counts_upperquartile <- readRDS("data/filtered_counts.rds")$counts
#normalized and filtered counts
filtered_upperquartile <- readRDS("data/norm_filtered_counts.rds")
#metadata
sampleinfo <- readRDS("data/Sample.info.RNAseq.reordered.rds")

Fitting unwanted variation

Before fitting the linear mixed model to quantify contributions of different variables to variation in the dataset, we are also interested in fitting one factor of unwanted variation to include in the model using RUVg. For the RUVg control genes, we use the top 100 (~1% of the genes considered) LEAST variable genes as control genes.

#input data consists of raw filtered data (filtered for lowly expressed genes)
#compute CV (stdev/mean) and rank from least to most; pick 100 least variable
cv <- apply(counts_upperquartile, 1, function(x) sd(x)/mean(x))
least_var_genes <- names(head(sort(cv), 100))

Apply RUVseq with the set of least variable genes to determine a single factor of unwanted variation.

#load data into expressionset
set <- newSeqExpressionSet(as.matrix(counts_upperquartile),phenoData = data.frame(sampleinfo, row.names=colnames(counts_upperquartile)))
set

SeqExpressionSet (storageMode: lockedEnvironment)
assayData: 10486 features, 17 samples 
  element names: counts, normalizedCounts, offset 
protocolData: none
phenoData
  sampleNames: 18855_3_S 19160_3_S ... 18855_3_U (17 total)
  varLabels: Sample_ID Individual ... LibSize (8 total)
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:  

#normalization
set <- betweenLaneNormalization(x = set, which = "upper")
#run RUVg
set1 <- RUVg(set, least_var_genes, k=1)
sample_info <- pData(set1)
sample_info$Replicate <- as.factor(sample_info$Replicate)
sample_info$LibraryPrepBatch <- as.factor(sample_info$LibraryPrepBatch)

Specify variables to consider in LMM

Here, we specify two separate models (one including the single factor of unwanted variation computed above and one without this factor).

# Specify variables to consider
form <- ~ (1|Individual) + (1|treatment) + (1|Replicate) + (1|LibraryPrepBatch) + W_1

# No ruv
form_no_ruv <- ~ (1|Individual) + (1|treatment) + (1|Replicate) + (1|LibraryPrepBatch)

Partition the variance and plot the result

# Fit model and extract results
# 1) fit linear mixed model on gene expression
# If categorical variables are specified,
# a linear mixed model is used
# If all variables are modeled as fixed effects,
# a linear model is used
# each entry in results is a regression model fit on a single gene
# 2) extract variance fractions from each model fit
# for each gene, returns fraction of variation attributable
# to each variable
# Interpretation: the variance explained by each variables
# after correcting for all other variables
varPart_1_int <- fitExtractVarPartModel( filtered_upperquartile, form, sample_info )

Dividing work into 100 chunks...

Total: 29 s

# sort variables (i.e. columns) by median fraction
# of variance explained
vp <- sortCols(varPart_1_int)
# violin plot of contribution of each variable to total variance
plotVarPart( vp )

Past versions of function_1-1.png

#noRuv
VarPart_noRUV <- fitExtractVarPartModel( filtered_upperquartile, form_no_ruv, sample_info )

Dividing work into 100 chunks...

Total: 28 s

vp_noRUV <- sortCols(VarPart_noRUV)
plotVarPart( vp_noRUV )

Past versions of function_1-2.png

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
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] stats4    parallel  stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] dplyr_1.0.2                 RUVSeq_1.18.0              
 [3] edgeR_3.26.5                EDASeq_2.18.0              
 [5] ShortRead_1.42.0            GenomicAlignments_1.20.1   
 [7] SummarizedExperiment_1.14.1 DelayedArray_0.10.0        
 [9] matrixStats_0.57.0          Rsamtools_2.0.0            
[11] GenomicRanges_1.36.1        GenomeInfoDb_1.20.0        
[13] Biostrings_2.52.0           XVector_0.24.0             
[15] IRanges_2.18.3              S4Vectors_0.22.1           
[17] BiocParallel_1.18.1         variancePartition_1.14.1   
[19] Biobase_2.44.0              BiocGenerics_0.30.0        
[21] scales_1.1.1                foreach_1.4.4              
[23] limma_3.40.6                ggplot2_3.3.3              

loaded via a namespace (and not attached):
 [1] minqa_1.2.4            colorspace_2.0-0       hwriter_1.3.2         
 [4] ellipsis_0.3.1         colorRamps_2.3         rprojroot_2.0.2       
 [7] fs_1.3.1               farver_2.0.3           bit64_0.9-7           
[10] AnnotationDbi_1.46.0   codetools_0.2-16       splines_3.6.1         
[13] R.methodsS3_1.8.1      doParallel_1.0.14      DESeq_1.36.0          
[16] geneplotter_1.62.0     knitr_1.23             workflowr_1.6.2       
[19] nloptr_1.2.2.1         pbkrtest_0.4-7         annotate_1.62.0       
[22] R.oo_1.24.0            compiler_3.6.1         httr_1.4.2            
[25] Matrix_1.2-18          later_1.1.0.1          htmltools_0.5.0       
[28] prettyunits_1.1.1      tools_3.6.1            gtable_0.3.0          
[31] glue_1.4.2             GenomeInfoDbData_1.2.1 reshape2_1.4.3        
[34] Rcpp_1.0.5             vctrs_0.3.6            gdata_2.18.0          
[37] nlme_3.1-140           rtracklayer_1.44.0     iterators_1.0.12      
[40] xfun_0.8               stringr_1.4.0          lme4_1.1-21           
[43] lifecycle_0.2.0        gtools_3.8.1           XML_3.98-1.20         
[46] zlibbioc_1.30.0        MASS_7.3-52            aroma.light_3.14.0    
[49] hms_0.5.3              promises_1.1.1         RColorBrewer_1.1-2    
[52] yaml_2.2.1             memoise_1.1.0          biomaRt_2.40.1        
[55] latticeExtra_0.6-28    stringi_1.4.6          RSQLite_2.1.1         
[58] genefilter_1.66.0      GenomicFeatures_1.36.3 caTools_1.17.1.2      
[61] boot_1.3-23            rlang_0.4.10           pkgconfig_2.0.3       
[64] bitops_1.0-6           evaluate_0.14          lattice_0.20-41       
[67] purrr_0.3.4            labeling_0.4.2         bit_4.0.4             
[70] tidyselect_1.1.0       plyr_1.8.6             magrittr_2.0.1        
[73] R6_2.5.0               gplots_3.0.1.1         generics_0.0.2        
[76] DBI_1.1.0              pillar_1.4.7           whisker_0.3-2         
[79] withr_2.3.0            survival_2.44-1.1      RCurl_1.98-1.1        
[82] tibble_3.0.4           crayon_1.3.4           KernSmooth_2.23-15    
[85] rmarkdown_1.13         progress_1.2.2         locfit_1.5-9.1        
[88] grid_3.6.1             blob_1.2.0             git2r_0.26.1          
[91] digest_0.6.27          xtable_1.8-4           httpuv_1.5.1          
[94] R.utils_2.9.0          munsell_0.5.0