to_paf.py codes the CIGAR incorrectly

[freeseek]

Generate a toy chain file between GRCh38 <-> T2T-CHM13v2.0:

$ echo -e "chain 2089 chr1 248956422 + 123985788 123986987 chr1 248387328 + 123418353 123419553 2482\n340 700 701\n159" > input.chain
$ cat input.chain
chain 2089 chr1 248956422 + 123985788 123986987 chr1 248387328 + 123418353 123419553 2482
340 700 701
159

Assuming I have GRCh38.fa and chm13v2.0.fasta, if I use chain2paf as follows I get:

$ chain2paf -i input.chain -f GRCh38.fa chm13v2.0.fasta
chr1	248387328	123418353	123419553	+	chr1	248956422	123985788	123986987	462	1200	255	cg:Z:5=1X18=1X36=1X9=1X12=1X15=1X38=1X18=2X1=1X47=1X24=1X6=1X44=1X4=1X46=1X1=701I700D2X2=1X4=1X1=1X8=1X22=1X18=2X7=1X2=1X3=1X2=1X6=2X17=1X12=1X13=1X2=2X11=1X8=

while if I use to_paf.py as follows I get:

$ to_paf.py -c input.chain -t GRCh38.fa -q chm13v2.0.fasta
chr1	248387328	123418353	123419553	+	chr1	248956422	123985788	123986987	9	709	255	cg:Z:5=1X18=1X36=1X9=1X12=1X15=1X38=1X18=2X1=1X47=1X24=1X6=1X44=1X4=1X46=1X1=700X1I2X2=1X4=1X1=1X8=1X22=1X18=2X7=1X2=1X3=1X2=1X6=2X17=1X12=1X13=1X2=2X11=1X8=	NM:i:738

The CIGAR outputs are actually different. If I want to highlight the CIGAR differences I get:

$ diff <(chain2paf -i input.chain -f GRCh38.fa chm13v2.0.fasta | cut -f13 | sed 's/\([MIDX=]\)/\1\n/g') <(to_paf.py -c input.chain -t GRCh38.fa -q chm13v2.0.fasta | cut -f13 | sed 's/\([MIDX=]\)/\1\n/g')
32,33c32,33
< 701I
< 700D
---
> 700X
> 1I

This is quite an issue because it causes chain information to be lost:

$ to_paf.py -c input.chain -t GRCh38.fa -q chm13v2.0.fasta > input.paf && paf2chain -i input.paf
chain	255	chr1	248956422	+	123985788	123986987	chr1	248387328	+	123418353	123419553	0
1040	0	1
159

which is different from the original chain file and it introduces a lot of matching bases which might not be matching at all.

Basically using to_paf.py all chain simultaneous gaps in both target and query are removed. This is actually what happened in the GRCh38 <-> T2T-CHM13v2.0 chains now distributed by UCSC and generated as explained here:

nextflow run main.nf --source GRCh38.fa --target chm13v2.0.fasta --outdir dir -profile local --aligner minimap2
python chaintools/src/split.py -c input.chain -o input-split.chain
python chaintools/src/to_paf.py -c input-split.chain -t GRCh38.fa -q chm13v2.0.fasta -o input-split.paf
awk '$1==$6' input-split.paf | rb break-paf --max-size 10000  | rb trim-paf -r | rb invert | rb trim-paf -r | rb invert > out.paf
paf2chain -i out.paf > out.chain
python chaintools/src/invert.py -c out.chain -o out_inverted.chain

This can be quickly verified as follows:

$ wget -O- http://hgdownload.cse.ucsc.edu/goldenPath/hs1/liftOver/hg38-chm13v2.over.chain.gz | zcat | awk 'NF==3 && $2>0 && $3>0' | wc -l
0

Because to_paf.py was used instead of chain2paf all the simultaneous chain gaps in both target and query were removed. I cannot fully understand what happens downstream as a result, but maybe the GRCh38 <-> T2T-CHM13v2.0 chain files should be regenerated as a result.

[milkschen]

Hi @freeseek ,

Thank you for reporting the issue with great details. This is highly appreciated. I will look into this.

[milkschen]

Hi,

chaintools would prefer 1M (whether match or mismatch) than a deletion plus an insertion (i.e., 1I+1D) at a locus. Thus, it is a design to not show a record like 340 700 701.

Therefore, for the grch38 <-> chm13v2 chains, it intended to avoid simultaneous chain gaps. We would be happy to update the chain files if there are clear examples showing improvements using another strategy.

However, I would agree that the default to_paf behavior does not fully translate the chain information. At least an improved doc would avoid the confusion you raised.

[freeseek]

Do you still have the original input.chain file generated by the nf-LO pipeline? I am worried that you might have a situation where a small inversion could have taken place, split.py would not have broken the chain at that locus and to_paf.py would have forced the bases within the inversion to match while another chain with the inverse orientation would have taken care of the small inversion instead. It would be some work to find such an example, but do you have a reason for why something like this wouldn't actually happen?

[freeseek]

Here is an example:

$ wget https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/chain/v1_nflo/grch38-chm13v2.chain
$ cat grch38-chm13v2.chain | awk 'NR>=255691 && NR<=255723'
chain	255	chr17	83257441	+	21975454	21992058	chr17	84276897	+	21947772	21960209	0
425	5	0
202	15	0
278	5	0
127	0	73
115	0	120
125	5	0
5	0	5
460	5	0
87	10	0
3219	0	10
27	5	0
438	10	0
391	10	0
46	0	10
150	15	0
466	0	5
1352	65	0
220	0	10
439	0	5
1335	4245	0
59	0	5
66	5	0
34	15	0
128	0	5
230	10	0
679	0	5
370	10	0
67	10	0
107	0	5
253	0	10
68	0	10
191
$ cat grch38-chm13v2.chain | awk 'NR>=255691 && NR<=255723' > input.chain
$ chain2paf -i input.chain -f hg38.fa hs1.fa
chr17	84276897	21947772	21960209	+	chr17	83257441	21975454	21992058	11131	16604	255	cg:Z:425=5D56=1X145=15D278=5D127=73I20X17=1X77=120I12=1X112=5D5=5I32=1X427=5D38=1X48=10D1583=1X155=1X62=1X44=1X297=1X321=1X751=10I27=5D438=10D8=1X80=1X67=1X214=1X18=10D46=10I143=1X6=15D39=1X205=1X116=1X19=1X83=5I1352=65D37X183=10I439=5I23=1X349=553X1=30X1=5X1=2X1=8X1=14X1=7X1=36X1=2X1=12X1=7X1=25X1=3X2=10X1=13X1=10X1=14X1=7X1=15X1=2X1=12X1=1X1=7X1=35X1=3X1=8X1=10X1=3X1=30X1=34X1=6X1=2X1=3X4245D55=1X3=5I9=1X32=1X15=1X7=5D29=1X4=15D31=2X22=2X5=1X24=1X15=1X21=1X2=5I7=1X1=1X13=1X6=1X199=10D38=1X10=1X568=1X60=5I271=1X98=10D67=10D107=5I253=10I68=10I191=
$ chain2paf -i input.chain -f hg38.fa hs1.fa | cut -f13 | sed 's/\([MIDX=]\)/\1\n/g' | awk 'NR>=80 && NR<=149'
5I
23=
1X
349=
553X
1=
30X
1=
5X
1=
2X
1=
8X
1=
14X
1=
7X
1=
36X
1=
2X
1=
12X
1=
7X
1=
25X
1=
3X
2=
10X
1=
13X
1=
10X
1=
14X
1=
7X
1=
15X
1=
2X
1=
12X
1=
1X
1=
7X
1=
35X
1=
3X
1=
8X
1=
10X
1=
3X
1=
30X
1=
34X
1=
6X
1=
2X
1=
3X
4245D

My best guess is that to_paf.py in this particular case replaced the chain gap:

373	5231	986

with the chain gap:

1359	4245	0

And it caused a huge amount of base pairs to mismatch for sequence that is completely unrelated

Of course, this problem is very sporadic, but it is probably widespread in the chain file and it is specific to the chain files processed with to_paf.py and it is not present in the chain files generated by UCSC (e.g. hg38ToHs1.over.chain.gz)

[milkschen]

Thanks for the informative example!

I added an option in to_paf.py (#17) to incorporate your suggestion.

In the CHM13-GRCh{37,38} chain, we used the reduction technique for chain records with a break point (simultaneous indel) <= 1000 base pairs. A chain record with a larger break points were split and considered independently.

[maximilianh]

Does this mean that we should update our chm13/hg38 chain files at UCSC ?

…

On Mon, Oct 23, 2023 at 06:58 Nae-Chyun Chen ***@***.***> wrote: Thanks for the informative example! I added an option in to_paf.py (#17 <#17>) to incorporate your suggestion. In the CHM13-GRCh{37,38} chain, we used the reduction technique for chain records with a break point (simultaneous indel) <= 1000 base pairs. A chain record with a larger break points were split and considered independently. — Reply to this email directly, view it on GitHub <#11 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AACL4TLCAUCV2VOHKQXXUPLYAX2OPAVCNFSM6AAAAAA5KMKDQSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONZUGQZTONJTGY> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[freeseek]

I do not understand why it would have to be an option with the original behavior being the default behavior. It is quite clear that this is a bug. The previous example causes the GRCh38 region:

chr17	21984176	21985511

which includes a gap in GRCh38, to match the region:

$ echo -e "chr17\t21984176\t21985511" | ./liftOver /dev/stdin hg38-chm13v2.over.chain.gz /dev/stdout /dev/stderr
Reading liftover chains
Mapping coordinates
chr17	21956582	21957917

Using the T2T chain file SNP rs1165854003 with GRCh38 position chr17:21985490 gets mapped to T2T position chr17:21957896. Now this makes no sense at all as the two regions have no homology whatsoever:

$ samtools faidx hg38.fa chr17:21985480-21985500
>chr17:21985480-21985500
GGAATGGAATGGAATGGAATG
$ samtools faidx hs1.fa chr17:21957886-21957906
>chr17:21957886-21957906
atggaatggatttgaatcaac

The behavior of to_paf.py is incorrect. It is not a matter of preference. And my take is that UCSC should not redistribute the T2T chain files until they are fixed