Read deduplication with extremely high coverage and variant calling

[Mike2188]

Hello,

I'm currently in a class, and the professor is having us run through some NGS practice data from HIV. I tried to find an exact answer in the forums but couldn't. I feel like the professor may be wrong with what he is asking us to do.

He gave us fastq files, and told us to clean them up, align them to the genome, remove duplicates with picard tools, and then do variant calling. The problem is, the coverage is insanely high (~200,000x) so every position has many many mapping reads.

If I understand correctly, Picard tools will select one read (is it the best mapping quality read, or most highly represented read?) at each mapping location and discard the rest.

I feel like this strategy isn't right. If the program is selecting for the best mapping read, I feel like this would select for variants that match the genome reference at that location even if they are a minor variant. If the program is using the most abundant read, I feel like every minor variant would definitely be lost. Am I wrong in my thinking??

[HESmith]

1) Are you sure you have 200K-fold genome coverage, or 200K reads?

2) Picard ranks duplicates by summed base-quality scores (see the documentation).

3) Reads are randomly distributed across the genome, so multiple reads will span every nucleotide. So even minor variants will be represented by a subset of reads (although the default parameters of your variant caller might need to be modified to detect low-frequency variants).

[Mike2188]

Yes, coverage is extremely high, had several million reads from a miseq run covering a small viral genome of ~10kb.

Oh, okay. If it takes the overall highest quality read, and doesn't consider mapping quality or the number of instances of each read, then each position will be sampled multiple times and variant calling will work. Reads were a long length (250 bases, SE) so each position would be sampled randomly ~250 times.

That makes sense to me. Thanks for clarifying.