Seqanswers Leaderboard Ad



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

  • DESeq2 - factor base level changes the DE genes


    I have an RNAseq experiment from two strains of mice, that have been exposed to two different environments (three biological replicates of each combination). So my experimental setup is

    > colData
    sample genetics environment
    black2 B6 black
    black3 B6 black
    black5 B6 black
    black1 B6 yellow
    black4 B6 yellow
    black6 B6 yellow
    agouti1 129 yellow
    agouti4 129 yellow
    agouti6 129 yellow
    agouti2 129 black
    agouti3 129 black
    agouti5 129 black

    I want to test which genes change their expression when the environment changes, controlling for the genetic background. I am using DESeq2 for this, with the formula ~ genetics + environment + genetics:environment.

    However, the choice of the base level for the factors in my experimental design changes the genes that are significant. If I leave the default reference levels (129 for genetics and black for env) I get more DE genes that if I set as reference B6 for genetics and yellow for environemnt.
    I thought the reference level is used to decide how to calculate the fold change and is necessary to correctly interpret the coefficients, but shouldn't the results be the same?

    Any thoughts on what am I missing?

    For cases where you have a control and mutant, it makes sense to use the control as the base level and compare everything to that, but in this case, I don't have any valid reason to chose one strain over the other. My results are different depending on this, so what is the correct form of doing this sort of analysis?

    Thank you very much in advance.

    My code is below.

    ## 1 ##
    > levels(colData$genetics)
    [1] "129" "B6"
    > levels(colData$environment)
    [1] "black" "yellow"

    > dds <- DESeqDataSetFromMatrix(countData = data[,8:19], colData = colData, design = ~ genetics + environment + genetics:environment)
    > dds <- DESeq(dds)

    > resG <- results(dds, alpha=0.05, name="genetics_B6_vs_129")
    > resE <- results(dds, alpha=0.05, name="environment_yellow_vs_black")
    > resGE <- results(dds, alpha=0.05, name="geneticsB6.environmentyellow")

    > dim(subset(resG, resG$padj < 0.05))
    [1] 3675 6
    > dim(subset(resE, resE$padj < 0.05))
    [1] 20 6
    > dim(subset(resGE, resGE$padj < 0.05))
    [1] 24 6

    ## 2 ##
    > colData2 <- colData
    > colData2$genetics <- relevel(colData2$genetics, ref="B6")
    > colData2$environment <- relevel(colData2$environment, ref="yellow")
    > levels(colData2$genetics)
    [1] "B6" "129"
    > levels(colData2$environment)
    [1] "yellow" "black"

    > dds2 <- DESeqDataSetFromMatrix(countData = data[,8:19], colData = colData2, design = ~ genetics + environment + genetics:environment)
    > dds2 <- DESeq(dds2)

    > res2G <- results(dds2, alpha=0.05, name="genetics_129_vs_B6")
    > res2E <- results(dds2, alpha=0.05, name="environment_black_vs_yellow")
    > res2GE <- results(dds2, alpha=0.05, name="genetics129.environmentblack")

    > dim(subset(res2G, resG$padj < 0.05))
    [1] 3127 6
    > dim(subset(res2E, resE$padj < 0.05))
    [1] 2 6
    > dim(subset(res2GE, resGE$padj < 0.05))
    [1] 0 6


    > g1 <- subset(resG, resG$padj < 0.05)
    > g2 <- subset(res2G, res2G$padj < 0.05)

    > length(intersect(rownames(g1), rownames(g2)))
    [1] 2640

    -- output of sessionInfo():

    R version 3.0.2 (2013-09-25)
    Platform: x86_64-apple-darwin10.8.0 (64-bit)

    [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

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

    other attached packages:
    [1] DESeq2_1.2.10 RcppArmadillo_0.4.320.0 Rcpp_0.11.2
    [4] GenomicRanges_1.14.4 XVector_0.2.0 IRanges_1.20.7
    [7] BiocGenerics_0.8.0

    loaded via a namespace (and not attached):
    [1] annotate_1.40.1 AnnotationDbi_1.24.0 Biobase_2.22.0 DBI_0.2-7
    [5] genefilter_1.44.0 grid_3.0.2 lattice_0.20-23 locfit_1.5-9.1
    [9] RColorBrewer_1.0-5 RSQLite_0.11.4 splines_3.0.2 stats4_3.0.2
    [13] survival_2.37-7 tools_3.0.2 XML_3.95-0.2 xtable_1.7-3

  • #2
    At a guess, I'd say it has to be your interaction term. When you change the parameters, you change the interaction term, and the fact that the results are different is telling you there is some difference in the interaction of 129&black versus B6&yellow.

    DESeq2 includes some mechanisms of looking at those specific interaction affects if I remember the vingette correctly? So you might want to read through that and try to pull out the interaction affect in both your contrasts. I think it is potentially telling you something about your actual data and biological responses, so you need to determine if that interaction term is significant or not. You should also try the other, alternate combinations of interaction so you get a complete overview of which interactions are or are not significant. Significant interaction affects will really confound trying to determine the significance of individual factors if you wish to stick with a single combined data set and analysis. You may have to separate the data so you can analyze a single factor at a time without any interaction.

    Or again if the interaction is significant, you could look at some alternate normalizations to subtract one factor before significance testing. For example a TMM or KDMM normalization with the 129 mean genotype responses subtracted as baseline.
    Last edited by mbblack; 08-13-2014, 12:06 PM.
    Michael Black, Ph.D.
    ScitoVation LLC. RTP, N.C.


    • #3
      hi xibarra,

      Note that version 1.2 is from October 2013. It's admittedly a hassle to keep R and Bioc updated, but version 1.2 was still only the second release of DESeq2. Some users reported to us last year that for certain designs, changing the base level could change the LFCs and p-values, usually only slightly, due to an asymmetry of the prior on LFCs. We then developed the methods such that changes of base level would have no effect on the results, and these were released with version 1.4 in March 2014. An additional change was to not shrink the main effect terms in designs with interaction terms.

      So you can either update R and Bioc, or for the old version (1.2) you can set betaPrior=FALSE, which should produce symmetric results tables.


      • #4
        Thank you for your replies.
        I suspected I needed to update to the latest version, which I'll do. Still nice to know I'm not misunderstanding something or doing the wrong analysis.


        Latest Articles


        • 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





        Topics Statistics Last Post
        Started by seqadmin, 07-19-2024, 07:20 AM
        0 responses
        Last Post seqadmin  
        Started by seqadmin, 07-16-2024, 05:49 AM
        0 responses
        Last Post seqadmin  
        Started by seqadmin, 07-15-2024, 06:53 AM
        0 responses
        Last Post seqadmin  
        Started by seqadmin, 07-10-2024, 07:30 AM
        0 responses
        Last Post seqadmin