RNA-seq time course experiment - LRT design hep
1
0
Entering edit mode
@maithe-barros-13840
Last seen 5 months ago

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!

deseq2 lrt timecourse • 461 views
ADD COMMENT
0
Entering edit mode

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?

ADD REPLY
1
Entering edit mode

I moved this from an Answer to a Comment.

Yes that’s what I meant by combining variables and using contrast.

ADD REPLY
0
Entering edit mode

Awesome! Thank you ever so much for your quick reply and for helping me out, Michael! 

ADD REPLY
1
Entering edit mode
@mikelove
Last seen 4 hours ago
United States

It sounds like you could accomplish your desired comparisons by combing time and condition into a single factor (see vignette section on interactions). So then you’d use ~horse + group, and make pairwise comparisons with results() and contrast.

ADD COMMENT

Login before adding your answer.

Traffic: 493 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6