We are doing whole exome sequencing on trios: an affected child and its parents, where we are looking for de novo variants that the child has but not the parents. The expected number of these de novo variants for each trio is just 1-3. Most our data comes from HiSeq2000 sequencing, but our sequencing center recently bought a HiSeq2500 machine. We tested a sample on the new HiSeq2500 (only the kids, we don't have HiSeq2500 yet for the parents), and overall everything looks much better. Higher coverages, increased base quality etc. However, where we previously found 1 de novo variant, we now find 27.
The 26 extra de novo variants are all transversions (T-->G & A-->C). When looking in the BAM file, we can furthermore see that all these 26 extra variants have 100% strand bias. That is: all the T-->G variants are found on the forward strand, whereas all the A-->C variants are found on the reverse strand. Since A-->C is the reverse complement of T-->G, I think I can make the assumption that all these 26 variants feature the same change (T-->G).
What is even more striking is that when I then fetch a region of 20bp around the variants from Ensembl, I can find a very clear, albeit small, motif: TCAT, where the first T is the affected base. This is unlikely to be something of our pipeline, since it is found in the very earliest BAM file (aka directly post-alignment).
Is there anything known about extreme stand bias with Illumina HiSeq2500, or can anything be done to prevent it? (Sure, I could filter it out, but that's not a very nice way of handling this imho).
EDIT: the chemistry type used was V4 for the HiSeq2500, V3 for HiSeq2000.
EDIT2: After acquiring more data, from more runs, I find the exact same variants in different runs with the same type of strand biases. It's a reproducible artifact.
Above do you mean that all forward reads would indicate a SNP and all reverse reads over the same area would not. What is the coverage? Or do you mean that some, but only forward reads indicate the SNP and no reverse reads have the SNP.
I mean the latter. WHen there is a T-->G SNP, G's are only found on forward reads. Not all forward reads have the variant. When there is a A-->C SNP, C's are only found on the reverse reads.
Coverage is 26-250. It varies.
Have you checked the base quality plots for your reads? Maybe a sequencing cycle got bad for a part of your reads?
Of course. Those look very good, in fact.
This problem needs more clarifications. It almost feels like some sort of reporting or display option that got turned on.
A SNP of T->G is indeed an A->C if you were looking it from a reverse strand' view.
When we perform an alignment with the typical tools we "force" all reads to be displayed relative to the forward strand even those that aligned to the reverse strand.
If one were to display reads on their "original" and "true" alignment these would indeed be displayed as all forward reads having complementary bases relative to the reverse ones. We don't do that because it would quite confusing - but it could be used as a sanity check.
Excuse me:
Could you please tell me how you find a de novo variant? What software or scrips did you use and how you did it? Hoping for your reply.
Thx in advance!