Pair end merging and statistical output

Question

[jomo018]

1. I am looking for a paired end merging utility similar to FLASH or PEAR that also outputs statistical information, mainly number of pair mismatches.

2. Is this type of information available from a SAM file after alignment with BOWTIE, BWA or some other aligner?

[Titus]

Hi, I think your information about number of pair mismatches (if you speak about your reads merged) is the consequences of your alignment then step 2 :) You can simply got the information if you compare the number of reads aligned between merged reads and no merged reads.

[jomo018]

I am looking for mismatches in base resolution, not read resolution. Pairs can be merged (or aligned concordantly) even if some overlapping bases disagree. I am looking for the number or rate of these mismatching bases.

[Titus]

Do you mean variant calling ?

[jomo018]

You can call it variant calling where one mate declares different base than the other.

Answer 1

[Brian Bushnell]

BBMerge can do this, with the "ecco" flag, which rather than combining the reads, just does error-correction via overlap:

bbmerge.sh in=reads.fq ecco mix out=corrected.fq
Total time: 1.890 seconds.

Pairs:                  1000000
Joined:                 182539          18.254%
Ambiguous:              817461          81.746%
No Solution:            0               0.000%
Too Short:              0               0.000%
Errors Corrected:       6994

Avg Insert:             159.9
Standard Deviation:     21.5
Mode:                   187

Insert range:           100 - 191
90th percentile:        186
75th percentile:        178
50th percentile:        164
25th percentile:        145
10th percentile:        128

[GenoMax]

In this example reads must be interleaved. @Brian: Is that a requirement?

[jomo018]

I am just reading the manual. They also allow in1 and in2 paired inputs.

[Brian Bushnell]

Yep, the syntax can also be:

bbmerge.sh in1=r1.fq in2=r2.fq ecco mix out1=corrected1.fq out2=correct2.fq

...but I normally show the interleaved version of the command for conciseness.

[jomo018]

Errors Corrected are number of pairs corrected rather than number of bases corrected. Right?

[Brian Bushnell]

No, it is the total number of bases corrected.

[jomo018]

Thank you Brian. Can you clarify the tag trimq=xx (as opposed to qtrim...). For example, suppose you have a low quality base in the middle of a read, is it considered an N or something else?

[Brian Bushnell]

"qtrim" tells the program which end to trim, while "trimq" specifies the quality threshold. Only ends can be trimmed (and mainly the right end is the important one for trimming with respect to merging). "qtrim=r trimq=15" will trim the bases on the right end such that the region trimmed has an average quality below 15, while the remaining region has an average quality at least 15. For example, if the last 5 base qualities were "20, 0, 17, 19, 16", then the last 4 bases would be trimmed as that region has an average quality below 15 (not that 0 is an N) but the 20 would not be trimmed. It's a little hard to calculate by eye because the scores are first transformed to the probability scale (where 16 is roughly 2.5% chance of error) before being averaged. Low-quality bases are not considered N. Rather, if two reads mismatch at a location and one is higher quality than the other, the base with the higher quality is assumed correct and the resuting quality score is higher-lower.