Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • DESeq Normalization Question

    I'm sure this issue has come up before, but I couldn't find an appropriate thread or answer either here or on the Bioconductor mailing list.

    What feature of the data or the distribution of counts among my samples can cause the sizeFactors to vary much more than the raw counts / library sizes?

    More detail: I'm using DESeq to analyze RNA-seq data mapped with STAR, counted with htseq-count. Comparing the "doubleTerm" samples to the "wt" samples, there are many genes that appear downregulated. While these samples were sequenced, on average, to a similar sequencing depth, the normalization factors are much smaller for WT, resulting in much larger normalized counts, resulting in more apparently downregulated genes in doubleTerm vs WT.
    Code:
    > cds <- newCountDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory)
    > cds <- estimateSizeFactors(cds)
    > cds <- estimateDispersions(cds)
    > data.frame(sizefactors=sizeFactors(cds), rawcounts=colSums(counts(cds, normalized=FALSE)))
                    sizefactors rawcounts
    S01_wt1           0.9016089  23466349
    S02_wt2           0.7679168  22428603
    S03_wt3           0.7952564  19841959
    S04_wt4           0.7839629  18363384
    S05_pten8w1       1.0301769  20859853
    S06_pten8w2       0.9949514  16809588
    S07_pten8w3       0.9425865  16731071
    S08_pten22w1      1.0826846  18906329
    S09_pten22w2      1.1640354  20164026
    S10_pten22w3      1.0111748  17306468
    S11_double8w1     0.7949001  17671986
    S12_double8w2     1.4509978  23673557
    S13_double8w3     1.1703853  22127841
    S14_doubleterm2   1.0786455  19063010
    S15_doubleterm4   1.1265935  19279814
    S16_doubleterm6   1.3059472  22750403
    Thanks.

    Stephen

    Code:
    > sessionInfo()
    R version 3.0.0 (2013-04-03)
    Platform: x86_64-apple-darwin10.8.0 (64-bit)
    
    locale:
    [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
    
    attached base packages:
    [1] parallel  stats     graphics  grDevices utils     datasets
    methods   base
    
    other attached packages:
    [1] DESeq_1.12.0         lattice_0.20-15      locfit_1.5-9
    Biobase_2.20.0
    [5] BiocGenerics_0.6.0   edgeR_3.2.3          limma_3.16.2
    BiocInstaller_1.10.1
    
    loaded via a namespace (and not attached):
     [1] annotate_1.38.0      AnnotationDbi_1.22.3 DBI_0.2-6
    DESeq2_1.0.9
     [5] genefilter_1.42.0    geneplotter_1.38.0   GenomicRanges_1.12.2
    grid_3.0.0
     [9] IRanges_1.18.0       RColorBrewer_1.0-5   RSQLite_0.11.3
    splines_3.0.0
    [13] stats4_3.0.0         survival_2.37-4      tools_3.0.0
    XML_3.95-0.2
    [17] xtable_1.7-1

  • #2
    One thing that could cause this is extreme overrepresentation of one or a handful of genes, like hemoglobin in whole-blood samples for instance. Or maybe rubisco in plants. Partial/failed rRNA removal? I would look for something that seems to dominate the RNA pool.

    Comment


    • #3
      Dear Stephen,

      how does the 'pairs' (or LSD::heatpairs) plot look like, or the pairwise MA plots? You could also try the arrayQualityMetrics report. From this, problems such as suggested by kopi-o might become apparent.

      Best wishes
      Wolfgang.
      Wolfgang Huber
      EMBL

      Comment


      • #4
        I have what is probably a naive question about DESeq normalization. The manual says that by dividing the count by the size factor, one makes samples comparable. Does this include comparable for, say, eQTL analysis? I randomly pulled some genes from counts(cds,normalized=TRUE) and plotted their distributions, and they are not all normalliy distributed. Some look normal-ish (bottom row), but I still don't know if I'd consider them normal; others are very clearly following other distributions (top row). If I want to identify eQTLs from my RNAseq data (as I also have genotype data on those individuals), then all the gene read counts need to be transformed to follow a normal distribution before I can test for eQTLs.

        So my question is, is the command counts(cds,normalized=TRUE) designed to transform the raw reads counts into expression levels that follow a normal distribution? If so, why do some of my genes not look normally distributed? If not, how could I transform my countDataSet so that each gene follows a normal distribution? If my dataset consists of expression data from three time points, could that be messing up what would otherwise be a normally distributed gene? Figure of some distributions for genes from counts(cds,normalized=TRUE) are below. Thanks for any suggestions!

        Comment


        • #5
          I would agree with kopi-o. The normalisation isn't based on the total counts, it's based on the distribution of counts across all genes. If you want a slightly more accurate comparison to size factors (but fairly easy to calculate), try splitting the count data into quantiles, or look at (say) the 75th percentile count instead of the total count.

          Comment


          • #6
            I came across this lovely thing: http://www.bios.unc.edu/research/gen..._eQTL/faq.html

            "Outliers in expression data are usually harder to deal with. The accepted remedy by the GTEx consortium is the transformation of the measurements for each gene into normally distributed while preserving relative rankings. The target distribution may be the standard normal distribution or the normal distribution the mean and spread of the original measurements. Here is the code for such transformation:"

            for( sl in 1:length(gene) ) {
            mat = gene[[sl]];
            mat = t(apply(mat, 1, rank, ties.method = "average"));
            mat = qnorm(mat / (ncol(gene)+1));
            gene[[sl]] = mat;
            }
            rm(sl, mat);

            I used my normalized DESeq count data as input, then used the program to transform each gene to a normal distribution of expression. Comparing before and after transformation for a few genes, they certainly look normal.



            The program claims to have been used successfully to identify eQTLs in RNAseq data. Whether or not my using this approach for eQTLs turns out to be biologically relevant, informative, or correct is yet to be determined. Has anyone tried other transformations of RNAseq for eQTL analysis?

            Comment

            Latest Articles

            Collapse

            • seqadmin
              Strategies for Sequencing Challenging Samples
              by seqadmin


              Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
              03-22-2024, 06:39 AM
            • seqadmin
              Techniques and Challenges in Conservation Genomics
              by seqadmin



              The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

              Avian Conservation
              Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
              03-08-2024, 10:41 AM

            ad_right_rmr

            Collapse

            News

            Collapse

            Topics Statistics Last Post
            Started by seqadmin, Yesterday, 06:37 PM
            0 responses
            7 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, Yesterday, 06:07 PM
            0 responses
            7 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 03-22-2024, 10:03 AM
            0 responses
            49 views
            0 likes
            Last Post seqadmin  
            Started by seqadmin, 03-21-2024, 07:32 AM
            0 responses
            66 views
            0 likes
            Last Post seqadmin  
            Working...
            X