Why of number of F3 and F5 reads differ a lot of my solid data

[chenyao]

Dear all,

we recently performed the paire-end solid RNA-seq data on several samples. However, after mapping by Tophat, I found the number of two reads differ a lot. As you can see, the flagstat of one sample:

67834275 in total
0 QC failure
0 duplicates
67834275 mapped (100.00%)
67834275 paired in sequencing
17769061 read1
50065214 read2
16068482 properly paired (23.69%)
17273682 with itself and mate mapped
50560593 singletons (74.54%)
-----------------
I don't know why, If the two kinds of read are not in similar number, why using pair-end data. How can I deal with so many singletons?

Need help!

[sdvie]

Did you check the number of reads in your original 2 reads files (F3.csfasta and F5.csfasta) as well? Do you observe the same difference there?

[chenyao]

Did you check the number of reads in your original 2 reads files (F3.csfasta and F5.csfasta) as well? Do you observe the same difference there?

good idea! Any command can check the number of reads in this color-space data?
Another question, Should the two reads of pair end data have similar numbers?

[sdvie]

You could try a simple unix line count command:

try:

wc -l reads_F3.csfasta

The associated reads_F5.csfasta file should have the same number of reads.

Another thing you could try, if you are not sure whether your experiment was done correctly is submit your reads to a quality control software like FastQC

Keep in mind that in the csfasta, each read has 2 lines. So to know how many reads are really in your sample, you would have to divide the number by 2.

Good luck!

[chenyao]

You could try a simple unix line count command:

try:

wc -l reads_F3.csfasta

The associated reads_F5.csfasta file should have the same number of reads.

Another thing you could try, if you are not sure whether your experiment was done correctly is submit your reads to a quality control software like FastQC

Keep in mind that in the csfasta, each read has 2 lines. So to know how many reads are really in your sample, you would have to divide the number by 2.

Good luck!

Thank you. Can "wc" work? I have tried wc -l on "bam" file, the number of reads is wrong. Alos FastQC can not work on Color-space data, I have tried.

[DZhang]

wc only works on text files.

[DZhang]

Chenyao,

As to your original question, you need to discuss with your sequencing provider why the read2 reads have so many ones are _NOT_ mapped in proper pairs. I am not sure what it could be but very interested in learning that.

[kopi-o]

You have to be a bit careful with interpreting TopHat output. TopHat will flag reads that are too far away from their mate, according to its criteria, as not properly paired. However, the default criteria used by TopHat are rather strict. If you would use BWA with default parameters, I am betting that you would get a higher number of properly paired reads. The parameters for this can be adjusted in TopHat, though.

The thing that jumps out at me is

17769061 read1
50065214 read2

suggesting that read2 mapped much better than read1. I wonder what that could be about (usually F5, which I assume is read2, maps worse)

[DZhang]

It could be that the read2 reads have wide-spread miscallings and/or Ns so they have a much worse mapping rate.

[chenyao]

Chenyao,

As to your original question, you need to discuss with your sequencing provider why the read2 reads have so many ones are _NOT_ mapped in proper pairs. I am not sure what it could be but very interested in learning that.

That really bothers me. Here is their explanation:

There are on average ~35M F3 reads per sample that, individually, map uniquely to the genome, and approximately 8M F5 reads that, individually, map uniquely to the genome. Please also note that for the
"Uniquely-Mapping Pairs", the number is smaller than the number of "F3 reads
uniquely mapped", but larger than the number of "F5 reads uniquely mapped". It is smaller than the former mostly because many of the F5 reads don't map uniquely and cannot be rescued from a set of possible locations based on proximity to the F3 location (or they don't map to the mouse genome at all). It is larger than the latter because many F5 reads that map to multiple sites in the genome are rescued because one of those locations is near to the location of the F3 read on its bead.
---------------------
I don't get it. I suspect because F3 is 50bp and F5 is 25bp.

[chenyao]

You have to be a bit careful with interpreting TopHat output. TopHat will flag reads that are too far away from their mate, according to its criteria, as not properly paired. However, the default criteria used by TopHat are rather strict. If you would use BWA with default parameters, I am betting that you would get a higher number of properly paired reads. The parameters for this can be adjusted in TopHat, though.

The thing that jumps out at me is

17769061 read1
50065214 read2

suggesting that read2 mapped much better than read1. I wonder what that could be about (usually F5, which I assume is read2, maps worse)

Yes, I think the first read I used in tophat is F3, and this result is opposite with the short report from the sequencing provider. Strangely, I got similar results (F5 mapped better than F3) for 5 differnt samples.

[chenyao]

Strangely, I check the original "csfasta" file by "wc -l". The F3 and F5 reads have equal numbers. So I am more confused.

[pmiguel]

There is some history to this, but the sad truth is that at a least as of v4, SOLiD F5 reads are of much lower quality than F3 reads. It is not surprising that F5 reads map much more poorly than F3.

We immediately purchased the new paired end chemistry as soon as ABI released it for the SOLiD (with the release of v3 chemistry?) Unlike Illumina PE chemistry, no new template strand is created. Ligases are more or less 3'->5'/5'->3' blind. So conceptually this had always seemed like a strong point for the SOLiD. However, as with many things, the implementation proved more difficult than the concept would have suggested.

Nevertheless our first experiences with the new chemistry were very positive. Cost for the new PE kits was very reasonable. Mapping percentages were about what we thought we would expect for the shorter reads.

The only problem was that orders of the PE chemistry took a very long time to be filled. After many bitter complaints, I was told that an outside provider of a reagent critical for the the chemistry had failed to keep up with ABI's demand for this reagent. ABI had finally in-housed production of this reagent. This corresponded with a fairly major increase in the price of the kits. But, after a few months, the kits were actually being shipped.

Anyway, it was also around this time, or maybe a little earlier, that the quality of the F5 reads really tanked. Not sure if things have improved recently with the 5500XL release. But given that, initially anyway, the new ECC chemistry is not even offered for PE reads, my guess is that the answer is "no".

My advice is that if you use the PE chemistry, treat the F5 reads as a minor "bonus" to an otherwise single read data set. Don't expect much more than 10-20% mapping rates. Then you won't be disappointed.

By the way, typically no "pass filter" removal of poor quality reads is done with SOLiD data sets. So looking at your raw data files won't tell you what is going on unless you use FastQC. Even that depends on the quality values generated by the instrument being accurate. However there is an error correction utility, "SAET" that can be run on a data set. We always run it. But it is pretty resource intensive. Nevertheless, it uses the embedded error detection ability of color space to improve the sequence accuracy. So, if you have the computational resources available, I would recommend it.

--
Phillip

[chenyao]

There is some history to this, but the sad truth is that at a least as of v4, SOLiD F5 reads are of much lower quality than F3 reads. It is not surprising that F5 reads map much more poorly than F3.

We immediately purchased the new paired end chemistry as soon as ABI released it for the SOLiD (with the release of v3 chemistry?) Unlike Illumina PE chemistry, no new template strand is created. Ligases are more or less 3'->5'/5'->3' blind. So conceptually this had always seemed like a strong point for the SOLiD. However, as with many things, the implementation proved more difficult than the concept would have suggested.

Nevertheless our first experiences with the new chemistry were very positive. Cost for the new PE kits was very reasonable. Mapping percentages were about what we thought we would expect for the shorter reads.

The only problem was that orders of the PE chemistry took a very long time to be filled. After many bitter complaints, I was told that an outside provider of a reagent critical for the the chemistry had failed to keep up with ABI's demand for this reagent. ABI had finally in-housed production of this reagent. This corresponded with a fairly major increase in the price of the kits. But, after a few months, the kits were actually being shipped.

Anyway, it was also around this time, or maybe a little earlier, that the quality of the F5 reads really tanked. Not sure if things have improved recently with the 5500XL release. But given that, initially anyway, the new ECC chemistry is not even offered for PE reads, my guess is that the answer is "no".

My advice is that if you use the PE chemistry, treat the F5 reads as a minor "bonus" to an otherwise single read data set. Don't expect much more than 10-20% mapping rates. Then you won't be disappointed.

By the way, typically no "pass filter" removal of poor quality reads is done with SOLiD data sets. So looking at your raw data files won't tell you what is going on unless you use FastQC. Even that depends on the quality values generated by the instrument being accurate. However there is an error correction utility, "SAET" that can be run on a data set. We always run it. But it is pretty resource intensive. Nevertheless, it uses the embedded error detection ability of color space to improve the sequence accuracy. So, if you have the computational resources available, I would recommend it.

--
Phillip

Thank you. the statistic of data provider also shows F5 has bad mapping quality. However, the "flagstat" shows read 2 maps better than read 1, but I put the F3 read as the first input. I don't know why. Does the read1 and read2 didn't follow the order of my input?

[pmiguel]

Thank you. the statistic of data provider also shows F5 has bad mapping quality. However, the "flagstat" shows read 2 maps better than read 1, but I put the F3 read as the first input. I don't know why. Does the read1 and read2 didn't follow the order of my input?

Possibly not. Standard protocol for the PE read is that F5 is done first, then F3. So, if that is the information stored in the bam file, that would explain your results.

Are you doing the mapping yourself? Or did your provider do the mapping?
I have never used Top hat. I am not sure how it assigns the read designations. Can you look at your bam file in a viewer like IGV? Might be possible to see some of the mapped reads and determine whether the long (F3) reads are called "read1" or "read2".

If you are getting high mapping rates from your F5 reads, then that is great news, in a way. It would indicate that ABI has fixed their issues with the F5 chemistry. The failure to map F3 reads might be the result of short inserts in your amplicons then. That is reads continuing into the adapter on the other side of the insert. Depending on the methodology that Top hat uses for mapping, this might prevent the read from being mapped.
--
Phillip