Unexpected liftover behavior when lifting PacBio long genomic reads between two haplotypes

[minkinaa]

Hi Nae-Chyun,

I have been using leviosam2 to perform liftover between two haplotypes in the GM12878 cell line using PacBio long read data, and I have run into a few cases of unexpected liftover behavior I'll describe below. Though by and large liftover appears to happen correctly, some liftover events are not entirely consistent with the underlying chain file with various outcomes, ranging from small local misalignment to those impacting an entire read or set of reads at a particular location in the genome. I am unsure whether these are due to the same underlying issue. The optional realignment module unfortunately is not a great solution to these problems and I'll explain why below as well.

For the purpose of these examples, we are lifting reads which natively map to haplotype 2 (I'll call this the native haplotype) to haplotype 1's coordinates (I'll call this "target" haplotype). No realignment was performed. I did not run the full workflow -- just 'leviosam index' followed by 'leviosam lift'.

I have included the parameters & test files which can be used to recreate each issue (files linked at the bottom of post).

(1) Indels are not lining up after liftover

In the top panel in the example above, the top track contains lifted over reads (hap2 to hap1) & the bottom track shows the genomic alignment of hap2 onto hap1 from which the chain file was made. We can see that the chain contains a 3 & 291 base insertion within this region in hap2 relative to hap1. We see these insertions in the reads as well, but their position does not line up with the chain. This behavior seems unexpected -- shouldn't the liftover follow the chain exactly? (As a sanity check, the bottom panel shows these same reads aligned to their native genome, with the inverted chain at the bottom showing where the 291bp indel lies. There is no obvious variation in the native alignment in these reads -- it seems to be introduced in the liftover step).

Read issue 1 files & params for testing:
chain: hap2-to-hap1coords.inverted.chain
fasta: hap1.fa
fasta index: hap1.fa.fai
input_bam: reads_issue_1.bam
-G 100000
-f $fasta_path
-m
-S mapq:0
-T 1000
-S aln_score:1000

(2) Liftover appears unstable at homopolymer regions

Homopolymers consistently seem to cause liftover issues as shown in the top panel. Bases adjacent to the homopolymer are shifted in a subset of reads, but the misalignment often resolves with an indel at seemingly a random position. Interestingly, this behavior at homopolymers seems read-dependent, so if a longer misalignment happens at homopolymer X in read 1, we observe longer misalignments at other homopolymers in that read. See the boxed region of the bottom panel for examples -- for an individual read, the misalignments at the four homopolymer regions shown are all of consistent length.

Read issue 2 files & params for testing:
chain: hap2-to-hap1coords.inverted.chain
fasta: hap1.fa
fasta index: hap1.fa.fai
input_bam: reads_issue_2.bam
-G 100000
-f $fasta_path -m
-S mapq:0 -T 1000
-S aln_score:1000

(3) Entire reads get mis-lifted when their start position falls within a chain indel

There is an issue that arises when a read begins in a location in hap2 (native haplotype) which does not exist in hap1, either appearing as an insertion in the chain file (bottom panel), or a break between two parts of the chain (top panel). Unsurprisingly, if we look at the distribution of these reads across a single chromosome (bottom right), we observe peaks corresponding to these indel/chain break events. These reads are lifted over to the generally correct location (see top right inset), but don't align perfectly to the reference.

Depending on the chain used, thus far I've seen this to impact .1-.3% of reads, except in cases like the one below.

Read issue 3 files & params for testing:
chain: hap2-to-hap1coords.inverted.chain
fasta: hap1.fa
fasta index: hap1.fa.fai
input_bam: reads_issue_3.bam
-G 100000
-f $fasta_path -m
-S mapq:0 -T 1000
-S aln_score:1000

(4) Gross mis-lifting impacting all reads on one chromosome arm

Shown above is an example where almost all the reads on the long arm of chromosome 6 were lifted over improperly. I suspect that something about the complex chain structure at the centromere is causing every read that lifts over to a location beyond it to be mis-lifted. It does not appear to be a simple off by X error, as the reads are not all shifted the same number of bases from their proper alignment (though they tend to be shifted just a small number of bases from their correct location for the most part). This does seem to be an artifact specific to this chain. For example, we can lift hap2 reads to the hap1 coordinate system using a chain which was constructed by aligning the hap2 assembly to the hap1 assembly (and then inverting, chain 2 above). Or, we can instead align the hap1 assembly to hap2 (chain 1 above). In this case, chain 2 produced good results across chromosome 6, while only the short arm of chr6 has accurate liftover events with chain 1. I have verified that the chain files themselves are accurate.

Read issue 4 files & params for testing:
chain: hap1-to-hap2coords.chain (the liftover works ok for me when run with the other chain I provided here)
fasta: hap1.fa
fasta index: hap1.fa.fai
input_bam: reads_issue_4_good_read.bam (reads on the short arm of chr6 which lift over ok); reads_issue_4_bad_read.bam (reads on the long arm of chr6 which are being mislifted)
-G 100000
-f $fasta_path -m
-S mapq:0 -T 1000
-S aln_score:1000

Why adding the optional realignment module does not work

In order for reads to get realigned via the -x parameter, we have to set a parameter which would funnel these problematic reads into the deferred category. I did this by setting the max NM tag value to 200. While a subset of the reads were realigned as expected, this introduced a new issue when lifting over from one haplotype to another, where we expect some variation between the haplotypes. For example, if one haplotype contains a 300bp insertion relative to the other and a read overlaps this insertion, this read will be funneled into the deferred category. Thus, using this parameter with long read data, we lose a lot of perfectly fine reads.

Parameters used for leviosam2 lift with realignment:
-f ${fasta}
-G 100000
-m
-T 1000
-S mapq:0
-S aln_score:1000
-S hdist:200
-x config/pacbio_all.yaml

Referenced input files:
Remove the '.txt' extension from these to run. These were added to get around github format constraints.

• the .fa file is too large to upload here but I can provide via upon request. It is not actually necessary for reproducing the errors above

reads_issue_4_good_read.bam.txt
reads_issue_4_good_read.bam.bai.txt
reads_issue_4_bad_read.bam.txt
reads_issue_4_bad_read.bam.bai.txt
reads_issue_3.bam.txt
reads_issue_3.bam.bai.txt
reads_issue_2_homopolymer.bam.txt
reads_issue_2_homopolymer.bam.bai.txt
reads_issue_1.bam.txt
reads_issue_1.bam.bai.txt
hap2-to-hap1coords.inverted.chain.txt
hap1.fa.fai.txt
hap1-to-hap2coords.chain.txt

[minkinaa]

A few other minor issues & comments:

• When a read falls entirely within a region that does not exist in the target haplotype, I get an error that the cigar string does not match the read length for that read when I try to sort or run any samtools command on the output from leviosam lift. Removing supplementary alignments before lifting over solved this in my case (since these are shorter and more likely to exhibit this property), but I am not sure if this is exclusively an issue impacting supplementary alignments.
• If the original .bam file of reads already contains an NM and MD tag, leviosam lift does not rewrite, but adds a second NM and MD tag to the read line. I imagine this may cause issues downstream, and I am curious whether the optional realignment step is impacted by this, as it uses the NM tag.
• For PacBio long reads & lifting between haplotypes, we are finding -G 100000 to work well, as we don't want to remove reads which overlap a large indel which exists in one haplotype but not the other. mapq=0 is also useful if the original alignment was done to the diploid genome, as reads are likely to map reasonably well to two position in that genome and thus have a low mapq score.

[minkinaa]

A small update on my end about issue 4, with widespread misalignment impacting a large portion of a chromosome: I strongly suspect this is occurring because of overlaps in the chain file, but it only appears to be an issue when there are overlaps on the target side of the paf file used to make the chain (columns 6-8 in the paf format). I use the paf2chain tool to create the chain files, with some trimming to remove overlaps but only on the query side. Overlaps on the query side of the paf file do not result in this issue. Note, the terms target & query in this context are flipped relative to what I think of as the target & query during liftover. For example, if we are lifting hap2 (query) to hap1 (target) coordinates, the paf/chain files used to produce these have hap1 information in the query positions in the paf format, and hap2 information in the target positions. Overlaps in the chain in the target fields of the file (hap2 in this case) seem to cause this widespread misalignment issue.

[milkschen]

Thank you @minkinaa for sharing these issues! They might be related to the core implementation and I'll need some time to get to the bottom of the issues.

[milkschen]

I added this commit to address the overlapping chain file issue #92 The leviosam2 data structure does not handle cases where two intervals on the source genome, the commit validates the input chain and errors out when such conditions exist. It should avoid a lot of confusion for users.

[minkinaa]

Thanks, Nae-Chyun! We transitioned to a graph-based approach for our liftovers, but hopefully this helps future users!