Gaps in alignment

Question

[mikegrootde]

I have aligned 100bp paired-end reads from 25x coverage sequencing experiment to a reference genome. When I view the alignment I sometimes see these large gaps where almost no reads align whatsoever:

Why does this happen and can I do something about it?

[Carlo Yague]

Hi, on what reference genome did you align the reads ? What kind of sequencing data is it ? How were the sequencing libraries prepared ?

[mikegrootde]

Hi Carlo,

Reference genome: Clint_PTRv2/panTro6 (Chimpanzee) https://hgdownload.soe.ucsc.edu/goldenPath/panTro6/bigZips/panTro6.fa.gz)

Sequencing: Illumina Hiseq 2000 instrument, PCR-based library https://www.nature.com/articles/nature12228#Sec2

Sample: https://www.ncbi.nlm.nih.gov/sra?LinkName=biosample_sra&from_uid=1920536

Alignment: bwa mem -T 0 (output no aligments score < 0)

[Carlo Yague]

One possibility is that the gaps correspond to DNA sequences of repetitive nature. Since you have excluded multi-mapping reads (whose alignment score = 0), repeated regions are expected to have 0 coverage. To validate this, you could either rerun bwa without the option -T, or investigate if the gaps overlap annotated repeated regions.

[mikegrootde]

Hi Carlo,

Thanks a lot for your suggestion!

From the bwa manual:

-T INT Don’t output alignment with score lower than INT. This option only affects output. [30]

With -T set to 0: If reads with score lower than 0 are omitted does't this mean reads with score exactly 0 are still in the output?

Also if I leave -T altogether won't bwa mem resort to default value [30] which would output no reads with a score < 30, leaving the same problem?

How do I allow multi-mapped reads in bwa mem?

It is very important for my downstream analysis to keep reads aligning to repetitive regions as I am doing variant calling on repeats.

[Carlo Yague]

You have no choice but allow multi-mapped reads if you are interested in repeated regions. Also make sure that you are using an unmasked reference genome. I can not tell you much more than that because I'm not an expert on analysis on CGG short tandem repeats. From what I read, the number of repeats vary from individual to individual, so it might requires a more specific mapping strategy to capture the reads supporting CGG repeats. See this paper for instance.

[mikegrootde]

Im still unsure how to allow multi-mapped reads with bwa mem, see my previous comment. (Sorry I changed it, still getting used to biostars site).

Another option is -a:

-a Output all found alignments for single-end or unpaired paired-end reads. These alignments will be flagged as secondary alignments.

Isn't that the way to do it?

Answer 1

[Istvan Albert]

First, you probably don't need to set the -T flag, I believe that to be the alignment score (the sum of reward + penalties).

If the score were zero, the alignment did not really work. I don't think -T =0 does anything. I think -T it is more of a way to select high scoring alignments rather than low-scoring ones. bwa won't even look for low-scoring alignments to begin with.

In general, you should aim to produce more/all alignments, then filter your BAM file later, rather than forcing the aligner to make various decisions. There is little to be gained in most cases.

The mapping quality filter is -q, (in samtools) would be a way to filter out multi mapped reads as bwa assigns a mapQ=0 for multimapped reads.

Now as for the gaps, it is a bit curious, but what region is this, show more information, what regions is this, also download some of the SNP/annotation data they have on the website, see if anything overlaps here

https://www.biologiaevolutiva.org/greatape/data.html

the human genome took many tries to complete, the chimp has nowhere that much effort put into it.

[mikegrootde]

Thank you for your reply. I believe the gaps are regions with CGG repeats and overall high GC content. I pasted the sequences corresponding to some gaps in a GC content calculator and get GC% of around 70.

Could it be possible that there are just no high GC% reads in this chimp data because of the bias in the PCR library amplification step used in sequencing?

[Istvan Albert]

Yes that is possible,

In addition assembling these regions is also much more difficult and error-prone.

Pick an additional newer data, or perhaps one with long-read technologies, and see what happens.

Answer 2

[mikegrootde]

I tried a PCR-free dataset and it works like a charm. We need more PCR-free WGS!! Thank you for your comments Istvan and Carlo