Dealing With Read Counts Under Pe And Se Scenarios

Question

[Arun]

Hi,

I am unsure how to deal with this case to go about analysing RNA-seq data.

Suppose that you have a control and treatment setup with 4 biological replicates each. However, two in control and two in treatment were pooled together and paired-end sequenced and the other 4 were single-end sequenced. That is:

Condition    Sample_No    Sequence_Type
Control        1,2          paired-end
Control        3,4          single-end
Treatment      1,2          paired-end
Treatment      3,4          single-end

Now, suppose I'd want to perform a differential gene expression analysis, how would you take care of the difference in the reads (due to PE and SE) within conditions?

i) You map the reads as such - PE as PE and SE as SE libraries. You count the total number of reads that fall under each sample. You then normalize using edgeR's TMM method for difference in the library size and then perform the differential expression analysis.
ii) You map the reads as such - PE as PE and SE as SE libraries. You count the total number of reads first in pair (meaningful for PE samples) that fall under each sample. You then normalize using edgeR's TMM method for difference in the library size and then perform the differential expression analysis.
iii) You discard the second pair altogether and treat them as two conditions with 8 SE libraries.

Understanding the inherent mess-up in the experimental setup and the possibility of bias and accepting them, which one would you prefer and why?

I personally prefer (i) because, I think PE mappings are much more efficient than SE mappings. So, I don't want to lose this info (rules out iii). And this allows us for a better estimate of gene expression in two of the four replicates under each condition. Then by counting the total number of reads, and normalising for the difference in the library sizes, we attempt to compensate for the difference. The SE library would have the same pattern of expression genome-wide within each condition and therefore we just normalise for library size (just a scaling), which I guess makes sense. And I don't think there is a difference between (i) and (ii) except maybe (i) has more power during statistical testing due to more confidence, although not quite sure if its desirable. What do you guys think?

Answer 1

[Istvan Albert]

One possible problem could arise from reads mapping to different locations in the two modes - that in turn could lead to different coverages simply due to the PE and SE modes. I don't know for sure if this is a real problem but it could be worth keeping in mind. Perhaps you could do the entire analysis separately for PE and SE modes and compare/consolidate the results only.

[Arun]

Istvan, thank you. Yes, I thought of the issue you mention as well. However, I couldn't think of getting around it. Since I have more than 1 replicate in both cases, your solution to consolidate results by comparing PE/SE with PE/SE makes total sense! I'll try them out. Thanks a lot! I'll wait a while to accept the answer.

[Arun]

Istvan, just a follow up question. The sequenced data was from Illumina GAII sequencer and had lot of samples multiplexed. So, the throughput wasn't quite high. So, I pooled the PE samples together and SE samples together within conditions for mapping (using tophat) and after mapping demultiplexed them. For SE reads, I provide the mapping junctions file (junctions.bed) obtained from PE run as input to single end files. So, I run the tophat pipeline 2 times and provide a junctions file as input to the SE case. I'd suppose the mapping is better than not providing it in the SE case. I hope this should compensate for the coverage due to SE and PE differences to some extent, the concern you raised?

[Istvan Albert]

Sounds good. In general I found that begin defensive helps avoid many/most of the problems.