Hi

I am comparing two different types of sequencing method for RNAseq. I am using the same set of samples (3 control, 3 treated) to prepare library using illumina mRNA kit and lexogen 3 prime library kit. The idea is to compare both library preparation method with the aim of using the three prime one for future experiments. we were hoping lexogen will generate roughly similar numbers of DGE, but contrary to our expectation three prime generated only 60% of the DGE identified in the illumina library. I have identified the genes that were not called significant by DEseq2 from lexogen dataset. How can I find why the genes are not called DEs? Any suggesion is much appreciated. The basic code is as follows.

Many thanks Nurun

dds <- DESeqDataSetFromMatrix(countData = Standard,
                              colData = SampleData.standard,
                              design = ~Treatment+Subject)

dds <- estimateSizeFactors(dds)

vsd <- vst(dds)

plotPCA(vsd, intgroup = "Treatment")

dds <- DESeq(dds)

res <- results(dds, lfcThreshold = 2, alpha = 0.05, contrast=c("Treatment","LPS","Control"))

summary(res)

sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

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

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

other attached packages:
 [1] forcats_0.4.0               stringr_1.4.0               dplyr_0.8.0.1              
 [4] purrr_0.3.2                 readr_1.3.1                 tidyr_0.8.3                
 [7] tibble_2.1.1                tidyverse_1.2.1             ggforce_0.2.2              
[10] VennDiagram_1.6.20          futile.logger_1.4.3         org.Hs.eg.db_3.8.2         
[13] AnnotationDbi_1.46.0        DESeq2_1.24.0               SummarizedExperiment_1.14.0
[16] DelayedArray_0.10.0         BiocParallel_1.18.0         matrixStats_0.54.0         
[19] Biobase_2.44.0              GenomicRanges_1.36.0        GenomeInfoDb_1.20.0        
[22] IRanges_2.18.0              S4Vectors_0.22.0            BiocGenerics_0.30.0        
[25] reshape2_1.4.3              ggplot2_3.1.1              

loaded via a namespace (and not attached):
 [1] nlme_3.1-140           bitops_1.0-6           lubridate_1.7.4        bit64_0.9-7           
 [5] httr_1.4.0             RColorBrewer_1.1-2     tools_3.6.1            backports_1.1.4       
 [9] R6_2.4.0               rpart_4.1-15           Hmisc_4.2-0            DBI_1.0.0             
[13] lazyeval_0.2.2         colorspace_1.4-1       nnet_7.3-12            withr_2.1.2           
[17] tidyselect_0.2.5       gridExtra_2.3          bit_1.1-14             compiler_3.6.1        
[21] cli_1.1.0              rvest_0.3.4            formatR_1.6            htmlTable_1.13.1      
[25] xml2_1.2.0             scales_1.0.0           checkmate_1.9.3        genefilter_1.66.0     
[29] digest_0.6.18          foreign_0.8-71         XVector_0.24.0         base64enc_0.1-3       
[33] pkgconfig_2.0.2        htmltools_0.3.6        readxl_1.3.1           htmlwidgets_1.3       
[37] rlang_0.3.4            rstudioapi_0.10        RSQLite_2.1.1          farver_1.1.0          
[41] generics_0.0.2         jsonlite_1.6           acepack_1.4.1          RCurl_1.95-4.12       
[45] magrittr_1.5           GenomeInfoDbData_1.2.1 Formula_1.2-3          Matrix_1.2-17         
[49] Rcpp_1.0.1             munsell_0.5.0          stringi_1.4.3          MASS_7.3-51.4         
[53] zlibbioc_1.30.0        plyr_1.8.4             blob_1.1.1             crayon_1.3.4          
[57] lattice_0.20-38        haven_2.1.0            splines_3.6.1          annotate_1.62.0       
[61] hms_0.4.2              locfit_1.5-9.1         knitr_1.22             pillar_1.4.0          
[65] geneplotter_1.62.0     futile.options_1.0.1   XML_3.98-1.19          glue_1.3.1            
[69] latticeExtra_0.6-28    modelr_0.1.4           lambda.r_1.2.3         data.table_1.12.2     
[73] tweenr_1.0.1           cellranger_1.1.0       gtable_0.3.0           polyclip_1.10-0       
[77] assertthat_0.2.1       xfun_0.6               xtable_1.8-4           broom_0.5.2           
[81] survival_2.44-1.1      memoise_1.1.0          cluster_2.1.0 

Assuming both protocols were performed well (skilled person at the bench, enough RNA input, comparable RIN scores, etc), my guess would be that the lower counts / gene you are getting in the quantseq prep should probably explain the majority of the differences ... features with more counts are more easily called as DGE.

If I were poking around to these data to better sort this out, here are some things I might try:

1. I think I'd rather analyze these data in independently from each other, ie. create two separate DESeqDataSets, one for each sequencing protocol used, and run the analysis / DGE pipeline separately for each. Perhaps the dispersion estimate is better fit separately for the protocols. You can then also see differences between plotDispEsts on the quantseq vs the other Illumina protocol you are using.

   This should also allow you to light up the genes on the dispersion plot that are not identified in as DGE in each protocol ... do the tend to be lowly expressed, have high dispersion, both?

2. Compare the MA plots for the replicates within each protocol ... anything funky there?

3. You only have 3 vs 3, so perhaps not that informative, but you should verify that the replicates at least group together reasonably enough via some clustering approach or PCA.

4. Plot the logFC's calculated from each library prep vs each other. Are the results largely concordant (ie. on the 45° line?) and it's just the statistical power that's lower in the quantseq data? If so, perhaps sequencing the quantseq library deeper (with, perhaps, requring more input RNA to begin with).

Also: are you using the UMI kits with your quantseq prep? Given the tagging nature of the protocol, you definitely should consider using them to remove PCR duplicates ... which also should make you think about how much PCR you're doing with the illumina prep. It's likely less of a problem there, but still ... these are all good things to know.

Also: I just noticed you are testing against a threshold, ie. results(dds, lfcThreshold = 2, ...). Start your initial analysis by dropping the lfcThreshold parameter (leaving it set to 0).

What will people downstream do with the list of DE genes? If you do pathway analysis or look at gene enrichment, and you get the same pathways and gene groups returned with both preps, the it might be okay that you lost some genes with lower expression.