The support.bioconductor.org editor has been updated to markdown! Please see more info at: Tutorial: Updated Support Site Editor

Question: DESeq2 with interaction using SVA
3
gravatar for kevin.kingsland
12 months ago by
kevin.kingsland30 wrote:

Hello all,

I tried searching for this specific problem but could not find the answer.  I will happily accept a link to a previous question if such exists.

To keep this simple, I can go into more details in response, I am using DESeq2 to analyze the effect of heat stress on expression in two closely-related fish species.  I want to determine if the species vary in their response to the heat stress. So the model I used in DESeq2 was:

dds=DESeqDataSetFromMatrix(countData=countmatrix, colData=coldata, design= ~species+treatment+species:treatment)
dds$test <-factor(paste0(dds$species,dds$treatment))
design(dds) <- ~test
dds <-DESeq(dds)

I also want to eliminate any potential hidden batch effects using SVA.  So I followed the code in the DESeq2 walkthrough.

I created a new DESeqDataSet and included the two surrogate variables in the model and the grouping variable from the original DESeqDataSet:

ddssva <- dds
ddssva$SV1 <- svseq$sv[,1]
ddssva$SV2 <- svseq$sv[,2]
design(ddssva) <- ~ SV1 + SV2 + test
ddssva <- DESeq(ddssva)

My question is, is this the correct model to use to test if the species vary in their response to the heat stress, considering the effects of the surrogate variables?


I'll include the relevant code below. Thanks much, and forgive mistakes in etiquette as this is my first question post.

Sincerely,

Kevin

countmatrix<-read.csv(file="C:\\Kevin\\UniqueCounts.csv" ,header=TRUE, row.names=1)
as.matrix(countmatrix[,-1])
colnames(countmatrix) <- NULL
coldata = data.frame(row.names = c('pallidC-1', 'pallidC-2', 'pallidC-3', 'pallidC-4', 'pallidC-5', 'shovelC-1', 'shovelC-2',  'shovelC-3','pallidH-1','pallidH-2','pallidH-3','pallidH-4','pallidH-5','shovelH-1','shovelH-2','shovelH-3'),species=c("pallid","pallid","pallid","pallid","pallid","shovelnose","shovelnose","shovelnose","pallid","pallid","pallid","pallid","pallid","shovelnose","shovelnose","shovelnose"),treatment=rep(c("control","heat"),each=8))
dds=DESeqDataSetFromMatrix(countData=countmatrix, colData=coldata, design= ~species+treatment+species:treatment)
dds$test <-factor(paste0(dds$species,dds$treatment))
design(dds) <- ~test
dds <-DESeq(dds)

dat <- counts(dds, normalized=TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx,]
mod <- model.matrix(~ test, colData(dds))
mod0 <- model.matrix(~ 1, colData(dds))
svseq <- svaseq(dat, mod, mod0, n.sv=2)
    Number of significant surrogate variables is:  2
    Iteration (out of 5 ):1  2  3  4  5  
svseq$sv
par(mfrow=c(2,1),mar=c(3,5,3,1))
stripchart(svseq$sv[,1] ~ dds$species,vertical=TRUE,main="SV1")
abline(h=0)
stripchart(svseq$sv[,2] ~ dds$species,vertical=TRUE,main="SV2")
abline(h=0)
ddssva <- dds
ddssva$SV1 <- svseq$sv[,1]
ddssva$SV2 <- svseq$sv[,2]
design(ddssva) <- ~ SV1 + SV2 + test
ddssva <- DESeq(ddssva)

Note: the "UniqueCounts.csv" file contains 49,126 rows of data, one for each gene/transcript in my assembled transcriptome, and 16 columns of data corresponding to the 16 RNA libraries used in the study.

ADD COMMENTlink modified 12 months ago • written 12 months ago by kevin.kingsland30
Answer: DESeq2 with interaction using SVA
1
gravatar for Michael Love
12 months ago by
Michael Love22k
United States
Michael Love22k wrote:

Yes that’s the correct code. 

ADD COMMENTlink written 12 months ago by Michael Love22k

Thank you Michael.  I'll be sure to submit any additional questions that may pop up.

ADD REPLYlink written 12 months ago by kevin.kingsland30
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 319 users visited in the last hour