Inquiry Regarding PECAT Phasing and Assembly Interpretation

[HLHsieh]

Hi,

Thank you for developing such a helpful tool. I’m planning to use it for local assembly of a repetitive region and would appreciate your insights on a couple of questions related to its usage and feasibility for my case.

1. Does PECAT provide the phase result of each read in dual-mode assembly? Specifically, I would like to know whether it assigns each read to a specific haplotype during the assembly process.

2. For my test, I simulated reads with lengths shorter than the known repetitive region I intend to assemble. Given this input, can PECAT reliably assemble the entire repetitive region? Below are the settings I used for assembly:

project=P8_dual
reads=./P8_extracted_region.fasta
genome_size=300000
threads=8
cleanup=1
grid=local
prep_min_length=3000
prep_output_coverage=50
corr_iterate_number=1
corr_block_size=4000000000
corr_correct_options=
corr_filter_options=--filter0=l=5000:al=2500:alr=0.5:aal=5000:oh=3000:ohr=0.3
corr_rd2rd_options=-x ava-ont
corr_output_coverage=80
align_block_size=4000000000
align_rd2rd_options=-X -g3000 -w30 -k19 -m100 -r500
align_filter_options=--filter0=l=5000:aal=6000:aalr=0.5:oh=3000:ohr=0.3 --task=extend --filter1=oh=300:ohr=0.03
asm1_assemble_options=
phase_method=
phase_rd2ctg_options=-x map-ont -c -p 0.5 -r 1000
phase_use_reads=1
phase_phase_options=
phase_clair3_command=singularity exec -B `pwd -P`:`pwd -P` -B /tmp:/tmp clair3_v0.1-r12.sif /opt/bin/run_clair3.sh
phase_clair3_use_reads=0
phase_clair3_options=--platform=ont --model_path=/opt/models/ont_guppy5/  --include_all_ctgs
phase_clair3_rd2ctg_options=-x map-ont -w10 -k19 -c -p 0.5 -r 1000
phase_clair3_phase_options=--coverage lc=30 --phase_options icr=0.1:icc=6:sc=10 --filter i=70
phase_clair3_filter_options=--threshold=2500 --rate 0.05
asm2_assemble_options=--contig_format dual,prialt
polish_map_options=-x map-ont
polish_filter_options=--filter0 oh=1000:ohr=0.1
polish_cns_options=
polish_medaka=
polish_medaka_command = singularity exec -B `pwd -P`:`pwd -P`  -B /tmp:/tmp medaka_v1.7.2.sif medaka
polish_medaka_map_options=-x map-ont
polish_medaka_cns_options=--model r941_prom_sup_g507

3. I noticed that PECAT can be run without using Clair3 or Medaka. Could you clarify the differences in phasing performance between the built-in method and using Clair3? Additionally, how much improvement can be expected by enabling Medaka polishing (polish_medaka=1)?

4. I would appreciate some clarification on how to interpret the file located in 5-assemble/haplotype_1_tiles. Specifically, for my use case, I am performing a local assembly on a specific region (~50 kb), and I noticed that the assembly generates several contigs. Could you advise on how to interpret this output, particularly in the context of the generated contigs?

ctg7 edge=9364:B~5840:E
ctg7 edge=5840:E~978:B
ctg3 edge=2237:E~1756:E
ctg3 edge=1756:E~3895:E
ctg3 edge=3895:E~8776:B
ctg3 edge=8776:B~10451:B
ctg1 edge=11619:E~2988:B
ctg1 edge=2988:B~1562:E
ctg1 edge=1562:E~11796:B
ctg12 edge=8148:E~15022:B
ctg12 edge=15022:B~10968:B
ctg12 edge=10968:B~10451:B
ctg0 edge=12123:B~10881:E
ctg0 edge=10881:E~12073:B

Thank you in advance for your time and assistance. I look forward to your feedback.

Best regards,
Hsin

[lemene]

@HLHsieh

1. PECAT doesn't assign each read to a specific haplotype. Instead, it identifies the inconsistent overlaps which connect reads from different haplotype.
2. In theory, PECAT does not assume a specific copy number for the repetitive region; however, the performance of PECAT in solving the repetitive region with multiple copies remains unknown.
3. The use of Clair3 has a impact on the assembly of low heterozygosity datasets. PECAT may overcorrects reads in low heterozygosity regions (<0.0005), so it uses clair3 to identify SNPs from raw reads for better results. For nanopore data, medaka is recommended to improve base-level quality.
4. *_tiles records how the contigs are constructed. Each contig is a path in the string graph. edge=9364:B~5840:E is an edge in the string graph. 9364 and 5840 are ids of the reaads. B and E indicate the two ends of the reads. More details about the string graph can be found at https://academic.oup.com/bioinformatics/article/21/suppl_2/ii79/227189

[HLHsieh]

Hi,

Thank you for your prompt response. I have several follow-up questions after testing with some datasets:

1. Regarding question 1, is it possible to assign each read to a specific haplotype based on inconsistent information on my own? What is your opinion on this? Additionally, how should I interpret the results under the /4-phase/fsa folder?

$inconsistent
4837:9197|1|39325,905|0|2320,16449|0|2331,7628|1|39364,4542|1|39398,6173|1|39314,
16449:10898|0|-2660,16082|0|-2621,4837|0|-2331,11721|1|39180,5323|0|-2482,4948|0|-2367,19406|0|-2659,16704|0|-2650,18764|1|39144,
16082:9197|1|39619,905|0|2610,16449|0|2621,7628|1|39655,4542|1|39712,6173|1|39608,
905:10898|0|-2841,16082|0|-2610,4837|0|-2320,11721|1|39000,5323|0|-2471,4948|0|-2423,19406|0|-2803,16704|0|-2831,18764|1|38966,
9197:10898|1|39843,16082|1|39619,4837|1|39325,11721|0|-1999,5323|1|39479,4948|1|39422,19406|1|39804,16704|1|39832,18764|0|-1964,
18764:16449|1|39144,4542|0|1866,905|1|38966,7628|0|1937,9197|0|1964,6173|0|1961,
6173:16082|1|39608,4837|1|39314,11721|0|-1996,5323|1|39468,10898|1|39845,4948|1|39425,19406|1|39806,16704|1|39835,18764|0|-1961,
4948:905|0|2423,16449|0|2367,7628|1|39468,4542|1|39522,9197|1|39422,6173|1|39425,
5323:9197|1|39479,905|0|2471,16449|0|2482,7628|1|39517,4542|1|39553,6173|1|39468,
19406:9197|1|39804,905|0|2803,16449|0|2659,7628|1|39848,4542|1|39902,6173|1|39806,

$readinfos
ctg5 16704 0 -12888  3961|-1|16849|3 5821|-1|18706|2 6804|-1|19690|2 8309|-1|21196|2 8562|-1|21450|2 10611|-1|23409|2 11545|-1|24275|2 13595|-1|26319|0 15773|2|28438|2 16290|1|28956|1 16380|3|29046|3 16390|0|29056|0 16404|1|29070|1 16433|1|29099|1 16447|0|29113|0 16465|3|29131|3 16487|3|29153|3 16498|1|29162|1 16589|2|29249|2 16591|0|29251|0 16633|3|29291|3 16636|1|29294|1 16638|1|29296|1 16641|1|29299|1 16682|3|29339|3 16698|1|29355|1 16732|3|29387|3 16743|1|29397|1 16745|0|29399|0 16748|3|29402|3 16751|0|29405|0 16753|1|29407|1 16797|2|29446|2 16809|0|29458|0 16812|0|29461|0 16817|1|29466|1 16837|0|29485|0 16846|0|29494|0 16863|-1|29510|3 16875|3|32364|3 16880|0|32370|0 16888|0|32378|0 16891|1|32381|1 16895|2|32385|2 16907|3|32397|3 16954|1|32444|1 16958|2|32448|2 16976|3|32465|3 16979|1|32468|1 16982|1|32471|1 16984|3|32473|3 16992|0|32482|0 16995|2|32485|2 17004|2|32494|2 17034|3|32524|3 17051|2|32540|2 17078|1|32564|1 17093|1|32584|1 17096|1|32587|1 17108|2|32598|2 17120|2|32609|2 17129|0|32618|0 17136|1|32625|1 17140|3|32629|3 17142|3|32631|3 17150|1|32639|1 17162|3|32651|3 17165|3|32654|3 17175|1|32666|1 17185|3|32676|3 17189|3|32680|3 17195|1|32686|1 17204|1|32695|1 17221|1|32712|1 17223|0|32714|0 17235|3|32726|3 17252|1|32743|1 17264|1|32755|1 17267|2|32758|2 17269|3|32760|3 17273|1|32764|1 17285|2|32776|2 17313|1|32803|1 17326|3|32816|3 17342|1|32832|1 17351|2|32841|2 17356|1|32845|1 17366|0|32855|0 17371|0|32860|0 17387|2|32876|2 17391|1|32880|1 17395|3|32884|3 17415|2|32904|2 17417|1|32906|1 17428|3|32917|3 17436|1|32925|1 17443|3|32932|3 17460|1|32948|1 17478|1|32966|1 17503|3|32991|3 17527|3|33015|3 17529|3|33017|3 17533|1|33021|1 17544|3|33032|3 17589|3|33077|3 17598|0|33086|0 17600|1|33088|1 17617|3|33105|3 17645|3|33133|3 17658|2|33146|2 17676|1|33164|1 17681|3|33169|3 17722|1|33210|1 17730|2|33218|2 17756|2|33244|2 17785|0|33273|0 17799|2|33287|2 17821|1|33378|1 17829|0|33386|0 17853|3|33410|3 17873|2|33430|2 17876|0|33433|0 17885|2|33442|2 17895|3|33452|3 17903|1|33460|1 17914|2|33471|2 17952|2|33508|2 17955|1|33511|1 17986|3|33542|3 17993|3|33549|3 18015|2|33571|2 18072|2|33628|2 18077|0|33633|0 18081|0|33637|0 18083|1|33639|1 18095|1|33651|1 18110|1|33666|1 18114|1|33670|1 18119|1|33675|1 18125|3|33681|3 18132|0|33688|0 18139|1|33695|1 18176|0|33732|0 18183|2|33739|2 18192|0|33748|0 18201|2|33757|2 18205|1|33761|1 18240|0|33796|0 18251|1|33807|1 18254|1|33810|1 18282|1|33841|1 18290|0|33849|0 18307|2|33866|2 18309|1|33868|1 18352|3|33911|3 19285|1|36956|1 19312|3|36983|3 19344|1|37015|1 19440|0|37111|0 19683|1|37354|1 19751|2|37422|2 19813|1|37484|1 20280|-1|37951|0 20345|3|38016|3 20451|-1|38122|2 20469|2|38140|2 20723|2|38394|2 20766|0|38437|0 21448|2|39122|2 21624|2|39299|2 21800|3|39475|3 21890|0|39565|0 21901|3|39576|3 22261|3|39936|3 22300|0|39975|0 22474|2|40149|2 22523|1|40198|1 22597|2|40272|2 22813|2|40489|2 22977|3|40652|3 22992|3|40667|3 23282|0|40957|0 23404|2|41079|2 23725|3|41413|3 23795|1|41483|1 23816|-1|41504|3 23838|2|41526|2 23872|2|41560|2 23924|1|41612|1 24060|3|41748|3 24088|3|41776|3
ctg5 4542 1 24054  3961|-1|20093|3 5821|-1|18234|2 6804|-1|17251|2 8309|-1|15749|2 8562|-1|15495|2 10611|-1|13470|1 11545|-1|12549|2 13595|-1|10502|2 15773|-1|8379|2 16290|-1|-1|-1 16380|-1|7748|2 16390|-1|7738|0 16404|-1|7724|1 16433|-1|7695|1 16447|-1|7681|0 16465|-1|7663|3 16487|-1|7641|3 16498|-1|7632|1 16589|-1|7546|2 16591|-1|7544|0 16633|-1|7504|3 16636|-1|7501|1 16638|-1|7499|1 16641|-1|7496|1 16682|-1|7456|3 16698|-1|-1|-1 16732|-1|7407|3 16743|-1|7397|1 16745|-1|7395|0 16748|-1|7392|3 16751|-1|7389|0 16753|-1|7387|1 16797|-1|7343|2 16809|-1|-1|-1 16812|-1|7329|0 16817|-1|-1|-1 16837|-1|7306|2 16846|-1|7297|1 16863|-1|7281|3 16875|-1|7267|1 16880|-1|7262|3 16888|-1|7254|3 16891|-1|7251|0 16895|-1|7247|1 16907|-1|7237|0 16954|-1|7190|0 16958|-1|7186|3 16976|-1|7168|1 16979|-1|7165|3 16982|-1|7162|2 16984|-1|7160|2 16992|-1|7152|1 16995|-1|7149|0 17004|-1|7140|0 17034|-1|7110|0 17051|-1|7093|0 17078|-1|7066|3 17093|-1|7051|3 17096|-1|7048|3 17108|-1|7036|0 17120|-1|7024|3 17129|-1|7015|3 17136|-1|7008|3 17140|-1|7004|2 17142|-1|7002|2 17150|-1|6994|0 17162|-1|6982|0 17165|-1|6979|2 17175|-1|6969|0 17185|-1|6959|0 17189|-1|6955|2 17195|-1|6949|2 17204|-1|6940|3 17221|-1|6923|0 17223|-1|6921|1 17235|-1|6909|1 17252|-1|6892|0 17264|-1|6880|3 17267|-1|6877|0 17269|-1|6875|1 17273|-1|6871|3 17285|-1|6859|1 17313|-1|6831|3 17326|-1|6818|1 17342|-1|6802|3 17351|-1|6793|3 17356|-1|6788|2 17366|-1|6779|1 17371|-1|6774|1 17387|-1|6758|1 17391|-1|6754|3 17395|-1|6750|1 17415|-1|6730|3 17417|-1|6728|3 17428|-1|6717|1 17436|-1|6709|3 17443|-1|6702|1 17460|-1|6685|3 17478|-1|6667|2 17503|-1|6642|1 17527|-1|6618|1 17529|-1|6616|0 17533|-1|6612|3 17544|-1|6601|1 17589|-1|6556|2 17598|-1|6547|3 17600|-1|6545|0 17617|-1|6528|1 17645|-1|6500|2 17658|-1|6487|3 17676|-1|6469|3 17681|-1|6464|1 17722|-1|6423|2 17730|-1|6415|3 17756|-1|6389|3 17785|-1|6360|3 17799|-1|6346|1 17821|-1|6324|3 17829|-1|6316|3 17853|-1|6292|2 17873|-1|6272|0 17876|-1|6269|3 17885|-1|6260|0 17895|-1|6250|0 17903|-1|6242|3 17914|-1|6231|3 17952|-1|6193|0 17955|-1|6190|3 17986|-1|6159|1 17993|-1|6152|0 18015|-1|6130|3 18072|-1|6073|1 18077|-1|6068|2 18081|-1|6064|1 18083|-1|6062|2 18095|-1|6050|0 18110|-1|6035|3 18114|-1|6031|2 18119|-1|6026|2 18125|-1|6020|2 18132|-1|6013|2 18139|-1|6006|3 18176|-1|5969|1 18183|-1|5962|1 18192|-1|5953|3 18201|-1|5944|0 18205|-1|5940|3 18240|-1|5905|2 18251|-1|5894|2 18254|-1|5891|3 18282|-1|5862|3 18290|-1|5854|1 18307|-1|5837|0 18309|-1|5835|0 18352|-1|5792|1 19285|-1|4864|3 19312|-1|4837|0 19344|-1|4805|3 19440|-1|4709|1 19683|-1|4466|2 19751|-1|4398|0 19813|-1|4336|0 20280|-1|3869|2 20345|-1|3804|2 20451|-1|3698|2 20469|-1|3680|0 20723|-1|3426|0 20766|-1|3383|1 21448|-1|2699|1 21624|-1|2523|3 21800|-1|2347|2 21890|-1|2257|3 21901|-1|2246|0 22261|-1|1886|1 22300|-1|1847|3 22474|-1|1673|3 22523|-1|1624|2 22597|-1|1550|3 22813|-1|1334|1 22977|-1|1171|2 22992|-1|1156|2 23282|-1|866|2 23404|-1|744|0 23725|-1|421|1 23795|-1|351|3 23816|-1|330|1 23838|-1|308|3 23872|-1|274|3 23924|-1|222|3 24060|-1|86|1 24088|-1|58|1

2. Regarding question 3, I am wondering how PECAT differentiates haplotype_1 and haplotype_2 in dual mode without considering SNPs. From my understanding, PECAT considers inconsistent overlaps as a core component. However, I am having a hard time understanding the essential difference between SNPs and inconsistent overlaps. Are inconsistent overlaps derived from SNPs? Besides, if phase_method=2 is enabled, are Clair3 and the built-in method executed sequentially or in parallel? Specifically, in my case, I encountered errors when using phase_method=2 (details of the errors are provided below). However, it works correctly when phase_method=0 or phase_method=1 is selected.

2025-01-30 23:04:08 [Info] Run script /scratch/hsinlun/P2_filter.1/scripts/pol_filter_hap.sh
2025-01-30 23:04:08 [Info] Sumbit command: /scratch/hsinlun/P2_filter.1/scripts/pol_filter_hap.sh 2>&1 | tee /scratch/hsinlun/P2_filter.1/scripts/pol_filter_hap.sh.log
2025-01-30 23:04:08 [Error] Failed to run filter overlaps in which the reads are different haplotype

Looking forward to your clarification and suggestions.

Best regards,
Hsin

[lemene]

@HLHsieh
1.
inconsistent files record the the inconsistent overlaps. For example 4837:9197|1|39325 means that the overlap between read 4837 and read 9197 is inconsistent. 1 means a reverse complement alignment. 39325 is the position where the start of read 9197 is mapped to read 4837 based on the overlap.

readinfos record the SNP alleles in the reads.
ctg5 16704 0 -12888 3961|-1|16849|3
ctg5 and 16704 are the ids of the contig and the read.
0 or 1 indicates whether it is a reverse complement.
-12888 is the position where the start of read 16704 is mapped to ctg5.
3961|-1|16849|3 are SNP alleles and their positions in the read and the contig

Inconsistent overlaps are derived from SNPs. If two reads contain different SNP alleles, they are considered to have an inconsistent overlap, assuming there are no sequencing errors. Clair3 and the built-in method executed sequentially. pol_filter_hap.sh failed to run, and its log is recorded in scripts/pol_filter_hap.sh.log that can help us locate the error.