Hi,
I'm trying to identify DEG between different conditions at two time points.
Here is my design :
group,ind,condition
4h,4h_cond1_1,cond1_4h
4h,4h_cond1_2,cond1_4h
4h,4h_cond1_3,cond1_4h
4h,4h_cond1_4,cond1_4h
4h,4h_cond2_1,cond2_4h
4h,4h_cond2_2,cond2_4h
4h,4h_cond2_3,cond2_4h
4h,4h_cond2_4,cond2_4h
8h,8h_cond1_1,cond1_8h
8h,8h_cond1_2,cond1_8h
8h,8h_cond1_3,cond1_8h
8h,8h_cond1_4,cond1_8h
8h,8h_cond2_1,cond2_8h
8h,8h_cond2_2,cond2_8h
8h,8h_cond2_3,cond2_8h
8h,8h_cond2_4,cond2_8h
I tried to use the "Group-specific condition effects, individuals nested within groups" vignette from the documentation to include two batch effects (group and ind).
Here is the code I used
> raw_counts <- read.csv("counts_raw.csv",sep=",",header=T)
> file_des <- read.csv("design_batch.csv", sep=",", header=T)
> cts <- as.matrix(raw_counts)
> coldata <- DataFrame(grp=factor(factor(file_des$group), ind=factor(file_des$ind), cnd=factor(file_des$condition))
> dds <- DESeqDataSetFromMatrix(countData = cts,colData = coldata,design = ~ cnd)
> coldata$ind.n <- factor(rep(rep(1:4,each=1),4))
> colData(dds) = coldata
> model.matrix(~ grp + grp:ind.n + grp:cnd, coldata)
----------------
-----PCA-----
----------------
> vsd <- vst(dds)
> assay(vsd) <- limma::removeBatchEffect(assay(vsd), vsd$grp, vsd$ind.n)
> plotPCA(vsd, intgroup = c("cnd"))
----------------
-----DEG----
----------------
> dds <- DESeq(dds)
> results(dds, name=c("cnd_cond1_4h_vs_cond1_8h"))
sessionInfo( )
R version 4.3.1 (2023-06-16)
other attached packages:
[1] DESeq2_1.40.2 SummarizedExperiment_1.30.2
[3] Biobase_2.60.0 MatrixGenerics_1.12.2
[5] matrixStats_1.0.0 GenomicRanges_1.52.0
[7] GenomeInfoDb_1.36.1 IRanges_2.34.1
[9] S4Vectors_0.38.1 BiocGenerics_0.46.0
[11] limma_3.58.1
When I remove the two batch effects for the PCA, I can see that there is a difference with and without the batch effect, but in the differential analysis there is no difference in the results.
Is my approach correct? Or is there another way to include two batch effects in the analysis?
Thanks for your help.
Rania