Hi,
I am having some problems with my DEseq2 experimental design (RNA-seq raw read counts). Basically, I have two treatment (sf, ls). Within each treatment, there are six cows with 3 time points(14,7 and 28). Time point 14 is the control. I want to find DE genes, 28 vs 14, 7 vs 14.
I first tried a design: ~ treat + treat:cowf + treat:day
Here is the script to find DE genes of 28 vs 14 in ls treatment:
I am not sure if the above design is correct. Since the two treatment are independent, i used data from each treatment separately for DE analysis, which is easier as I can remove the treatment factor.
Below is for ls treatment with design ~cow + day
However, results from the two methods above produced different results. Could anybody explain me which one I should use? Thanks
I am having some problems with my DEseq2 experimental design (RNA-seq raw read counts). Basically, I have two treatment (sf, ls). Within each treatment, there are six cows with 3 time points(14,7 and 28). Time point 14 is the control. I want to find DE genes, 28 vs 14, 7 vs 14.
I first tried a design: ~ treat + treat:cowf + treat:day
Code:
cow treat day cowf 14LS2 2 ls 14 1 14LS63 63 ls 14 2 14LS66 66 ls 14 3 14LS67 67 ls 14 4 14LS73 73 ls 14 5 14LS74 74 ls 14 6 28LS2 2 ls 28 1 28LS63 63 ls 28 2 28LS66 66 ls 28 3 28LS67 67 ls 28 4 28LS73 73 ls 28 5 28LS74 74 ls 28 6 7LS2 2 ls 7 1 7LS63 63 ls 7 2 7LS66 66 ls 7 3 7LS67 67 ls 7 4 7LS73 73 ls 7 5 7LS74 74 ls 7 6 14SF5355 5355 sf 14 1 14SF61 61 sf 14 2 14SF62 62 sf 14 3 14SF70 70 sf 14 4 14SF71 71 sf 14 5 14SF72 72 sf 14 6 28SF5355 5355 sf 28 1 28SF61 61 sf 28 2 28SF62 62 sf 28 3 28SF70 70 sf 28 4 28SF71 71 sf 28 5 28SF72 72 sf 28 6 7SF5355 5355 sf 7 1 7SF61 61 sf 7 2 7SF62 62 sf 7 3 7SF70 70 sf 7 4 7SF71 71 sf 7 5 7SF72 72 sf 7 6
Here is the script to find DE genes of 28 vs 14 in ls treatment:
Code:
countdata <- read.csv("readcounts_matrix.csv", sep=",",stringsAsFactors=FALSE, row.names=1) coldata <- read.csv("coldata.csv", sep=",",stringsAsFactors=TRUE, row.names=1) coldata$cow <- factor(coldata$cow) coldata$day <- factor(coldata$day) coldata$cowf <- factor(rep(c(1,2,3,4,5,6),times=6)) coldata$day <-relevel(coldata$day, "14") dds <- DESeqDataSetFromMatrix(countData = countdata, colData = coldata, design = ~ treat + treat:cowf + treat:day) dds <- estimateSizeFactors(dds) nc <- counts(dds, normalized=TRUE) filter <- rowSums(nc >= 10) >= 12 dds <- dds[filter,] dds <- DESeq(dds,modelMatrixType="standard") ls_d28vsd14 <- results(dds, name = "treatls.day28")
Below is for ls treatment with design ~cow + day
Code:
cow treat day cowf 14LS2 2 ls 14 1 14LS63 63 ls 14 2 14LS66 66 ls 14 3 14LS67 67 ls 14 4 14LS73 73 ls 14 5 14LS74 74 ls 14 6 28LS2 2 ls 28 1 28LS63 63 ls 28 2 28LS66 66 ls 28 3 28LS67 67 ls 28 4 28LS73 73 ls 28 5 28LS74 74 ls 28 6 7LS2 2 ls 7 1 7LS63 63 ls 7 2 7LS66 66 ls 7 3 7LS67 67 ls 7 4 7LS73 73 ls 7 5 7LS74 74 ls 7 6
Code:
ls_countdata = countdata[,1:18] ls_coldata = coldata[1:18,] ls_dds <- DESeqDataSetFromMatrix(countData = ls_countdata, colData = ls_coldata, design = ~ cow + day) ls_dds <- estimateSizeFactors(ls_dds) nc <- counts(ls_dds, normalized=TRUE) filter <- rowSums(nc >= 10) >= 6 ls_dds <- ls_dds[filter,] ls_dds$day <- relevel(ls_dds$day, "14") ls_dds <- DESeq(ls_dds) ls_res_28vs14 <- results(ls_dds,contrast=c("day", "28","14"))
Comment