Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Significant up-regulated and downregulated genes with DESeq

    Hi,

    I'm trying to identify the significant up-regulated and down-regulated genes in 6 different moments of the life cycle of "my" animal.

    I have 11 libraries (2 replicates/stage, except for one stage that i couldn't sequence a replicate). What i'm doing is to compare one stage with the stage immediately after using DESeq.

    This is my script:


    Code:
    #Load data
    Cele_SPvsEB = read.csv (file.choose(), header=TRUE, row.names=1)
    CeleDesign <- data.frame(
      row.names = colnames(Cele_SPvsEB),
      condition = c("SP", "SP", "EB", "EB"))
    
    #to load conditions w/out pData file: conditions
    conds = factor(c("SP", "SP", "EB", "EB"))
    
    #Make cds file
    cds <-newCountDataSet(Cele_SPvsEB, conds )
    
    #Est size factor = normalize for library size
    cds<-estimateSizeFactors(cds) 
    
    #estimate dispersion from the expected
    cds <- estimateDispersions(cds, method=c("pooled"), fitType=c("local"))
    str(fitInfo(cds))
    
    #plot dispersion
    plotDispEsts <- function( cds )
    {
       plot(
          rowMeans( counts( cds, normalized=TRUE ) ),
          fitInfo(cds)$perGeneDispEsts,
          pch = '.', log="xy", ylab="dispersion", xlab="mean of normalized counts" )
       xg <- 10^seq( -.5, 5, length.out=300 )
       lines( xg, fitInfo(cds)$dispFun( xg ), col="red" )
    }
    plotDispEsts( cds )
    
    #call differential expression
    res <- nbinomTest( cds, "SP", "EB" )
    
    #there is no expression for a gene accross all samples... need to remove those genes
    whichIsNA<-which(is.na(res[,8]))
    res<-res[-whichIsNA,]
    
    #plot volcano
    plotDE <- function( res )
       plot(
          res$baseMean,
          res$log2FoldChange,
          log="x", pch=20, cex=.3,
          col = ifelse( res$pval < .05, "red", "black" ) )
    
    plotDE( res )
    
    #filter for upregulated and downregulated genes
    
    resSig <- res[ res$pval < .05, ]
    up <- subset(resSig, foldChange > 2)
    down <- subset(resSig, foldChange < 2)
    
    #list most sig genes
    head( resSig[ order(res$pval), ] )
    By using this script i found around 600 genes in both up-regulated and down-regulated categories for this example (i.e. SP vs EB).

    But, if i change my filtering to:

    Code:
    resSig <- res[ res$padj < .1, ]
    up <- subset(resSig, foldChange > 2)
    down <- subset(resSig, foldChange < 2)
    Then i have 0 genes up or down-regulated.

    If i just keep going only filtering by pval < 0.05, is that enough to consider my genes significantly differentially expressed?

    Do i have any error in my script?

    Any suggestions?

    Thanks!!!
    Last edited by alisrpp; 10-15-2013, 11:31 AM.

  • #2
    If i just keep going only filtering by pval < 0.05, is that enough to consider my genes significantly differentially expressed?
    No. No one will (or at least should) believe those are differentially expressed without additional evidence. Also, keep in mind also that that's probably log2(fold-change), so your up/down regulated may not mean what you think (I haven't used DESeq1 in a while, but I remember this being the case). So, things to think about trying are as follows:
    1. Try DESeq2, which has a number of very nice improvements.
    2. Use the genefilter package to help properly perform independent filtering. This can drastically up your power.


    These two options aren't mutually exclusive, of course.

    Comment


    • #3
      Can you calculate fold change from log2FC and the filter out above say 1.5?

      You set it at 2, why not lower? cant you set it at 1.5?

      Comment


      • #4
        Originally posted by sindrle View Post
        Can you calculate fold change from log2FC and the filter out above say 1.5?
        Yes, that works.

        You set it at 2, why not lower? cant you set it at 1.5?
        The thresholds are all semi-arbitrary, you could set them to whatever makes sense for your experiment.

        Comment


        • #5
          How do I know "it makes sense". Say you have a sign. diff. exp. gene at about 30 FPKM and a FC of 1.5. Maybe not much? But its biological function makes a lot of sense in the context.
          What then do you do?

          Comment


          • #6
            There's no rule I can give you for determining "what makes sense", it'll depend entirely on the system you're using, your end goals, the resources (financial and otherwise) at hand, and how risk averse you are. At some point you just have to become comfortable with contextualizing results in terms of biological implications (I only filter by adjusted p-value for this reason). Back in the early microarray days, using fold-change thresholds was pretty common, seemingly due to the methods not being that great. These days, I haven't seen as many people doing that.

            Comment


            • #7
              Yeah, that was what i thought!

              Question about "filtering", does it make sense only to keep the intersect of Cuffdiff, edger and deseq2? Or how to know which results to use?

              Comment


              • #8
                Presumably you're going to do some validations via qPCR or something. So just include a few unique hits from each in an initial screen to get an idea which is better fitting your experiment. I would expect either edgeR or DESeq2 to produce the best results, but it's easy enough to check.

                Comment


                • #9
                  Yes, I have calculated the result of one gene so far with qPCR. And edgeR is in the lead!

                  Make that 3/3 so far.
                  Last edited by sindrle; 10-29-2013, 09:19 AM.

                  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, 03-27-2024, 06:37 PM
                  0 responses
                  12 views
                  0 likes
                  Last Post seqadmin  
                  Started by seqadmin, 03-27-2024, 06:07 PM
                  0 responses
                  11 views
                  0 likes
                  Last Post seqadmin  
                  Started by seqadmin, 03-22-2024, 10:03 AM
                  0 responses
                  53 views
                  0 likes
                  Last Post seqadmin  
                  Started by seqadmin, 03-21-2024, 07:32 AM
                  0 responses
                  69 views
                  0 likes
                  Last Post seqadmin  
                  Working...
                  X