Hello,
I am doing RNA-seq analysis of DE genes from mare's endometrial biopsies. We collected our samples at an abattoir due to restrictions with samples collected from live horses.
Thus, I have three time points: alive_0h (when the mare was just slaughtered, tissue representing the uterus of a live mare still). Then we took biopsies from the uterus and used them in a tissue culture. After culturing the explants we had 2 time points (24h and 48h) for both control and LPS challenge.
I am having a hard time trying to figure out the design I should use since I have this samples at 0h. We have 5 biological replicates (horses).
My LRT design is the following one:
ddsLRT_exp3 <- DESeqDataSet(se_exp3, design = ~ horse + treatment)
ddsLRT_exp3 <- estimateSizeFactors(ddsLRT_exp3)
ddsLRT_exp3 <- estimateDispersions(ddsLRT_exp3)
ddsLRT_exp3 <- nbinomLRT(ddsLRT_exp3, reduced = ~ horse)
resLRT_exp3 <- results(ddsLRT_exp3)
My colData is:
samples,treatment,horse
sample_55.bam,alive_0h,7B
sample_56.bam,control_24h,7B
sample_57.bam,LPS_24h,7B
sample_58.bam,control_48h,7B
sample_59.bam,LPS_48h,7B
sample_60.bam,alive_0h,13D
sample_61.bam,control_24h,13D
sample_62.bam,LPS_24h,13D
sample_63.bam,control_48h,13D
sample_64.bam,LPS_48h,13D
sample_65.bam,alive_0h,14B
sample_66.bam,control_24h,14B
sample_67.bam,LPS_24h,14B
sample_68.bam,control_48h,14B
sample_69.bam,LPS_48h,14B
sample_70.bam,alive_0h,15A
sample_71.bam,control_24h,15A
sample_72.bam,LPS_24h,15A
sample_73.bam,control_48h,15A
sample_74.bam,LPS_48h,15A
sample_75.bam,alive_0h,15C
sample_76.bam,control_24h,15C
sample_77.bam,LPS_24h,15C
sample_78.bam,control_48h,15C
sample_79.bam,LPS_48h,15C
Our aims are to compare:
1) Control vs LPS at 24h and 48h
2) LPS 24h vs 48h
3) Alive 0h vs Control 24h vs Control 48h
Many thanks indeed!
Thank you so much for the quick reply! I did think about it and I already ran something like that:
bam_exp3 <-c("sample_55.bam", "sample_56.bam", "sample_57.bam","sample_58.bam", "sample_59.bam", "sample_60.bam", "sample_61.bam", "sample_62.bam", "sample_63.bam", "sample_64.bam", "sample_65.bam", "sample_66.bam", "sample_67.bam", "sample_68.bam", "sample_69.bam", "sample_70.bam", "sample_71.bam", "sample_72.bam", "sample_73.bam", "sample_74.bam", "sample_75.bam", "sample_76.bam", "sample_77.bam", "sample_78.bam", "sample_79.bam")
library(Rsamtools)
bamfiles <- BamFileList(bam_exp3)
library(GenomicFeatures)
txdb_exp3 <- makeTxDbFromGFF("Equus_caballus.EquCab2.89.gtf", format="gtf")ebg_exp3 <- exonsBy(txdb_exp3, by="gene")
se_exp3 <- summarizeOverlaps(features=ebg_exp3, reads=bamfiles,
mode="Union",
singleEnd=FALSE,
ignore.strand=TRUE,
fragments=TRUE )
dim(se_exp3)
assayNames(se_exp3)
colSums(assay(se_exp3))
rowRanges(se_exp3)
sampleTableExp3 <- read.csv("exp3_metadata.csv")
colData(se_exp3) <- DataFrame(sampleTableExp3)
library(DESeq2)
dds_exp3 <- DESeqDataSet(se_exp3, design = ~ horse + treatment)
dds_exp3 <- DESeq(dds_exp3)
resultsNames(dds_exp3)
plotDispEsts(dds_exp3, main="Estimation of dispersion", sub="For DESeq2 object")
res_exp3_control24h_LPS24h <- results(dds_exp3, contrast=c("treatment", "control_24h", "LPS_24h"))
res_exp3_LPS24h_LPS48h <- results(dds_exp3, contrast=c("treatment", "LPS_24h", "LPS_48h"))
res_exp3_control48h_LPS48h <- results(dds_exp3, contrast=c("treatment", "control_48h", "LPS_48h"))
Would you give me your opinion? Does that look correct to you?
I moved this from an Answer to a Comment.
Yes that’s what I meant by combining variables and using contrast.
Awesome! Thank you ever so much for your quick reply and for helping me out, Michael!