Question

RNA-Seq analysis edgeR

0

Entering edit mode

eb0906 • 0

@eb0906-7337

Last seen 9.7 years ago

United States

We have a time course infection experiment with 12 samples (2 technical replicates and 2 biological replicates per sample). The samples were sequenced on a HiSeq over 2 lanes, 1 x 50 bp reads.

24hr_infection1_lane3
24hr_infection1_ lane4
24hr_infection2_lane3
24hr_infection2_lane4
24hr_control1_lane3
24hr_control1_lane4
24hr_control2_lane3
24hr_control2_lane4
48hr_infection1_lane3
48hr_infection1_lane4
48hr_infection2_lane3
48hr_infection2_lane4

I aligned the samples using Tophat and then converted the .bam files to .sam files and created a counts file using HTSeq for each sample. My questions now are:

Should the technical replicates be collapsed (add gene counts for technical replicates)?
How should the design matrix be setup in edgeR?

edgeR rna-seq replicates • 2.5k views

ADD COMMENT • link updated 10.0 years ago by James W. MacDonald 67k • written 10.0 years ago by eb0906 • 0

score 2 · Accepted Answer · 2015-02-06

2

Entering edit mode

James W. MacDonald 67k

@james-w-macdonald-5106

Last seen 1 day ago

United States

Note that Tophat is intended for finding splice junctions, usually with the goal that you would then infer transcript abundances and presumably then infer differential transcript splicing. This is different from what one normally does with edgeR, which is to infer differential gene expression (which by default ignores differences in the form of the transcripts).

The former is really hard to do, and requires both high read depth and lots of biological replication. The latter is much less difficult, as it is a simpler problem. My point being that Tophat spends a lot of time figuring out things that you are planning to ignore, so you are likely better off using something like bowtie2 or subread or bwa, etc to do the aligning. subread in particular is very fast for this sort of thing, and there is a Bioconductor Rsubread package that makes it simple to use.

Anyway, the answers to your question are:

1.) Yes

2.) It depends on what you are trying to detect. I don't see a 48 hour control, so presumably you are planning to compare both time points back to the 24 hour control. In which case you would want to do something like

design <- model.matrix(~factor(rep(c("inf_24","control","inf_48"), each = 2)))

which parameterizes the second and third coefficient as 24h infection - control, and 48h infection - control, respectively.

ADD COMMENT • link 10.0 years ago James W. MacDonald 67k

0

Entering edit mode

Thanks very much for your input.

Based on your response, we are thinking that we should redo the alignment using bowtie2 and use the .sam output for creating a counts file using HTSeq for each sample and then collapse technical replicates. Yes, you are correct in assuming that we will compare both time points to the 24 hr control since we do not have an additional control at 48 hrs.

To be clear, if one is only interested in differential gene expression between treatment groups, then the Tuxedo suite tools is not a good choice? Ideally, one would use a different aligner (BWA, Bowtie, STAR, etc) and then use either edgeR or DESeq2?

ADD REPLY • link 10.0 years ago eb0906 • 0

0

Entering edit mode

I am saying that Tophat is probably not the tool you want for the task at hand. It uses bowtie to do the alignment, and then goes back and tries to find junction splicing reads. If you don't care about the junction splicing reads (which are in a different file than your aligned reads), then you are doing an extra step, which takes that much more time, and then not using the data that you took the extra time to generate. I don't think that the accepted-hits.bam file contains the junction-spanning reads, so you may be losing those.

Using bowtie2 for the aligner is helpful because it is a gapped aligner, so you shouldn't lose reads that span exon-exon junctions, which you might lose with bowtie1. IIRC, bwa and star are gapped aligners as well. Subread just soft-clips the junction spanning reads, but that is a distinction without a difference. You will still get the same number of reads per gene (because subread says 'the read aligns to the end of this exon', and the gapped aligners say 'the read aligns to the ends of these two exons', but both result in that read being counted as overlapping a gene).

So to answer your question, yes, ideally you would just use one of the gapped aligners that you mention and then use either edgeR or DESeq2 to test for differential expression.

ADD REPLY • link 10.0 years ago James W. MacDonald 67k

0

Entering edit mode

And I forgot to mention that you can just use bowtie2 (or pretty much any aligner) to align the technical replicates together, so you don't have to collapse after the alignment. You just do something like

bowtie2-align -U techrep1.fastq,techrep2.fastq | samtools view -bS - > output.bam

or if you have paired end

bowtie2-align -1 techrep1-1.fastq,techrep2-1.fastq -2 techrep1-2.fastq,techrep2-2.fastq | samtools view -bS - > output.bam