Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • DESeq2 without biol replicates

    Hey all,

    i wanted to try DESeq2 in a pairwise comparison setting without replicates.

    unfortunately using the code below on a pasilla-test-case i can only find differentially expressed genes using DESeq and not DESeq2.
    (DESeq2 says "same number of samples and coefficients to fit, estimating dispersion by treating samples as replicates" so should be getting the task)

    i know it must be me, but where am I going wrong?

    thanks a lot!!


    #-------------------------- deseq2

    Code:
    library(DESeq2)
    library(pasilla)
    data("pasillaGenes")
    countData <- counts(pasillaGenes)
    countData<-countData[,c("treated1fb","untreated1fb")]
    colData <- pData(pasillaGenes)[c("treated1fb","untreated1fb"),c("condition","type")]
    dds <- DESeqDataSetFromMatrix(
           countData = countData,
           colData = colData,
           design = ~ condition)
    colData(dds)$condition <- factor(colData(dds)$condition,levels=c("treated","untreated"))
    dds
    dds <- DESeq(dds)
    res <- results(dds)
    res <- res[order(res$padj),]
    head(res)
    DataFrame with 6 rows and 5 columns
    baseMean log2FoldChange lfcSE pvalue padj
    <numeric> <numeric> <numeric> <numeric> <numeric>
    FBgn0000008 60.7056967 -0.09993564 0.36063169 0.7816935 0.9973207
    FBgn0000014 0.8344279 -0.03920709 0.08405437 0.6408940 0.9973207
    FBgn0000015 0.4172140 -0.02034815 0.06095374 0.7385084 0.9973207

    #-------------------------- deseq

    Code:
    library(DESeq)
    library(pasilla)
    data("pasillaGenes")
    countData <- counts(pasillaGenes)
    countData<-countData[,c("treated1fb","untreated1fb")]
    condition <- factor( c( "treated", "untreated") )
    cds <- newCountDataSet( countData, condition )
    cds <- estimateSizeFactors( cds )
    cds <- estimateDispersions( cds, method="blind", sharingMode="fit-only" )
    res <- nbinomTest( cds, "treated", "untreated") 
    res <- res[order(res$padj),]
    head(res)
    id baseMean baseMeanA baseMeanB foldChange log2FoldChange
    9696 FBgn0039155 513.8001 31.70826 995.89187 31.40796187 4.973058
    2605 FBgn0026562 20645.8226 6151.40275 35140.24243 5.71255758 2.514137
    29 FBgn0000071 241.5499 444.75009 38.34963 0.08622736 -3.535710
    pval padj
    9696 9.669445e-23 1.062188e-18
    2605 1.477349e-12 8.114339e-09
    29 3.097019e-10 1.134025e-06

  • #2
    To be honest, we couldn't yet be bothered to explain how to analyse such data in DESeq2. It's tricky to write up because too many people will misinterpret whatever I write as if it were actually possible to conduct a meaningful statistical analysis when comparing just two samples.

    So, if you promise to not use any such comparisons for actual science, here is how you do it:

    Start as above:

    Code:
    library(DESeq2)
    library(pasilla)
    data("pasillaGenes")
    countData <- counts(pasillaGenes)
    countData<-countData[,c("treated1fb","untreated1fb")]
    colData <- pData(pasillaGenes)[c("treated1fb","untreated1fb"),c("condition","type")]
    dds <- DESeqDataSetFromMatrix(
           countData = countData,
           colData = colData,
           design = ~ condition)
    Now use DESeq2's new "rlog transformation". This replaced the VST we had before. It transforms the average of the genes across samples to a log2 scale but "pulls in" those genes for which the evidence for strong fold changes is weak due to low counts.

    Code:
    rld <- rlogTransformation( dds )
    As this is a logarithm-like scale, the differences between the samples can be considered as a "regularized log2 fold change". Let's make a result data frame:
    Code:
    res <- data.frame(
       assay(rld), 
       avgLogExpr = ( assay(rld)[,2] + assay(rld)[,1] ) / 2,
       rLogFC = assay(rld)[,2] - assay(rld)[,1] )
    And now we have a ranking of genes by regularized fold changes:

    Code:
    > head( res[ order(res$rLogFC), ] )
                treated1fb untreated1fb avgLogExpr    rLogFC
    FBgn0011260   7.830359     6.627326   7.228842 -1.203033
    FBgn0001226  10.128636     8.929985   9.529311 -1.198652
    FBgn0034718   8.503006     7.315640   7.909323 -1.187366
    FBgn0003501   7.927864     6.743974   7.335919 -1.183889
    FBgn0033635  11.126300     9.973979  10.550139 -1.152321
    FBgn0033367  13.411814    12.269436  12.840625 -1.142378
    This ranking put ones genes on top which are strongly downregulated in the second sample compared to the first one. If you do this with normal log expressions, the weak genes will appear at the top because they are noisiest and hence tend to have exaggerated fold changes.

    The advantage of this procedure is that it does not produce any p values (which would be misleading anyway).
    Last edited by Simon Anders; 06-12-2013, 05:10 AM.

    Comment


    • #3
      thanks for your help Simon- you are the best!

      > So, if you promise to not use any such comparisons for actual science
      hihihi...
      you are totally right of course!

      Comment


      • #4
        hey,

        could you please help again, i am a bit as it seems that the DESeq and DESeq2 regularized LogFC calculations differ strongly if i stick to the recommended procedure

        looked for a potentially differentially expressed gene from pasilla (FBgn0039155 counts are 38 to 831), wanted to see what the DESeq2 regularized LogFC would be and tried to compared it to the DESeq rLogFC as depicted below.
        DESeq says rLogFC is 3.095034, DESeq2 says 0.1595547 (and the later doesn't change much if i try out blind rlogTransformation, estimateSizeFactors, estimateDispersions). somehow i believe the DESeq2 procedure i used isn't doing what i want.

        thanks!

        ----------------- deseq2

        Code:
        library(DESeq2)
        library(pasilla)
        data("pasillaGenes")
        countData <- counts(pasillaGenes)
        countData<-countData[,c("treated1fb","untreated1fb")]
        head(countData)
        colData <- pData(pasillaGenes)[c("treated1fb","untreated1fb"),c("condition","type")]
        dds <- DESeqDataSetFromMatrix(
               countData = countData,
               colData = colData,
               design = ~ condition)
        rld <- rlogTransformation( dds)#,blind = TRUE)
        res <- data.frame(
           assay(rld), 
           avgLogExpr = ( assay(rld)[,2] + assay(rld)[,1] ) / 2,
           rLogFC = assay(rld)[,2] - assay(rld)[,1] )
        res[c("FBgn0039155"),]
        treated1fb untreated1fb avgLogExpr rLogFC
        FBgn0039155 8.850547 9.010102 8.930324 0.1595547

        ----------------- deseq

        Code:
        library(DESeq)
        library(pasilla)
        data("pasillaGenes")
        countData <- counts(pasillaGenes)
        countData<-countData[,c("treated1fb","untreated1fb")]
        condition <- factor( c( "treated", "untreated") )
        cds <- newCountDataSet( countData, condition )
        cds <- estimateSizeFactors( cds )
        cds <- estimateDispersions( cds, method="blind", sharingMode="fit-only" )
        vsd <- varianceStabilizingTransformation( cds )
        res <- nbinomTest( cds, "treated", "untreated") 
        rlfc <- (exprs(vsd)[, c("untreated1fb")]  -exprs(vsd)[, c("treated1fb")] )
        res$rLogFC<-rlfc
        res[res$id=="FBgn0039155",]
        id baseMean baseMeanA baseMeanB foldChange log2FoldChange
        9696 FBgn0039155 513.8001 31.70826 995.8919 31.40796 4.973058
        pval padj rLogFC
        9696 9.669445e-23 1.062188e-18 3.095034
        Last edited by sisterdot; 06-12-2013, 11:10 PM.

        Comment


        • #5
          hi,

          I made a change to rlogTransformation in the development branch of DESeq2 (v 1.1.21), which might be relevant for your analysis.

          Previously, the rlog transformation by default would use a blind design formula and use the MAP dispersion estimates (the blue points in a plot produced by plotDispEsts). True, large log fold changes would have high MAP dispersions and be shrunk toward 0 (especially with no replicates). Instead we now (in the devel branch) use the fitted dispersion estimates, which shrinks log fold changes based on the general mean-dispersion trend line, and preserves true large log fold changes.

          It is described here:

          Comment


          • #6
            Hi Michael,
            Is the shrunken log fold change also applied in DESeq2 v 1.0.17?

            Comment


            • #7
              Yes, shrunken log fold changes were incorporated for the Wald test since the beginning of DESeq2.

              You can always check the methods of the current version with ?DESeq, vignette("DESeq2"), or for more specific updates: news(package="DESeq2").

              Comment


              • #8
                You are comparing DESeq(1) and DESeq2. DESeq(1) does not have and never had any shrinkage on fold changes.

                Comment


                • #9
                  Is there an implicit normalization across the samples in this?

                  Comment


                  • #10
                    There are loads of threads on here about DeSeq and experiments with no replicates. And Simon Anders has been very patient coming on here and repeatedly explaining why they are a bad idea.

                    So I have come along after an experiment has been done to have a second look at the data (some bioinformatics already done by the sequencing provider) and.. yes.. they didn't have any replicates!

                    The bioinformatics guys that came before used DeSeq and they DID find differentially expressed transcripts. Even with adjusted p-values.

                    So just to be clear - I should be telling these guys they can't trust the p-values anyway since they have no replicates and to just look at fold change? We're going for just taking the top 200 absaolute fold change (for now) and confirming with qPCR.

                    Comment


                    • #11
                      The section of the original DESeq paper might shed some light:


                      "Working without replicates
                      DESeq allows analysis of experiments with no biological replicates in one or even both of the conditions. While one may not want to draw strong conclusions from such an analysis, it may still be useful for exploration and hypothesis generation. If replicates are available only for one of the conditions, one might choose to assume that the variance-mean dependence estimated from the data for that condition holds as well for the unreplicated one. If neither condition has replicates, one can still perform an analysis based on the assumption that for most genes, there is no true differential expression, and that a valid mean-variance relationship can be estimated from treating the two samples as if they were replicates. A minority of differentially abundant genes will act as outliers; however, they will not have a severe impact on the gamma-family GLM fit, as the gamma distribution for low values of the shape parameter has a heavy right-hand tail. Some overestimation of the variance may be expected, which will make that approach conservative."

                      Comment


                      • #12
                        Originally posted by whataBamBam View Post
                        There are loads of threads on here about DeSeq and experiments with no replicates. And Simon Anders has been very patient coming on here and repeatedly explaining why they are a bad idea.

                        So I have come along after an experiment has been done to have a second look at the data (some bioinformatics already done by the sequencing provider) and.. yes.. they didn't have any replicates!

                        The bioinformatics guys that came before used DeSeq and they DID find differentially expressed transcripts. Even with adjusted p-values.

                        So just to be clear - I should be telling these guys they can't trust the p-values anyway since they have no replicates and to just look at fold change? We're going for just taking the top 200 absaolute fold change (for now) and confirming with qPCR.
                        I think it is fairly obvious that Simon has a particular interpretation or hypothesis in mind when he is making statements about p-values(IE for the hypothesis that it is possible to test with DEseq), which isn't necessarily "valid science" either. A p-value is just the probability of the observations to be as extreme or more extreme under the null hypothesis, but you have to be aware of the hypothesis that was tested(Were there genes expressed at a different level in one sample and not the other? Could be answered quite scientifically.). Though it may be suspect that DEseq would return significant p-values, when the p-values couldn't possibly be significant, is it really answering the question that we think it is testing?

                        Comment


                        • #13
                          I wouldn't say that it is suspect to return low p-values. The point of using most software specifically designed for microarray or sequencing data, rather than just running linear or generalized linear models row-by-row, is that information can be pooled across genes.

                          for example, if you observe two random Normal variables and want to know if they have different means, you're out of luck. But if I provide you the population variance then you can compute a probability of seeing such a difference. Here we have an in-between case, we can learn something about the variance over all genes, but are not provided the population variance for a given gene. The question is how much can we learn about the variance of a single gene by looking over all genes. Are most of the genes not differentially expressed, such that one can learn about the mean-variance relationship over the majority of genes? The paragraph I quote above points out that the behavior will be conservative, in overestimating the variance of all genes by including some number of differentially expressed genes.

                          Comment


                          • #14
                            Hello,

                            I deviate a bit from the discussion but my question fits well to this topic.
                            I am analyzing RNASeq data of 2 conditions and 7 replicates per condition. I performed the analysis with DESeq2, this is ok.

                            In one condition, two samples come from the same patient with a leukemia. At the second time, the patient is in relapse.
                            I would like to find what changed between these 2 time points. Of course, it can be due to additional somatic mutations. I want to study the differences on RNA as well.

                            I think about performing a DE analysis between these 2 "conditions" (remission, relapse). Do you think it is meaningful to perform this analysis, knowing there is only one patient?

                            How are you use to handle these cases? It is very patient-dependent I guess, I don't know if we could see general changes when using small groups of patients in remission and relapse.

                            Comment


                            • #15
                              Originally posted by Michael Love View Post
                              The section of the original DESeq paper might shed some light:


                              "Working without replicates
                              DESeq allows analysis of experiments with no biological replicates in one or even both of the conditions. While one may not want to draw strong conclusions from such an analysis, it may still be useful for exploration and hypothesis generation. If replicates are available only for one of the conditions, one might choose to assume that the variance-mean dependence estimated from the data for that condition holds as well for the unreplicated one. If neither condition has replicates, one can still perform an analysis based on the assumption that for most genes, there is no true differential expression, and that a valid mean-variance relationship can be estimated from treating the two samples as if they were replicates. A minority of differentially abundant genes will act as outliers; however, they will not have a severe impact on the gamma-family GLM fit, as the gamma distribution for low values of the shape parameter has a heavy right-hand tail. Some overestimation of the variance may be expected, which will make that approach conservative."
                              Great. Actually my original interpretation (before I posted this) was correct then. That the p values are perfectly valid (in fact conservative) and the problem of no replicates is actually low statistical power.

                              So basically you are saying that you have less statistical power because you have overestimated the variance. And if you see significant differences DESPITE this low statistical power then go for it.

                              To be fair it says in the vignette (or the paper I can't remember which) that there is simply low statistical power if you have no replicates.

                              Comment

                              Latest Articles

                              Collapse

                              • seqadmin
                                Exploring the Dynamics of the Tumor Microenvironment
                                by seqadmin




                                The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
                                07-08-2024, 03:19 PM
                              • seqadmin
                                Exploring Human Diversity Through Large-Scale Omics
                                by seqadmin


                                In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
                                06-25-2024, 06:43 AM

                              ad_right_rmr

                              Collapse

                              News

                              Collapse

                              Topics Statistics Last Post
                              Started by seqadmin, Yesterday, 07:20 AM
                              0 responses
                              23 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 07-16-2024, 05:49 AM
                              0 responses
                              38 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 07-15-2024, 06:53 AM
                              0 responses
                              44 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 07-10-2024, 07:30 AM
                              0 responses
                              41 views
                              0 likes
                              Last Post seqadmin  
                              Working...
                              X