-
Notifications
You must be signed in to change notification settings - Fork 2
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
to_paf.py codes the CIGAR incorrectly #11
Comments
Hi @freeseek , Thank you for reporting the issue with great details. This is highly appreciated. I will look into this. |
Hi, chaintools would prefer 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 |
Do you still have the original |
Here is an example:
My best guess is that
with the chain gap:
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 |
Thanks for the informative example! I added an option in 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. |
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: ***@***.***>
|
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:
which includes a gap in GRCh38, to match the region:
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:
The behavior of |
Generate a toy chain file between GRCh38 <-> T2T-CHM13v2.0:
Assuming I have
GRCh38.fa
andchm13v2.0.fasta
, if I usechain2paf
as follows I get:while if I use
to_paf.py
as follows I get:The CIGAR outputs are actually different. If I want to highlight the CIGAR differences I get:
This is quite an issue because it causes chain information to be lost:
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:This can be quickly verified as follows:
Because
to_paf.py
was used instead ofchain2paf
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.The text was updated successfully, but these errors were encountered: