Question: RNAseq analysis DeSeq2 strange MAplot
gravatar for alice.checcucci
7 days ago by
alice.checcucci0 wrote:


I am using DESeq2 for DEG analysis in my RNA-Seq experiment. I have 2 bacterial strains: one is the wild type and the other is a mutant (it lacks one plasmid). When I look for differentially expressed genes, plotting the results with MAplot I obtain a strange plot, which seems not to have the "classical" symmetrical shape (like the one that is reported in the tutorial).

Here is the codes:

txi.tx <- tximport(files, type = "salmon", txOut = TRUE)

samples <- data.frame(species = c("2011","2011","2011","delta","delta","delta")) # sample metadata 
ddsTxi <- DESeqDataSetFromTximport(txi.tx, colData = samples, design = ~species) 

res <- t(apply(counts(ddsTxi),1, function(x)by(x, samples$species, sum)))
keep <- res[,2] > 0
ddsTxi <- ddsTxi[keep,]

keep <- rowSums(counts(ddsTxi)) >= 10
ddsTxi <- ddsTxi[keep,]

ddsTxi$species<- factor(ddsTxi$species, levels = c("2011","delta"))

ddsTxi <- DESeq(ddsTxi, fitType = "local")
res <- lfcShrink(ddsTxi, coef = "species_delta_vs_2011",  type = "normal")

plotMA(res, ylim=c(-5,5)) MAplot

R version 3.4.3 (2017-11-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 8 (jessie)

Matrix products: default
BLAS: /usr/lib/libblas/
LAPACK: /usr/lib/lapack/

 [1] LC_CTYPE=it_IT.UTF-8       LC_NUMERIC=C               LC_TIME=it_IT.UTF-8        LC_COLLATE=it_IT.UTF-8    
 [5] LC_MONETARY=it_IT.UTF-8    LC_MESSAGES=it_IT.UTF-8    LC_PAPER=it_IT.UTF-8       LC_NAME=C                 

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

other attached packages:
 [1] DESeq2_1.18.1              SummarizedExperiment_1.8.1 DelayedArray_0.4.1         matrixStats_0.54.0        
 [5] Biobase_2.38.0             GenomicRanges_1.30.3       GenomeInfoDb_1.14.0        Biostrings_2.46.0         
 [9] XVector_0.18.0             IRanges_2.12.0             S4Vectors_0.16.0           BiocGenerics_0.24.0       
[13] tximport_1.6.0             BiocInstaller_1.28.0       ggplot2_3.1.0              phyloseq_1.22.3           

loaded via a namespace (and not attached):
 [1] jsonlite_1.6           bit64_0.9-7            splines_3.4.3          foreach_1.4.4          Formula_1.2-3         
 [6] latticeExtra_0.6-28    blob_1.1.0             GenomeInfoDbData_1.0.0 pillar_1.3.0           RSQLite_2.0           
[11] backports_1.1.2        lattice_0.20-35        digest_0.6.18          RColorBrewer_1.1-2     checkmate_1.9.1       
[16] colorspace_1.3-2       htmltools_0.3.6        Matrix_1.2-12          plyr_1.8.4             pkgconfig_2.0.1       
[21] XML_3.98-1.19          genefilter_1.60.0      zlibbioc_1.24.0        xtable_1.8-3           scales_1.0.0          
[26] BiocParallel_1.12.0    htmlTable_1.13.1       tibble_1.4.2           annotate_1.56.2        mgcv_1.8-23           
[31] withr_2.1.2            nnet_7.3-12            lazyeval_0.2.1         survival_2.41-3        magrittr_1.5          
[36] crayon_1.3.4           memoise_1.1.0          nlme_3.1-131           MASS_7.3-48            foreign_0.8-69        
[41] vegan_2.5-3            tools_3.4.3            data.table_1.11.8      stringr_1.3.1          munsell_0.5.0         
[46] locfit_1.5-9.1         cluster_2.0.6          AnnotationDbi_1.40.0   ade4_1.7-13            compiler_3.4.3        
[51] rlang_0.3.0.1          rhdf5_2.22.0           grid_3.4.3             RCurl_1.95-4.11        biomformat_1.6.0      
[56] iterators_1.0.10       rstudioapi_0.8         htmlwidgets_1.3        igraph_1.2.2           bitops_1.0-6          
[61] base64enc_0.1-3        multtest_2.34.0        gtable_0.2.0           codetools_0.2-15       DBI_0.7               
[66] reshape2_1.4.3         gridExtra_2.3          knitr_1.22             bit_1.1-12             Hmisc_4.2-0           
[71] permute_0.9-4          ape_5.2                stringi_1.2.4          Rcpp_1.0.0             geneplotter_1.56.0    
[76] rpart_4.1-12           acepack_1.4.1          xfun_0.6

Thank you in advance! Alice

deseq2 • 64 views
ADD COMMENTlink modified 7 days ago by Michael Love23k • written 7 days ago by alice.checcucci0
Answer: RNAseq analysis DeSeq2 strange MAplot
gravatar for Michael Love
7 days ago by
Michael Love23k
United States
Michael Love23k wrote:

Can you describe your samples? How many are there? Is the sequencing depth even across samples?

ADD COMMENTlink written 7 days ago by Michael Love23k

The samples are 6 (3 for the wild type and 3 for the mutant). The sequencing depth is even across samples.

> colSums(counts(ddsTxi))
[1] 5030644 3461589 5517470 5246512 1681959 4651938
> cbind(samples, colSums(counts(ddsTxi)))
  species colSums(counts(ddsTxi))
1    2011                 5030644
2    2011                 3461589
3    2011                 5517470
4   delta                 5246512
5   delta                 1681959
6   delta                 4651938
ADD REPLYlink written 7 days ago by alice.checcucci0

I'm not sure what's going on. It looks like there are global changes to gene expression distribution, and I don't trust the in silico normalization approaches here. Can you show the boxplots for log normalized counts?

ADD REPLYlink written 5 days ago by Michael Love23k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 244 users visited in the last hour