Reference bias occurs during read mapping and causes the mapped reads to more closely resemble the reference genome than it actually does. This happens because alignment programs inherently favor alignments with a greater number of matches and with read mapping we are aligning sequences that do not match. Conceptually, this is well defined. But other than Sarver et al. 2017, little has been done to directly quantify the effects of reference bias. However, we can layout some expectations:

1. Reference bias will lead to lower divergence estimates between the reference species and the mapped species.
2. When mapping multiple species to a single reference, all species will appear to be more similar with respect to sequence.
3. Reference bias can lead to false negatives in downstream analyses.

Briefly, reference bias occurs because some reads remain unmapped during the read mapping process. This can happen because they contain many sequencing errors (red below) OR because they contain many actual variants (green below):

Figure 5.1: Unmapped reads can be reads containing errors (dark red) and/or variants (green).

All of those unmapped reads with variants lead to our final consensus sequence looking more like the reference than it should.

Pseudo-it is software designed to mitigate the effects of reference bias by using an iterative mapping approach to generate pseudo-reference assemblies. It performs the following steps:

1. Map reads to a reference genome.
2. Call variants based on the mapped reads.
3. Generate a consensus sequence using the called variants.
4. Repeat the process, but instead of mapping reads to the reference, map them to the newly generated consensus sequence.

The idea here is to generate a sequence that incorporates the variation that has been mapped to allow mapping of the reads that were previously too divergent:

Figure 5.2: Iterative mapping incorporates mapped variation to allow for more reads to be mapped.

The first version of this software (Sarver et al. 2017) used BWA to map reads and GATK 3. to call variants. It was shown to increase estimated divergence in a set of mouse exomes:

Figure 5.3: Divergence increases (A) and bias decreases (B) over iterations of mapping in several mouse species when mapping to the same reference genome (Sarver et al. 2017).

The newest version that we'll be using today will still use BWA, but is upgraded to use GATK 4 and incorporates indels in the final pseudo-assembly. Let's try it out!

First, let's navigate back to the main project directory:

cd ..

And now we can run psuedo-it on D. pseudoobscura chromosome 2 with two iterations:

pseudo_it -ref dmel-3R-reference/Drosophila_melanogaster.BDGP6.28.dna.chromosome.3R.fa -pe1 dpse-chr2-reads/illumina/pse-chr2_1.fastq.gz -pe2 dpse-chr2-reads/illumina/pse-chr2_2.fastq.gz -i 2 -p 4 -o dpse-pseudo-it

This will perform the steps outlined above twice (mapping, variant calling, consensus) and will save the relevant files for each iteration, including a final consensus FASTA file as a pseudo-assembly. You can explore the output folder dpse-pesudo-it to see what's there. You can view its progress on the screen as it runs.

Note that this runs an entire gamut of mapping, variant calling, and consensus sequence generation, which can take a very long time (about 1 week for a mammalian whole genome with 4 iterations and 4 cores per chromosome). Even for our single chromosome from D. pseudoobscura this would take over an hour. So let's go ahead and cancel this run by hitting <ctrl>+c while the Terminal is selected.

In order to see how pseudo-it has reduced reference bias, we would first need to generate the consensus sequence after one iteration:

pseudo_it -ref dmel-3R-reference/Drosophila_melanogaster.BDGP6.28.dna.chromosome.3R.fa -pe1 dpse-chr2-reads/illumina/pse-chr2_1.fastq.gz -pe2 dpse-chr2-reads/illumina/pse-chr2_2.fastq.gz -i 1 -p 4 -resume dpse-pseudo-it

However, since this would also take some time, we'll just move forward with the provided expected output in dpse-pseudo-it.

Remember, reference bias makes mapped reads look more similar to the reference that expected, so we can outline a couple of expectations to test:

1. More reads should map to the second iteration.
2. The consensus sequence from the second iteration should appear to be more diverged.

01. More reads should map to the second iteration

This is fairly easy to check with samtools flagstat. First, we can run flagstat on the BAM file from the first iteration:

samtools flagstat expected-outputs/dpse-pseudo-it/iter-01/bam/merged-rg-mkdup-iter-01.bam.gz

This results in the following output:

404438 + 0 mapped (37.53% : N/A)

Now let's check the second iteration:

samtools flagstat expected-outputs/dpse-pseudo-it/iter-02/bam/merged-rg-mkdup-iter-02.bam.gz

420425 + 0 mapped (38.97% : N/A)

A modest increase in this example, but an increase nonetheless.

Another effect of more reads mapping is that we should have higher average read depth per iteration. Recall that with our original mappings (one iteration), we calculated an average read depth on this chromosome of 3.95, even though our expected coverage based on the Lander/Waterman equation was about 5. Let's see if our iterative mapping has improved read depth.

Then, like before, let's run samtools depth on our BAM file:

samtools depth expected-outputs/dpse-pseudo-it/iter-02/bam/merged-rg-mkdup-iter-02.bam.gz > dpse2-to-dmel3R-pi2-depth.tab

And let's calculate average depth with our custom R script:

Rscript scripts/depth_plot.r dpse2-to-dmel3R-pi2-depth.tab

This will generate a read depth histogram (depth-hist.png) and tell us the average read depth:

# Avg. read depth: 4.174893

Again, a modest increase from 3.95, but a step in the right direction, showing that iterative mapping improves read depth.

02. The consensus sequence from the second iteration should appear to be more diverged.

There are many ways to calculate divergence. We'll use the simplest one: the percent identity, or the number of sites that match between the alignment of two sequences. To do that we'll first need to generate some alignments.

To generate alignments between these two chromosomes, we'll actually be using a long-read mapper called minimap2. I've found this program useful for generating between species chromosome alignments.

Next, let's use minimap2 to generate an alignment between the first iteration of mapping and the original D. melanogaster 3R reference:

minimap2 -x asm10 dmel-3R-reference/Drosophila_melanogaster.BDGP6.28.dna.chromosome.3R.fa expected-outputs/dpse-pseudo-it/iter-01/fa/iter-01-final.fa > alignments/dpse-2-pi-1-to-dmel3.paf

minimap2 output can be formatted in two ways, the BAM format we are familiar with or PAF, the pairwise alignment format. I believe this format is unique to minimap2, but is simply a tab-delimited file that contains informationo about blocks that have been aligned, such as sequence ID, start position, number of matching bases, etc.

Let's take a look at this particular PAF file:

less -S alignments/dpse-2-pi-1-to-dmel3.paf

3R      32168558        50522   32168555        +       3R      32079331        50522   32079328        30591593        32177084        60      tp:A:P  cm:i:3034393    s1:i:30543658   s2:i:25089      dv:f:0.0028 rl:i:147655
3R      32168558        5       44485   +       3R      32079331        5       44485   23645   44480   60      tp:A:P  cm:i:2055       s1:i:23645      s2:i:3382       dv:f:0.0000     rl:i:147655

This shows that between the D. pesudoobscura pseudo-reference of chromosome 2 and the D. melanogaster reference of chromosome 3R there are two large blocks that could be confidently aligned (one block per line).

The following table explains the PAF columns explicitly:

Figure 5.4: The column definitions of PAF files. Source.

Conveniently, this tells us the number of matches (column 10) and total number of sites (column 11) per alignment block. We can use this to quickly calculate the percent identity from a PAF file. We've provided a custom python script, div_est.py that does this. Let's run this script on the alignments/dpse-2-pi-1-to-dmel3.paf file.

python3 scripts/div_est.py alignments/dpse-2-pi-1-to-dmel3.paf

This should simply print the percent identity to the screen:

# Percent matching bases: 95.015

95.015% identity after one iteration of mapping. What about after two iterations? Remember, reference bias makes mapped reads look more like the reference than expected, so we would expect % identity to DECREASE if we are improving on reference bias.

Next, run minimap2 on the second iteration of mapping:

minimap2 -x asm10 dmel-3R-reference/Drosophila_melanogaster.BDGP6.28.dna.chromosome.3R.fa expected-outputs/dpse-pseudo-it/iter-02/fa/iter-02-final.fa > alignments/dpse-2-pi-2-to-dmel3.paf

And then again let's calculate % identity:

python3 scripts/div_est.py alignments/dpse-2-pi-2-to-dmel3.paf

# Percent matching bases: 94.541

Compared to the previous 95.015%, the % identity has decreased with an iteration of mapping as expected. This means we are likely capturing more variation in this version of the pseudo-reference.

This brings us to the end of our workshop. Thanks for attending!