How to get multi-factor Deseq LRT test log P values for a particular condition
2
0
Entering edit mode
aishu.jp ▴ 30
@aishujp-15144
Last seen 5.2 years ago

Hi, My RNAseq dataset consist of two cultivars (R and S), three time points (1day, 3day and 7day) and two treatments (Control, Treated) with two replicates each. I have a total of 24 samples.

I am trying to find DEGs using LRT test for the following condition

> R Vs S due Treated Vs Control within 1 day
   R Vs S due Treated Vs Control within 3 day
   R Vs S due Treated Vs Control within 7 day

My code is 

> library(DESeq2)
   countdata <- read.table("Suscep_Resis_ctrl_trt_all_matrix.txt", header=TRUE, row.names=1)
   countdata
   countdata <- as.matrix(countdata)
   countdata
   colData<- read.table ("Metadata.txt", header=TRUE, row.names=1)
   colData
   dds <- DESeqDataSetFromMatrix(countData = countdata, colData = colData, design= ~ Genotype.S.R. +   Treatmenttime.C.T. +  Genotype..S.R.:Treatmenttime.C.T.)
    dds
    dds <- DESeqDataSetFromMatrix(countData = countdata, colData = colData, design= ~ Genotype.S.R. + Treatmenttime.C.T. + Genotype.S.R.:Treatmenttime.C.T.)
    dds
    dds <- dds[ rowSums(counts(dds)) > 10, ]
    dds
    dds <- DESeq(dds, test = "LRT", reduced = ~ Genotype.S.R. + Treatmenttime.C.T.)
    resultsNames(dds)

The result names I am getting are

> resultsNames(dds)
[1] "Intercept"                                                "Genotype.S.R._Susceptible_vs_Resistant"                  
[3] "Treatmenttime.C.T._ctrl3day_vs_ctrl1day"        "Treatmenttime.C.T._ctrl7day_vs_ctrl1day"       
[5] "Treatmenttime.C.T._trt1day_vs_ctrl1day"         "Treatmenttime.C.T._trt3day_vs_ctrl1day"        
[7] "Treatmenttime.C.T._trt7day_vs_ctrl1day"         "Genotype.S.R.Susceptible.Treatmenttime.C.T.ctrl3day"
[9] "Genotype.S.R.Susceptible.Treatmenttime.C.T.ctrl7day" "Genotype.S.R.Susceptible.Treatmenttime.C.T.trt1day" 
[11] "Genotype.S.R.Susceptible.Treatmenttime.C.T.trt3day"  "Genotype.S.R.Susceptible.Treatmenttime.C.T.trt7day"

I tried to extract "Genotype.S.R.SusceptiblevsResistant" and "Treatmenttime.C.T.trt1dayvsctrl1day" using the following command

> day1 <- results(dds, contrast=list("Genotype.S.R._Susceptible_vs_Resistant", "Treatmenttime.C.T._trt1day_vs_ctrl1day"))

which is giving me a log P values of "Genotype.S.R.SusceptiblevsResistantvsTreatmenttime.C.T.trt1dayvsctrl1day"

But I am trying to get ResistantvsSusceptiblevsTreatmenttime.C.T.trt1dayvs_ctrl1day. For this I tried the relevel option but I am not getting the desired results.

Can anybody please guide me for How to rewrite the code to get my desired output as Genotype.S.R.ResistantvsSusceptiblevsTreatmenttime.C.T.trt1dayvsctrl1day
Deseq2 LRT Factors RNAseq • 1.5k views
ADD COMMENT •
0
Entering edit mode

This is my actual colData

>                       Genotype.S.R. Treatmenttime.C.T.
        87_S_Ctl1_1   Susceptible      ctrl1day
        88_S_Ctl1_2   Susceptible      ctrl1day
        93_S_Trt1_1   Susceptible       trt1day
        94_S_Trt1_2   Susceptible       trt1day
        89_S_Ctl3_1   Susceptible      ctrl3day
        90_S_Ctl3_2   Susceptible      ctrl3day
        95_S_Trt3_1   Susceptible       trt3day
        96_S_Trt3_2   Susceptible       trt3day
        91_S_Ctl7_1   Susceptible      ctrl7day
        92_S_Ctl7_2   Susceptible      ctrl7day
        97_S_Trt7_1   Susceptible       trt7day
        98_S_Trt7_2   Susceptible       trt7day
        46_R_Ctl1_1     Resistant      ctrl1day
        47_R_Ctl1_2     Resistant      ctrl1day
        52_R_Trt1_1     Resistant       trt1day
        53_R_Trt1_2     Resistant       trt1day
        48_R_Ctl3_1     Resistant      ctrl3day
        49_R_Ctl3_2     Resistant      ctrl3day
        54_R_Trt3_1     Resistant       trt3day
        55_R_Trt3_2     Resistant       trt3day
        50_R_Ctl7_1     Resistant      ctrl7day
        51_R_Ctl7_2     Resistant      ctrl7day
        56_R_Trt7_1     Resistant       trt7day
        57_R_Trt7_2     Resistant       trt7day

The colData(dds) after relevel

> dds$Genotype.S.R.<- relevel(dds$Genotype.S.R., "Resistant")
> colData(dds)
DataFrame with 24 rows and 3 columns
           Genotype.S.R. Treatmenttime.C.T.        sizeFactor
               <factor>           <factor>         <numeric>
87_S_Ctl1_1   Susceptible      ctrl1day  1.08864786740592
88_S_Ctl1_2   Susceptible      ctrl1day    1.183334509422
93_S_Trt1_1   Susceptible       trt1day  1.30304310388832
94_S_Trt1_2   Susceptible       trt1day     1.34793137296
89_S_Ctl3_1   Susceptible      ctrl3day  1.55674229437058
...             ...                ...               ...
55_R_Trt3_2     Resistant       trt3day  1.14519994590583
50_R_Ctl7_1     Resistant      ctrl7day  1.22742368966149
51_R_Ctl7_2     Resistant      ctrl7day  1.27590660944462
56_R_Trt7_1     Resistant       trt7day 0.341976417848844
57_R_Trt7_2     Resistant       trt7day 0.270151634860834

I am not sure whether I have executed the relevel option correctly. Kindly help me to figure out this query to get my desired output as Genotype.S.R.ResistantvsSusceptiblevsTreatmenttime.C.T.trt1dayvsctrl1day

ADD REPLY •
1
Entering edit mode

Thanks - I can now see what is happening.

With a variable containing terms 'Susceptible' and 'Resistant', R will automatically set 'Resistant' as the reference level:

test <- c('Susceptible','Susceptible','Resistant','Resistant')
test
[1] "Susceptible" "Susceptible" "Resistant"   "Resistant"  

factor(test)
[1] Susceptible Susceptible Resistant   Resistant  
Levels: Resistant Susceptible

The left-most term ('level') is the reference level, i.e., Resistant

DESeq2 interprets this and then automatically concludes that you want Susceptible vs Resistant, i.e., Resistant as the denominator. Thus, what you actually need to do is:

dds$Genotype.S.R.<- relevel(dds$Genotype.S.R., ref = 'Susceptible')

Can you please try that?

ADD REPLY •
0
Entering edit mode

Thank you for your suggestion. I tried the command you suggested and it worked. I am attaching the code below

dds <- DESeq(dds, test = "LRT", reduced = ~ Genotype.S.R. + Treatmenttime.C.T.)
resultsNames(dds)
dds$Genotype.S.R.<- relevel(dds$Genotype.S.R., ref = 'Susceptible')
dds <- DESeq(dds, test = "LRT", reduced = ~ Genotype.S.R. + Treatmenttime.C.T.)
 resultsNames(dds)
 [1] "Intercept"                                              "Genotype.S.R._Resistant_vs_Susceptible"                
 [3] "Treatmenttime.C.T._ctrl3day_vs_ctrl1day"      "Treatmenttime.C.T._ctrl7day_vs_ctrl1day"     
 [5] "Treatmenttime.C.T._trt1day_vs_ctrl1day"       "Treatmenttime.C.T._trt3day_vs_ctrl1day"      
 [7] "Treatmenttime.C.T._trt7day_vs_ctrl1day"       "Genotype.S.R.Resistant.Treatmenttime.C.T._ctrl3day"
 [9] "Genotype.S.R.Resistant.Treatmenttime.C.T._ctrl7day" "Genotype.S.R.Resistant.Treatmenttime.C.T._trt1day" 
 [11] "Genotype.S.R.Resistant.Treatmenttime.C.T._trt3day"  "Genotype.S.R.Resistant.Treatmenttime.C.T._trt7day" 
  day1 <- results(dds, contrast=list("Genotype.S.R._Resistant_vs_Susceptible", "Treatmenttime.C.T._trt1day_vs_ctrl1day"))
  day1

log2 fold change (MLE): Genotype.S.R._Resistant_vs_Susceptible vs Treatmenttime.C.T._trt1day_vs_ctrl1day 
LRT p-value: '~ Genotype.S.R. + Treatmenttime.C.T. + Genotype.S.R.:Treatmenttime.C.T.' vs '~ Genotype.S.R. + Treatmenttime.C.T.' 
 DataFrame with 46821 rows and 6 columns
                  baseMean     log2FoldChange             lfcSE             stat              pvalue               padj
                 <numeric>          <numeric>         <numeric>        <numeric>           <numeric>          <numeric>
 LOC112783 134.308826562054  0.529460046050067  0.31507263423852 17.4734743533944 0.00368409961553274 0.0303657679001606
LOC112783702  234.852123618286  -0.41513990269723 0.372198899957457 4.77841472846441    0.44351693182079  0.609832013695552

I am getting the log P values for Genotype.S.R.ResistantvsSusceptiblevsTreatmenttime.C.T.trt1dayvsctrl1day

ADD REPLY •
0
Entering edit mode

Great - thanks for providing feedback

ADD REPLY •
0
Entering edit mode
Kevin Blighe ★ 4.0k
@kevin
Last seen 10 minutes ago
The Cave, 181 Longwood Avenue, Boston, …

So, you just want the order of the treatment factor to flip? What is the output of colData(dds)[,'Genotype.S.R.'] before and after you believe you have re-leveled this factor?
ADD COMMENT •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 957 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6