open2c / pairtools

Extract 3D contacts (.pairs) from sequencing alignments
MIT License
104 stars 32 forks source link

Pairtools dedup ValueError or incompatible reads? #177

Closed ddepierre closed 1 year ago

ddepierre commented 1 year ago

Hi,

I am trying to run dedup command, but the command is running super long, I get error msgs and then my command is stopped/killed.

Command: pairtools dedup --output-stats 04_PAIRTOOLS/stats/sample.04_dedup.stats --max-mismatch 500 04_PAIRTOOLS/sam/sample.03_selected.temp.sam --output 04_PAIRTOOLS/sam/sample.04_dedup.temp.sam with sample.03_selected.temp.sam being a sam/pairs output from pairtools select

Error:

File "/pairtools/pairtools-1.0.2_venv/lib/python3.9/site-packages/pairtools/lib/dedup.py", line 267, in _dedup_chunk
    df_mapped.loc[:, "clusterid"] = connected_components(a_mat, directed=False)[1]
ValueError: Buffer dtype mismatch, expected 'ITYPE_t' but got 'long'
ValueError: Buffer dtype mismatch, expected 'ITYPE_t' but got 'long'
Exception ignored in: 'scipy.sparse.csgraph._traversal._connected_components_undirected'

I tried the same command line on other samples with about the same number of pairs, it works well in less than 10min, while on some other samples, it takes >10hours before I get the value type error.

Do you think it could come from the reads themselves? Here are some pairs.sam on which it runs fast:


SRR8811899.24702503 chr1    3005633 chr1    3005818 +   -   UU  SRR8811899.24702503161chr1300563360151M=3005668186AGTCAGAGCTATCAATGAGACATTCAAATCATGCAGCCTATTGCTATTGCCTTTGGTTGTTCCCAAGAGGCCAAAGGTAAATTCCTATTCATAAAAATACCATGCACTTTAGACACAAAAGTCAGAGTCCCCTGAGCTAAGGGGGAAAAAA???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????NM:i:0MD:Z:151MC:Z:151MAS:i:151XS:i:20Yt:Z:UU   SRR8811899.2470250381chr1300566860151M=3005633-186GCCTATTGCTATTGCCTTTGGTTGGTCCCAAGAGGCCAAAGGTAAATTCCTATTCATAAAAATACCATGCACTTTAGACACAAAAGTCAGAGTCCCCTGAGCTAAGGGGGAAAAAAAGAGCCTACTTTTTGAGGTACCAGAAGGTGTCATG???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????NM:i:1MD:Z:24T126MC:Z:151MAS:i:146XS:i:19Yt:Z:UU
SRR8811899.25133564 chr1    3005633 chr1    3005818 +   -   UU  SRR8811899.25133564161chr1300563360151M=3005668186AGTCAGAGCTATCAATGAGACATTCAAATCATGCAGCCTATTGCTATTGCCTTAGGTTGTTCCCAAGAGGCCAAAGGTAAATTCCTATTCATAAAAATACAATGCACTTTAGACACAAAAGTCAGAGTCACCTGAGCTAAGGGGGAAAAAA???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????NM:i:3MD:Z:53T46C28C21MC:Z:151MAS:i:136XS:i:19Yt:Z:UU   SRR8811899.2513356481chr1300566860151M=3005633-186GCCTATTGCTATTGCCTTTGGTTGTTCCCAAGAGGCCAAAGGTAAATTCCTATTCATAAAAATACCATGCACTTTAGACACAAAAGTCAGAGTCCCCTGAGCTAAGGGGGAAAAAAAGAGCCTACTTTTTGAGGTACCAGTAGGTGTCATG???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????NM:i:1MD:Z:140A10MC:Z:151MAS:i:146XS:i:22Yt:Z:UU
SRR8811899.25598630 chr1    3005633 chr1    3005818 +   -   UU  SRR8811899.25598630161chr1300563360151M=3005670186ATTCAGAGCTATAAATGAGACATTCAAATCATGCAGCCTATTGCTATTGCCTTTGGTTGTTCCCAAGAGGCAAAAGGTAAATTCCTATTCATAAAAATACCATGAACTTTAGACACAAAAGTCAGAGTCCACTGAGCTAAGGTGGAAAAAA???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????NM:i:6MD:Z:1G10C58C32C25C11G8MC:Z:2S149MAS:i:124XS:i:23Yt:Z:UU  SRR8811899.2559863081chr13005670602S149M=3005633-186TACTATTGATATTGCCTTTGGTTGTTCCCAAGAGTCCAAAGGTAAATTCCTATTCATAAAAATACCATGCACTTTAGACACAAAAGTCAGAGTCCCCTGAGCTAAGGGGGAAAAAAAGAGCCTACTTTTTGAGGTACCAGAAGGTGTCATG???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????NM:i:2MD:Z:6C25G116MC:Z:151MAS:i:139XS:i:19Yt:Z:UU
SRR8811899.25639144 chr1    3005633 chr1    3005818 +   -   UU  SRR8811899.25639144161chr1300563360151M=3005668186AGTCAGAGCTATCAATGAGACATTCAAATCATGCAGCCTATTGCTATTGCCTTTGGTTGTTCCCAAGAGGCCAAAGGTAAATTCCTATTCATAAAAATACCATGCACATTAGACACAAAAGACAGAGTCCCCTGAGCAAAGGGGCAAAAAA???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????NM:i:4MD:Z:107T13T15T6G6MC:Z:151MAS:i:131XS:i:23Yt:Z:UU SRR8811899.2563914481chr1300566860151M=3005633-186GCCTATTGATATTGCCTTTGGTTGTTCCCAAGAGGCCAAAGGTAAATTCCTATTCATAAAAATACCATGCACTTTAGACACAAAAGTCAGAGTCCCCTGAGCTAAGGGGGAAAAAAAGAGCCTACTTTTTGAGGTACCAGAAGGTGTCATG???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????NM:i:1MD:Z:8C142MC:Z:151MAS:i:146XS:i:22Yt:Z:UU

Here are some pairs.sam on which it runs super slow and then outputs an error:


A01878:56:HKKK3DMXY:2:2411:18909:16423  chr1    3005837 chr1    3005945 +   -   UU  A01878:56:HKKK3DMXY:2:2411:18909:16423161chr1300583760110M=3005836109AAGCCACCACAATCCAATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGTAACTATGATTGGACAFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFF:F:FFFFFFFFFFFFFFFFFFNM:i:1MD:Z:109T0MC:Z:110MAS:i:109XS:i:30Yt:Z:UU    A01878:56:HKKK3DMXY:2:2411:18909:1642381chr1300583660110M=3005837-109TAAGCCACCACAATCCAATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGTAACTATGATTGGACFFFFFFFFFFF:FFFF,FF,FFFFFF,FFFFF:FFFFFFFFFFF:,FFF,,FFFF:FFFFFF:FFFFFFFF:FF::FFFFFFF:FFFFFFFF:FFFFF:FF:FFFFFFFFNM:i:1MD:Z:0C109MC:Z:110MAS:i:109XS:i:30Yt:Z:UU
A01878:56:HKKK3DMXY:2:2411:18096:23469  chr1    3005837 chr1    3005945 +   -   UU  A01878:56:HKKK3DMXY:2:2411:18096:23469161chr1300583760110M=3005836109AAGCCACCACAATCCAATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGTAACTATGATTGGACAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFF:F:FFFFFFFFFFFFFFFFFFFFNM:i:1MD:Z:109T0MC:Z:110MAS:i:109XS:i:30Yt:Z:UU    A01878:56:HKKK3DMXY:2:2411:18096:2346981chr1300583660110M=3005837-109TAAGCCACCACAATCCAATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGTAACTATGATTGGACFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFNM:i:1MD:Z:0C109MC:Z:110MAS:i:109XS:i:30Yt:Z:UU
A01878:56:HKKK3DMXY:2:2385:11641:5337   chr1    3005845 chr1    3006293 +   -   UU  A01878:56:HKKK3DMXY:2:2385:11641:533797chr1300584560110M=3006184449ACAATCCAATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGTAACTATGATTGGACTTACAATCCFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFNM:i:0MD:Z:110MC:Z:110MAS:i:110XS:i:30Yt:Z:UU    A01878:56:HKKK3DMXY:2:2385:11641:5337145chr1300618460110M=3005845-449ACACGCCTGCTCAAAATGCAGAGTTGTGAAGCCCAGTTACAACTGATATACCTATAACACAAATTCTACACCTAAACCTTGAGGACTATTGTGGAAGAAGGGGTAGAAAGFFFFFFFFFFFFFF,,FFFFFFF,FFF:FF:FFFFFF,FFFFFFFFFFF,FFFFFFFFFFFFF:FFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFNM:i:0MD:Z:110MC:Z:110MAS:i:110XS:i:23Yt:Z:UU
A01878:56:HKKK3DMXY:2:2385:11568:5494   chr1    3005845 chr1    3006293 +   -   UU  A01878:56:HKKK3DMXY:2:2385:11568:549497chr1300584560110M=3006184449ACAATCCAATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGTAACTATGATTGGACTTACAATCCFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFNM:i:0MD:Z:110MC:Z:110MAS:i:110XS:i:30Yt:Z:UU    A01878:56:HKKK3DMXY:2:2385:11568:5494145chr1300618460110M=3005845-449ACACGCCTGCTCAAAATGCAGAGTTGTGAAGCCCAGTTACAACTGATATACCTATAACACAAATTCTACACCTAAACCTTGAGGACTATTGTGGAAGAAGGGGTAGAAAGFF:FFFFFFF:FFFFFFFFFFFFFFF,FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFF:::FFFF:F:FFFF,,FF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFNM:i:0MD:Z:110MC:Z:110MAS:i:110XS:i:23Yt:Z:UU
A01878:56:HKKK3DMXY:2:2325:24930:31046  chr1    3005848 chr1    3005956 +   -   UU  A01878:56:HKKK3DMXY:2:2325:24930:3104697chr1300584860110M=3005847109ATCCAATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGTAACTATGATTGGACTTACAATCCATAFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FNM:i:1MD:Z:109T0MC:Z:110MAS:i:109XS:i:30Yt:Z:UU A01878:56:HKKK3DMXY:2:2325:24930:31046145chr1300584760110M=3005848-109TATCCAATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGTAACTATGATTGGACTTACAATCCATFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFNM:i:1MD:Z:0A109MC:Z:110MAS:i:109XS:i:30Yt:Z:UU
A01878:56:HKKK3DMXY:2:1136:10999:15029  chr1    3005852 chr1    3005928 +   -   UU  A01878:56:HKKK3DMXY:2:1136:10999:15029161chr130058526079M31S=300585277AATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFNM:i:0MD:Z:79MC:Z:33S77MAS:i:79XS:i:30Yt:Z:UU A01878:56:HKKK3DMXY:2:1136:10999:1502981chr130058526033S77M=3005852-77TGACTGGAGTTCAGACGTGTGCTCTTCCGATCTAATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFNM:i:0MD:Z:77MC:Z:79M31SAS:i:77XS:i:30Yt:Z:UU
A01878:56:HKKK3DMXY:2:2136:10881:20431  chr1    3005852 chr1    3005928 +   -   UU  A01878:56:HKKK3DMXY:2:2136:10881:20431161chr130058526079M31S=300585677AATCCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCAGATCGGAAGAGCGTCGTGGAGGGAAAGAGTGTFFFFFFFFF:F:FFF::FFFF:FFFFFFFF:FFF,FF:F:FFF:FFFF:FFFFFF,:FFFF,,,:F,:F,F:FF:FF:FF,FFF:,FFF,:FFFF,,FFFF,,:F:::F,NM:i:0MD:Z:79MC:Z:37S73MAS:i:79XS:i:30Yt:Z:UU A01878:56:HKKK3DMXY:2:2136:10881:2043181chr130058566037S73M=3005852-77TGACTGGAGTTCAGACGTGTGCTCTTCCGATCTAATGCAGCTGTGGCACCTGTGGACCACAACAAAGACTAGCTTGACACAATAACCCTAAGGGTATACTAGTGACATGCFF:,FFF,F,F,,FFF:F:FF,,,:F,F,FF,FF,F,,F,FFF:FF,F,FFF:,FF,,,:,F,,FFF:,,FFF:FFFF,,FFFF,,,:F:,F:FF:::,FF:FF:,,FFFNM:i:0MD:Z:73MC:Z:79M31SAS:i:73XS:i:30Yt:Z:UU

Also when I take only a sample of my sample.03_selected.temp.sam pairs, it works fine. So I don't exactly get where the problem comes from and why it taeks so long to process the second type of pairs/reads.

Best, David

Phlya commented 1 year ago

--max-mismatch 500 - this is a very high value for this argument, do you have a specific biological reason to use 500 there? More typical would be not more than 3... I suspect this is somehow causing this problem...

If you have to use it, you can try using --backend cython.

If you could share the file with the pairs that causes this error, it would help us a lot.

ddepierre commented 1 year ago

Ok it worked fine with --max-mismatch 3 and --backend cython

Can --max-mismatch 500 be used to filter out pairs with a distance <500bp (even if of course it is not the purpose of the dedup function) ? Or do I misunderstand this option?

Phlya commented 1 year ago

Generally, dedup considers pairs with exactly matching coordinates to be duplicates of each other. However there might occasionally be pairs shifted by a few base pairs that might also be duplicates. This option sets the threshold for this inexact matching (e.g. pairs with coordinates (3, 5), (3, 5) are also duplicates, but (3, 5), (3, 6) - only with --max-mismatch 1 or higher.

So in principle you could use a large number if you are using some unusual protocol that would generate such duplicates... But what can happen is that the algorithm tries to annotate all pairs as duplicates of each other and strange things can happen (e.g. it takes forever, it runs out of memory, etc). Cython backend wouldn't have a problem like that, but still you'd end up with no pairs in the end.

Glad this solved your problem!

agalitsyna commented 1 year ago

With --max-mismatch 3 I would recommend switching to --backend cython. Cython and the default engine produce the same results, but Cython is more appropriate for massive numbers of duplicates (usually happens in low-complexity libraries or single-cell analysis).