I have RNASeq data that I am analysing using DESeq2. I have 2 genotypes (WT and mutant). The RNASeq was run in 3 separate batches (Batch1, 2 and 3; all batches contained some WT and mutant samples). I ran DESeq using design=~genotype + batch, re-leveled to genotype=WT. I then tried including surrogate variables based on the batch using SVA as follows below. I get a different number of DE genes in each analysis.
When I run a PCA and sample distances, the samples cluster more by batch than by genotype (particularly the 3rd batch where I know the quality of the RNA was poor).
I'm not sure which results I should be actually using. It seems that using SVA is a better way to exclude technical variation. Am I correct? When applying SVA, should the batch be included in the original design, or just the genotype?
I appreciate any help! Thanks in advance!
dds <- estimateSizeFactors(dds)
dat <- counts(dds, normalized=TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx,]
mod <- model.matrix(~ genotype, colData(dds))
mod0 <- model.matrix(~ 1, colData(dds))
svseq <- svaseq(dat, mod, mod0, n.sv=2)
# output:
# Number of significant surrogate variables is: 2
ddssva <- dds
ddssva$SV1 <- svseq$sv[,1]
ddssva$SV2 <- svseq$sv[,2]
design(ddssva) <- ~ SV1 + SV2 + genotype
ddssva <- DESeq(ddssva)
ressva <- results(ddssva,contrast=c("genotype","HET","WT"))
•
written
by
Chris Gene
2 answers
Update: After looking at your colData.
Create matrix for svaseq like this:
mod <- model.matrix(~ dex, colData(dds))
mod0 <- model.matrix(~ 1, colData(dds))
Then do unsupervised svaseq
svseq <- svaseq(dat, mod, mod0, n.sv=2)
Now use this as input to DESeq design, along with genotype
dds$SV1 <- svseq$sv[,1]
dds$SV2 <- svseq$sv[,2]
design(dds) <- ~ genotype + SV1 + SV2
I wouldn't expect as much difference as you told, if you implement both methods properly..
Estimating batches by svaseq followed by providing covariates in the DESeq design is expected to perform either equally good or better than adding batch as a factor in glm.. So design matrix for deseq could be something like:
model.matrix(as.formula(paste0("~ genotype",svseq))
where svseq is the output from svseq. Then you can continue with DESeq ..
Different number of significant genes is expected since the covariate estimation in different in svseq vs blocking design (multifactor glm). I wouldn't expect too much difference though..
Log in to answer this question.
-
How to correct for known and surrogate variables in DESeq2?
written by Crystal571 •
-
SVA correction reduces ability to detect genes as differential expressed?
written by MB •
-
SVA and ruvseq error: system is computationally singular
written by Jasmin •
-
Partial NA coefficients warning message using svaseq and deseq
written by Jasmin •
-
Batch correction: Before and After Visualisation
written by prabin.dm
-
Design in DESeq: Can you combine explicit and implicit batch effect correction in (SVA+Metadata Var…
written by Kristin Muench
-
Can't make sva comparison despite very distinct batch effects
written by jjrin •
-
SVA + DESEQ2 : Looking for batcheffect in RNA-SEQ
written by ZheFrench
-
why my resultsNames(deseq) is not the same as model.matrix ?
written by unique379
-
Sva batch correction for rna seq
written by iside •
Thank you for your reply.
As I understand, you are suggesting I use
~ genotype,svseqinstead of genotype, colData in the model matrix as I did above?The different number is something like 40 to 400. The numbers aren't huge to start, but it still is a 10 fold difference. Is that unreasonable?
Many thanks for your help
best wishes
Chris
Hi Vivek,
Sorry to bother again, but when I try to do the above it won't let me add svseq to the design.
error:
My colData contains the information on sample, genotype and batch:
Yes, since I assumed a very simple colData design for demo. I updated my answer now. But you don't need to put "sample" as a column in there.. You wouldn't need even the "Run" column if you are estimating batches in unsupervised way.
Thank you Vivek. Very helpful.