This is a test version of Biostars. For the public version, visit https://www.biostars.org.
multi-factor Deseq LRT test analysis

Hi All, I need help to resolve the design to get differentially expressed gene with deseq2. I have 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 need to compare:

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

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"

My metadata condition

id  Genotype(S/R)   Treatmenttime(C/T)
S_Ctl1_1    Susceptible ctrl1day
S_Ctl1_2    Susceptible ctrl1day
S_Trt1_1    Susceptible trt1day
S_Trt1_2    Susceptible trt1day
S_Ctl3_1    Susceptible ctrl3day
S_Ctl3_2    Susceptible ctrl3day
S_Trt3_1    Susceptible trt3day
S_Trt3_2    Susceptible trt3day
S_Ctl7_1    Susceptible ctrl7day
S_Ctl7_2    Susceptible ctrl7day
S_Trt7_1    Susceptible trt7day
S_Trt7_2    Susceptible trt7day
R_Ctl1_1    Resistant           ctrl1day
R_Ctl1_2    Resistant           ctrl1day
R_Trt1_1    Resistant           trt1day
R_Trt1_2    Resistant           trt1day
R_Ctl3_1    Resistant           ctrl3day
R_Ctl3_2    Resistant           ctrl3day
R_Trt3_1    Resistant           trt3day
R_Trt3_2    Resistant           trt3day
R_Ctl7_1    Resistant           ctrl7day
R_Ctl7_2    Resistant           ctrl7day
R_Trt7_1    Resistant           trt7day
R_Trt7_2    Resistant           trt7day

Can someone please guide me through the code to get my desired output as I am not sure if this is the correct way of doing it.

rna-seq deseq2 lrt design
• 2,622 views • link
updated by tiago211287 tiago211287 • written by aiswaryabioinfo • aiswaryabioinfo

I don't quite understand what exact questions you are trying to answer, but I don't think an LRT is what you want.

• 0 views • link
written by swbarnes2 swbarnes2

1 answer

You can extract the LRT results with:

res <- results(dds)

and use the table or save it as an rds or xlsx file like:

openxlsx::write.xlsx(as.data.frame(res), "LRT.xlsx", row.names = T)
saveRDS(as.data.frame(res), "LRT.rds)

As note: The LRT does not output a fold change value, instead, it shows an adjusted pvalue for the gene that reflects how well it matches your experimental design or not.

One idea to further investigate the LRT results is to produce a heatmap that lets you able to cluster the significant genes and visualize the enriched groups.

Using the ComplexHeatmap package you could (untested code):

library(ComplexHeatmap)
library(dplyr)

res <- as.data.frame(res)
res <- res %>% arrange(padj) %>% filter(padj < 0.1)

mat <- assay(rld)
mat <- mat[rownames(res),]
mat <- scale(t(mat), center = T, scale = T)

set.seed(42)
HM <- Heatmap(t(mat), column_split = 2, row_split = 2)
print(HM)
• 0 views • link
written by tiago211287 tiago211287

Log in to answer this question.

More posts like this