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$
3
  • 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, 2017 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$ Jul 17, 2017 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, 2017 at 18:29

4 Answers 4

6
$\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.

$\endgroup$
3
  • $\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$ Jul 20, 2017 at 12:54
  • 2
    $\begingroup$ @KamilSJaron think about 150X coverage. $\endgroup$
    – user172818
    Jul 20, 2017 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$ Jul 21, 2017 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.

$\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.

$\endgroup$
1
  • $\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$ Jul 19, 2017 at 15:58
1
$\begingroup$

You might want to check out this paper:

They have collected a list of "dark" (regions that have few reads aligning to them) and "camouflaged" (regions that have ambiguous alignment) regions of the human genome:

Here, we identify regions with few mappable reads that we call dark by depth, and others that have ambiguous alignment, called camouflaged.

[. . .]

Other dark regions arise, not because the sequencing is inherently problematic, but because of bioinformatic challenges. Specifically, many dark regions arise from duplicated genomic regions, where confidently aligning short reads to a unique location is not possible; we term these regions as “camouflaged”. These camouflaged regions are generally either large contiguous tandem repeats (e.g., centromeres, telomeres, and other short tandem repeats), or a larger specific DNA region that has been duplicated (e.g., a gene duplication) either in tandem or in a more distal genome region. In fact, many genes in the human genome were duplicated over evolutionary time and are still transcriptionally and translationally active (e.g., heat-shock proteins) [3–9], while others have been duplicated, but are considered inactive (i.e., pseudogenes). Regardless of whether the duplication is active, however, any genomic region that has been nearly identically duplicated and is large enough to prevent sequencing reads from aligning unambiguously will be “dark”, because the aligner cannot determine which genomic region the read originated from.

The authors have made their scripts available on github, along with bed files for the relevant regions. The answer to your question, i.e. a specific percentage, is hard to give since it varies enormously with the genome build and sequencing technology, but this paper should have all the details you want. Have a look at Table 1:

Dark and camouflaged regions vary by genome build. We identified dark and camouflaged regions throughout the genome for three different builds, including GRCh37, GRCh38, and GRCh38+alt, across five different sequencing technologies (or read lengths for Illumina). Specifically, we measured dark regions for Illumina based on 100-nucleotide read lengths, Illumina based on 250-nucleotide read lengths, 10x Genomics, PacBio, and Oxford Nanopore Technologies (ONT). Here, the counts for dark and camouflaged regions are combined. We found that the number of dark regions and nucleotides, both within gene bodies (represented as GB in the table) and outside gene bodies, varies dramatically by build and technology. Overall, each technology has its respective strengths. GRCh38 including alternate contigs has > 3x more dark nucleotides than GRCh37, and more than 2x more dark regions. Results presented throughout the manuscript are based on GRCh38 (in gray)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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