Question

DEseq2 2x2 design using a grouping variable

0

Entering edit mode

weichengz ▴ 10

@weichengz-23557

Last seen 3.5 years ago

Melbourne, Australia

I use DEseq2 to analyse the RNA-seq data out of an animal experiment (two conditions and two gender). The conditions include HS and Control and I am also interested in the condition effect in each in gender (Males or Females), and also sex effect in each condition

As the interaction terms seem complicated, I use the alternative grouping variable ~group to compare specific group. However, there may be something wrong as I got exactly the same number of DEG (P<0.1) when calling for resconmale, ressexcontrol and ressex_HS (Please see scripts below).

Your correction and suggestions are very much appreciated!

Cheers, Wei

#import count file
countall<- read.table(file.choose(),header = F, row.names = 1)

countall <- as.matrix(countall)



#load libraries 
library(DESeq2)


#Assign condition

condition <- factor(c(rep("Control",5),rep("HS",8),rep("Control",7),rep("HS",9)))

sex <- factor(c("female","male","female","male","male",
                "male","male","female","female","female","female","male","female",
                "female","female","female","female","female","male","female",
                "female","female","male","female","male","female","female","male","female"))


#Create a data frame
coldata_all<- data.frame(condition=as.factor(condition),sex=as.factor(sex))


# Construct a DESeqDataSet

dds <- DESeqDataSetFromMatrix(countData=countall,colData=coldata_all, design=~sex+condition)

#combine the factors of interest into a single factor with all combinations of the original factors
dds$group <- factor(paste0(dds$sex, dds$condition))


#change the design to include just this factor, e.g. ~ group
design(dds) <- ~group



#run DEseq
dds <- DESeq(dds)

#Check available comparisons
resultsNames(dds)

#outputs:[1] "Intercept"           "group_femaleHS_vs_femaleControl"    "group_maleControl_vs_femaleControl" "group_maleHS_vs_femaleControl" 

#The effect of condition (HS vs. Control) in female
rescon_female = results(dds,contrast = c("group","femaleHS","femaleControl"))


##The effect of condition (HS vs. Control) in male
rescon_male = results(dds,contrast = c("group","maleHS","maleControl"))


#the sex effect (male vs. female) in Control group

ressex_Control = results(dds,contrast = c("group","maleControl","femaleControl"))


#the sex effect (male vs. female) in HS group
ressex_HS = results(dds,contrast = c("group","maleHS","femaleHS"))

deseq2 ~group • 895 views

ADD COMMENT • link 4.0 years ago weichengz ▴ 10

score 0 · Answer 1 · 2020-06-04

0

Entering edit mode

Kevin Blighe ★ 3.9k

@kevin

Last seen 24 days ago

Republic of Ireland

Can you show some sample output from each results table? - thank you.

ADD COMMENT • link 4.0 years ago Kevin Blighe ★ 3.9k

0

Entering edit mode

Thanks Kevin, outputs for each comparison beblow:

The follow up questions are: The overall comparison code is correct?

The call ressex_Control and ressex_HS seems giving me the same result table. Why?

Thanks, W

> ##The effect of condition (HS vs. Control) in male
> rescon_male = results(dds,contrast = c("group","maleHS","maleControl"))
> rescon_male_sig <- rescon_male[which(rescon_male$padj<0.1),]
> rescon_male_sig
log2 fold change (MLE): group maleHS vs maleControl 
Wald test p-value: group maleHS vs maleControl 
DataFrame with 52 rows and 6 columns
                    baseMean log2FoldChange     lfcSE      stat      pvalue        padj
                   <numeric>      <numeric> <numeric> <numeric>   <numeric>   <numeric>
ENSSSCG00000000084 1987.1604      -0.282329 0.0772782  -3.65341 2.58776e-04 0.090600299
ENSSSCG00000000839  483.5407       0.805278 0.2096528   3.84101 1.22530e-04 0.069388080
ENSSSCG00000000951  562.2993      -1.015217 0.1821303  -5.57412 2.48779e-08 0.000391951
ENSSSCG00000001638  287.4844      -0.309489 0.0787058  -3.93222 8.41645e-05 0.053040467
ENSSSCG00000002517   47.7341      -0.679071 0.1872047  -3.62743 2.86261e-04 0.091450756
...                      ...            ...       ...       ...         ...         ...
ENSSSCG00000038557   93.0551      -9.786734 2.6091998  -3.75086 1.76232e-04 0.077956398
ENSSSCG00000038913  258.0820      -0.356470 0.0943037  -3.78002 1.56815e-04 0.074867179
ENSSSCG00000039546   44.4305       0.883057 0.2350387   3.75707 1.71915e-04 0.077956398
ENSSSCG00000039573  311.4669      -1.412189 0.2641799  -5.34556 9.01392e-08 0.000710072
ENSSSCG00000043115   12.3646      -2.526871 0.6842601  -3.69285 2.21754e-04 0.085514570

#the sex effect (male vs. female) in Control group
> ressex_Control = results(dds,contrast = c("group","maleControl","femaleControl"))
> ressex_Control_sig <- ressex_Control[which(ressex_Control$padj<0.1),]
> ressex_Control_sig
log2 fold change (MLE): group maleControl vs femaleControl 
Wald test p-value: group maleControl vs femaleControl 
DataFrame with 35 rows and 6 columns
                    baseMean log2FoldChange     lfcSE      stat      pvalue        padj
                   <numeric>      <numeric> <numeric> <numeric>   <numeric>   <numeric>
ENSSSCG00000000951  562.2993       1.163727  0.172759   6.73613 1.62655e-11 3.19113e-07
ENSSSCG00000002329   75.4451       0.743164  0.177161   4.19484 2.73061e-05 2.55104e-02
ENSSSCG00000006415   24.6487       1.495953  0.398680   3.75226 1.75245e-04 9.93738e-02
ENSSSCG00000008164 1395.8176       0.480559  0.117042   4.10585 4.02826e-05 3.29294e-02
ENSSSCG00000009298  370.9642       0.511900  0.136530   3.74937 1.77281e-04 9.93738e-02
...                      ...            ...       ...       ...         ...         ...
ENSSSCG00000038128 103.49179       0.786437  0.194402   4.04543 5.22274e-05 0.040985982
ENSSSCG00000039573 311.46687       1.352061  0.250226   5.40336 6.54022e-08 0.000320782
ENSSSCG00000040260  31.98008       1.222878  0.313508   3.90063 9.59424e-05 0.068415461
ENSSSCG00000044223   1.45022       4.721301  1.248731   3.78088 1.56275e-04 0.095407727
ENSSSCG00000048182  20.63859       2.380074  0.610843   3.89638 9.76417e-05 0.068415461

#the sex effect (male vs. female) in HS group
> ressex_HS = results(dds,contrast = c("group","maleHS","femaleHS"))
> ressex_HS_sig <- ressex_HS[which(ressex_HS$padj<0.1),]
> ressex_Control_sig
log2 fold change (MLE): group maleControl vs femaleControl 
Wald test p-value: group maleControl vs femaleControl 
DataFrame with 35 rows and 6 columns
                    baseMean log2FoldChange     lfcSE      stat      pvalue        padj
                   <numeric>      <numeric> <numeric> <numeric>   <numeric>   <numeric>
ENSSSCG00000000951  562.2993       1.163727  0.172759   6.73613 1.62655e-11 3.19113e-07
ENSSSCG00000002329   75.4451       0.743164  0.177161   4.19484 2.73061e-05 2.55104e-02
ENSSSCG00000006415   24.6487       1.495953  0.398680   3.75226 1.75245e-04 9.93738e-02
ENSSSCG00000008164 1395.8176       0.480559  0.117042   4.10585 4.02826e-05 3.29294e-02
ENSSSCG00000009298  370.9642       0.511900  0.136530   3.74937 1.77281e-04 9.93738e-02
...                      ...            ...       ...       ...         ...         ...
ENSSSCG00000038128 103.49179       0.786437  0.194402   4.04543 5.22274e-05 0.040985982
ENSSSCG00000039573 311.46687       1.352061  0.250226   5.40336 6.54022e-08 0.000320782
ENSSSCG00000040260  31.98008       1.222878  0.313508   3.90063 9.59424e-05 0.068415461
ENSSSCG00000044223   1.45022       4.721301  1.248731   3.78088 1.56275e-04 0.095407727
ENSSSCG00000048182  20.63859       2.380074  0.610843   3.89638 9.76417e-05 0.068415461

ADD REPLY • link 4.0 years ago weichengz ▴ 10

0

Entering edit mode

Just note that, in the bottom table, you are actually showing ressex_Control_sig results, not ressex_HS_sig

ADD REPLY • link 4.0 years ago Kevin Blighe ★ 3.9k

0

Entering edit mode

Thanks for the correction. BTW, are the codes above correct for each specific comparison?

ADD REPLY • link 4.0 years ago weichengz ▴ 10

0

Entering edit mode

Yes, they seem fine; although, you are not doing any log-fold change shrinkage step. You can check the DESeq2 vignette for details on that: http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

ADD REPLY • link 4.0 years ago Kevin Blighe ★ 3.9k