Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • DESeq and edgeR up/down regulation

    Hi,

    I' am relatively new to bioinformatics.

    I' am assisting in analyzing 7 RNA-seq samples. Two conditions: (1) experimental - 6 samples, and (2) control - 1 sample.

    The input for DESeq consisted of read counts. See below for the R code for DESeq:

    ---------------------------------------------------------------
    library(DESeq)
    readCounts<-read.csv("readCounts_IndividualDigEx",row.names=1,header=TRUE)

    #which data are from MPSNT and Control
    group<-c(rep("MPNST",4),"Control",rep("MPNST",2))

    #annotate the data
    cds <- newCountDataSet(readCounts,group)

    #calculate library sizes and associated affect on DE testing
    cds <- estimateSizeFactors( cds )

    cds <- estimateVarianceFunctions( cds )

    #test for DE
    res <- nbinomTest( cds, "MPNST", "Control" )

    #there is no expression for a gene accross all samples... need to remove those genes
    whichIsNA<-which(is.na(res[,8]))
    res<-res[-whichIsNA,]

    #get list of DE genes, then display some of up and down-regulated genes
    resSig <- res[ res$pval < .05, ]
    #up-regulated
    head(resSig[order(resSig$foldChange,-resSig$baseMean),])
    #down-regulated
    head(resSig[order(-resSig$foldChange,-resSig$baseMean),])
    ---------------------------------------------------------------

    See the following for example data:

    a3 a4 a5 a6 a7 a8 a9
    MKL2 1521 991 1058 286 1965 1376 1026
    OOSP1 0 0 1 3 0 1 0
    CHL1 12 11 14 5154 3565 20 4


    My problem is the following. The list of DE genes severely under represents up-regulated genes (low expression in control, high in experiment). For example, with a pval<0.05 (FDR corrected): 215 genes were DE, 211 were down-regulated and 4 were up regulated.

    The issue does not seem related to the biology of the samples, but perhaps an error in setting up the tool. For example, I added two mock genes to test for DE:

    a3 a4 a5 a6 a7 a8 a9
    test1 10000 10000 10000 10000 1000 10000 10000
    test2 5000 5000 5000 5000 1000 5000 5000

    The mock genes were concatenated onto the other genes. The following are the p-values for each mock gene: test1: 0.25, test2: 0.30, and the adjusted p-values for both are 1.

    I understand that the tool confers much less confidence to expression measurements from conditions in which there is only one replicate, but the above outcome seems rather incorrect.

    Moreover, here is a sample gene which DESeq found to be significantly up-regluated (p-value = 0.0016):

    a3 a4 a5 a6 a7 a8 a9
    PAX7 4346 1 0 19 0 773 86


    Any advice would be greatly appreciated!

    Charlie
    Last edited by murphycj; 09-19-2011, 01:18 PM.

  • #2
    Do you mean you used RPKM values as input for DEseq? If so, you need to redo that analysis using raw count data. Tools like DEseq and DEGseq are expecting to have raw or simple count data, not normalized data.

    P.S. I use Bedtools' coverageBed function to get my raw counts from my mapping run BAM files (usually using a UCSC RefGene BED file as the reference).
    Last edited by mbblack; 09-19-2011, 12:07 PM.
    Michael Black, Ph.D.
    ScitoVation LLC. RTP, N.C.

    Comment


    • #3
      My apologies, the input data for DESeq were the read counts. I incorrectly said RPKM.

      Comment


      • #4
        You might have a problem with normalization. For each of your treatment library, plot the libraries individual fold w.r.t. the control library against the mean and compare with the size factors:
        Code:
        plot( ( counts(cds)[,1] + counts(cds)[,5] )/2, counts(cds)[,5] / counts(cds)[,1], log="xy", pch="." )
        abline( h = sizeFactors(cds)[5] / sizeFactors(cds)[1] )
        The horizontal line should go through the middle of the bulk of the genes, especially at the right-hand side of the plot. If it does not, the 'estimateSizefactor' function has not worked well. Maybe try this and let us know.

        Comment


        • #5
          I did what you suggested, and the plot seems to suggest that the 'estimateSizeFactors' works correctly. The horizontal line goes through the middle of the bulk of the genes. The plot looks very similar across all library comparisons (e.g. 1v5, 2v5, 3v5, etc..).

          See the attached .pdf for an example.

          Thanks for your reply.
          Attached Files
          Last edited by murphycj; 09-20-2011, 05:05 PM.

          Comment


          • #6
            Maybe the following two plots will help.

            (1) I have attached a histogram of the p-values. The distribution seems to indicate an error in the statistical test.

            (2) The MvA plot.

            Hope they help.

            Additionally, the mean fold change is -0.169 with a standard deviation of 1.5.
            Attached Files
            Last edited by murphycj; 09-20-2011, 06:20 PM.

            Comment


            • #7
              This seems all fine to me.

              The reason for the strong asymmetry in detection power depending on the direction of change is due to your extremely unbalanced design. To get the same precission of expression strength estimation, a gene has to be much stronger expressed in control than in treatment. Hence, if a gene is low in control and high in treatment, the fold change cannot be estimated as well than if it the other way around. DESeq takes this into account and so is more reluctant to call down-regulated genes significant.

              Comment


              • #8
                Seems you are right. To test if the asymmetry is due to the unbalanced design, I replicated the the control expression data set (i.e. I inputted into DESeq 8 samples: 6 experimental and 2 control). I of course randomly modified the second control expression values by multiplication by a random number, which had a mean of 1 and low standard deviation.

                What results is many fold more DE genes (due to the similarity of the two control data sets), but more importantly, roughly equal amounts of up and down-regulated genes.

                Thanks for your input on this issue. The reason for the biased DESeq results make sense at least.

                Our collaborators stated that they are in the process of adding more controls to the data.

                Comment

                Latest Articles

                Collapse

                • seqadmin
                  Current Approaches to Protein Sequencing
                  by seqadmin


                  Proteins are often described as the workhorses of the cell, and identifying their sequences is key to understanding their role in biological processes and disease. Currently, the most common technique used to determine protein sequences is mass spectrometry. While still a valuable tool, mass spectrometry faces several limitations and requires a highly experienced scientist familiar with the equipment to operate it. Additionally, other proteomic methods, like affinity assays, are constrained...
                  04-04-2024, 04:25 PM
                • 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

                ad_right_rmr

                Collapse

                News

                Collapse

                Topics Statistics Last Post
                Started by seqadmin, 04-11-2024, 12:08 PM
                0 responses
                18 views
                0 likes
                Last Post seqadmin  
                Started by seqadmin, 04-10-2024, 10:19 PM
                0 responses
                22 views
                0 likes
                Last Post seqadmin  
                Started by seqadmin, 04-10-2024, 09:21 AM
                0 responses
                17 views
                0 likes
                Last Post seqadmin  
                Started by seqadmin, 04-04-2024, 09:00 AM
                0 responses
                49 views
                0 likes
                Last Post seqadmin  
                Working...
                X