4
$\begingroup$

When doing Illumina 2x150bp sequencing of genomic DNA, and after aligning the reads to GRCh38, what percentage of the non-N fraction of the human genome is MAPQ=0? This is, what part corresponds to regions that can't be uniquely mapped with 2x150bp reads.

And, how many genes are affected by MAPQ=0 regions?

I presume the numbers will dance around depending on fragment sizes, read quality, etc. but I am happy with some starting numbers.

$\endgroup$
  • 1
    $\begingroup$ Are you expecting someone to have these numbers handy or do you want to know how to figure this out? :) $\endgroup$ – Devon Ryan Jul 17 '17 at 12:27
  • 2
    $\begingroup$ UCSC used to have downloadable mappability tracks. These don’t seem to exist any more; the best bet now is to generate them yourself, using GEM. $\endgroup$ – Konrad Rudolph Jul 17 '17 at 12:36
  • 1
    $\begingroup$ You can start with the repeat masked regions of the genome. That would reduce the search space if you want to look at 300bp regions. What is the insert size of your 2x150 library? That will play a large role in the mappability of the pair. $\endgroup$ – Bioathlete Jul 18 '17 at 18:29
5
$\begingroup$

Mapping quality is determined by the repetitiveness of the genome, the sequencing error rate, insert size, the capability of the mapper and the nasty heuristics behind the mapper. MAPQ=0 to one mapper is not necessarily MAPQ=0 to another.

That said, I get what you mean. You want to know the uniqueness/repetitiveness. It is still hard if you want to get a useful answer. For 150bp reads, each reference position is covered by 150 reads. What if 50 of them have no other exact hits elsewhere, but the rest 100 have? Is this position a repeat or not? In addition, what if the 50 are unique only because one mismatch? If there is a variant at that mismatch, the 50 would become repeats or mapped elsewhere.

My preference is to say a position is "unique" under k-long reads if over k/2 reads overlapping the position have no other perfect or 1-mismatch/1-gap hits elsewhere in the genome. Under this definition, 79.3% of human genome are unique for 35bp reads. 92.4% for 75bp reads. I don't have the number for 150bp reads. I guess it will be around 95%. Empirically, 94-95% of human genome is callable with 100bp paired-end reads.

As to other measurements, the most common one is the fraction of reads that has an exact hit elsewhere. You can use Fabio's method. It gives a good enough estimate once you simulate over 1 million reads. This fraction is around 86% for 35bp reads and 95% for 75bp reads, as I remember. The problem with this approach is this fraction is not very informative to variant calling due to the issues I talked about. Another way is to use RepeatMasker. It is worse. RepeatMasker masks 50% of human genome, but excludes segmental duplications where short reads can't be confidently mapped.

|improve this answer|||||
$\endgroup$
  • $\begingroup$ "For 150bp reads, each reference position is covered by 150 reads." sounds wrong to me. This is what I think should be in the sentence "For 150bp reads, each read mapped to reference has a 150 positions." $\endgroup$ – Kamil S Jaron Jul 20 '17 at 12:54
  • 2
    $\begingroup$ @KamilSJaron think about 150X coverage. $\endgroup$ – user172818 Jul 20 '17 at 21:49
  • $\begingroup$ Ah, all right, now I get it. From the sentence, I thought that this supposed to be a property of 150bp reads. $\endgroup$ – Kamil S Jaron Jul 21 '17 at 9:53
5
$\begingroup$

As has been said before, mappability to the 'human genome' depends on a number of factors, among these the reference version and type of reads, for which you are interested in GRCh38 and 2x150bp reads. Although I am not aware of numbers accounting for these particular reference and type of Illumina reads, the 1000 genomes project has provided the community with a similar and close estimation that you might be interested in considering regarding your inquiry.

Similar to your question, the 1K project estimated 'the proportion of the human genome that is less accessible to short reads'. In these estimates the human genome is GRCh37 and the types of reads in question are mostly 2x Illumina with a mixture of lengths with the longest being up to 250bp. In these estimates each base in the human genome is considered (and marked) 'less accesible' according to these criteria:

L - depth of coverage is much lower than average

H - depth of coverage is much higher than average

Z - too many reads with zero mapping quality overlap this position

Q - the average mapping quality at the position is too low

Each of these criteria has "standard" and "strict" thresholds for a base to be considered - or not - in each category. You can read more in the link below:

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/accessible_genome_masks/README.accessible_genome_mask.20140520

According to the strict thresholds, the human genome has about 16.8% of "Z" and 3.1% of "Q" bases, respectively. Considering the "Z" and "Q" criteria as a proxy for ~ mapq=0, about 19.9% of the human genome can not be uniquely mapped.

|improve this answer|||||
$\endgroup$
1
$\begingroup$

I would suggest you to run some simulation. You can use wgsim to simulate reads that are as much similar to what you want an answer for. I don't think there are faster methods.

|improve this answer|||||
$\endgroup$
  • $\begingroup$ I don’t think this is an appropriate approach to investigate mappability because the sample will be too small. Tools like GEM perform exhaustive calculations. $\endgroup$ – Konrad Rudolph Jul 19 '17 at 15:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.