Could not extend the phase block length.

[Modernism-01]

To tool managers,

Thank you very much for creating this helpful tool in phasing. But I can not get the effective results. this is my command:

$pomfret methphase --vcf $vcf $bam --write-bam -t 8 -o $out/AQ01.chr1.pomfret.methphase

the VCF file has been pahsed and BAM file includes the MM/ML and HP/PS tag.

Here I used one chromosome to test the tool, but I found there are no any changes comparing to original phased blocks.

the below is the part of log file:

[M::main] pomfret v0.1-r14
[M::main] CMD: /public/home/zenglingsen/04.software/48.Pomfret/down/pomfret/pomfret methphase --vcf /public/home/zenglingsen/01.data/01.ONT_data/04.Bam_ONT_for_SNP/05.vcf_p
[M::sancheck_cliopt] read coverage not provided, will estimate.
[M::load_intervals_from_file] at ref NC_010443.5
[T::load_intervals_from_file] used 0.7s
[M::blockjoin_parallel] input has 1 references
[M::blockjoin_parallel] loaded phase block gaps.

[M::estimate_read_coverage_dirtyfast] estimate read depths...
[M::estimate_read_coverage_dirtyfast] NC_010443.5 est. coverage is 20
[M::estimate_read_coverage_dirtyfast] NC_010444.4 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] NC_010445.4 est. coverage is 0

[dbg::haplotag_region_given_bam] loaded 0 reads (interval: NC_010443.5:8647889-8671464), left has #ref=42 (strict: 18), right has=46 (strict: 13)
[dbg::haplotag_region_given_bam] mmrconfig: k 3, k_span 5000, lo 100, hi 156; cov_known 0, cov_for_selection 3, cov_for_runtime 6; readlen thre 15000, min_mapq 10
[dbg::haplotag_region_given_bam] parsed meth sites NC_010443.5:8647889-8671464 (k=3 span=5000 cov=3 runtime_cov=6), fwd n=0, bwd n=0
[W::haplotag_region_given_bam] NC_010443.5:8647889-8671464 does not have methmer in both directions. Skipping.
[dbg::haplotag_region_given_bam] loaded 138 reads (interval: NC_010443.5:8776297-8886673), left has #ref=50 (strict: 17), right has=45 (strict: 17)
[dbg::haplotag_region_given_bam] mmrconfig: k 3, k_span 5000, lo 100, hi 156; cov_known 0, cov_for_selection 3, cov_for_runtime 6; readlen thre 15000, min_mapq 10
[dbg::haplotag_region_given_bam] parsed meth sites NC_010443.5:8776297-8886673 (k=3 span=5000 cov=3 runtime_cov=6), fwd n=997, bwd n=997
[M::haplotag_region_given_bam] left ref ratio: 1.41 (valid=0); right ref ratio: 9.00 (valid=0)
[dbg::haplotag_region1] has 45 ref reads
[dbg::evaluate_separation1] n is 17; 0 0 0 0 (hard coverage fail=0)
[dbg::haplotag_region2] (noperm) xxx undecided (un=1, dir1=0, dir2=0)
[dbg::haplotag_region1] has 50 ref reads
[dbg::evaluate_separation1] n is 17; 1 0 0 0 (hard coverage fail=0)
[dbg::haplotag_region2] (noperm) xxx undecided (un=1, dir1=0, dir2=0)
[dbg::haplotag_region_given_bam] fwd join -1, bwd join -1

[dbg::haplotag_region_given_bam] loaded 0 reads (interval: NC_010443.5:269966250-270082097), left has #ref=44 (strict: 9), right has=48 (strict: 16)
[dbg::haplotag_region_given_bam] mmrconfig: k 3, k_span 5000, lo 100, hi 156; cov_known 0, cov_for_selection 3, cov_for_runtime 6; readlen thre 15000, min_mapq 10
[dbg::haplotag_region_given_bam] parsed meth sites NC_010443.5:269966250-270082097 (k=3 span=5000 cov=3 runtime_cov=6), fwd n=0, bwd n=0
[W::haplotag_region_given_bam] NC_010443.5:269966250-270082097 does not have methmer in both directions. Skipping.
[dbg::haplotag_region_given_bam] loaded 0 reads (interval: NC_010443.5:272918145-273083346), left has #ref=42 (strict: 16), right has=42 (strict: 15)
[dbg::haplotag_region_given_bam] mmrconfig: k 3, k_span 5000, lo 100, hi 156; cov_known 0, cov_for_selection 3, cov_for_runtime 6; readlen thre 15000, min_mapq 10
[dbg::haplotag_region_given_bam] parsed meth sites NC_010443.5:272918145-273083346 (k=3 span=5000 cov=3 runtime_cov=6), fwd n=0, bwd n=0
[W::haplotag_region_given_bam] NC_010443.5:272918145-273083346 does not have methmer in both directions. Skipping.

[M::blockjoin_parallel] done, used 89.8s.
[M::output_gtf] wrote gtf (1 refs, total 536 blocks)
[M::main_blockjoin] gtf written.
[M::main_blockjoin] writing vcf...
[T::load_intervals_from_file] used 0.5s
[M::output_modify_vcf] wrote vcf output, (350052 ok + 0 dropped)/357888 lines modified
[M::main_blockjoin] vcf written.
[M::main_blockjoin] bam written. now indexing...
[M::main_blockjoin] bam index written.

So my results did not show any improvement in length of phase block. The pomfret produced the same result like the whatshap phase based on SNPs. I do not know how to do with this. could you help me to solve it?

Look forward to you reply.

Best,

Lingsen

[xfengnefx]

Hi Lingsen, it seems the run finished normally but pomfret was not confident enough to join any of the gaps. The heuristic is conservative as to avoid switch errors and in this dataset it decided to not close the gaps.

[Modernism-01]

Thanks for your reply. That means my data quality does not have high confidence to close the gaps.

But here I tested your released test data. But I did not get the same results in your test data. the below is the command and log file:

pomfret methphase --vcf variants.vcf -t 8 -o test.zls phased.bam

[M::main] pomfret v0.1-r14
[M::main] CMD: /public/home/zenglingsen/04.software/48.Pomfret/down/pomfret/pomfret methphase --vcf variants.vcf phased.bam -t 8 -o test.zls

[M::sancheck_cliopt] read coverage not provided, will estimate.
[M::load_intervals_from_file] at ref chr6
[T::load_intervals_from_file] used 0.0s
[M::blockjoin_parallel] input has 1 references
[M::blockjoin_parallel] loaded phase block gaps.

[M::estimate_read_coverage_dirtyfast] estimate read depths...
[M::estimate_read_coverage_dirtyfast] chr1 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr2 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr3 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr4 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr5 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr6 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr7 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr8 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr9 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr10 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr11 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr12 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr13 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr14 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr15 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr16 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr17 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr18 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr19 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr20 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr21 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chr22 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chrX est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chrY est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chrM est. coverage is 0
......

[M::estimate_read_coverage_dirtyfast] chrUn_GL000218v1 est. coverage is 0
[M::estimate_read_coverage_dirtyfast] chrEBV est. coverage is 0
[T::estimate_read_coverage_dirtyfast] used 0.0s
[W::blockjoin_one_chrom_callback] had to clamp n_candidates_per_iter (ref: chr6)
[dbg::blockjoin_one_chrom_callback] ref chr6 using: cov_for_selection=1, n_cand_per_iter=2
[dbg::haplotag_region_given_bam] loaded 199 reads (interval: chr6:11092382-11147866), left has #ref=74 (strict: 45), right has=62 (strict: 44)
[dbg::haplotag_region_given_bam] mmrconfig: k 3, k_span 5000, lo 100, hi 156; cov_known 0, cov_for_selection 1, cov_for_runtime 2; readlen thre 15000, min_mapq 10
[dbg::haplotag_region_given_bam] parsed meth sites chr6:11092382-11147866 (k=3 span=5000 cov=1 runtime_cov=2), fwd n=1516, bwd n=1516
[M::haplotag_region_given_bam] left ref ratio: 1.00 (valid=1); right ref ratio: 1.00 (valid=1)
[dbg::haplotag_region1] has 62 ref reads
[dbg::evaluate_separation1] n is 45; 0 21 0 24 (hard coverage fail=0)
[dbg::haplotag_region2] (noperm) xxx undecided (un=1, dir1=0, dir2=0)
[dbg::haplotag_region1] has 74 ref reads
[dbg::evaluate_separation1] n is 44; 24 0 20 0 (hard coverage fail=0)
[dbg::haplotag_region2] (noperm) xxx undecided (un=1, dir1=0, dir2=0)
[dbg::haplotag_region_given_bam] fwd join -1, bwd join -1

[M::blockjoin_parallel] done, used 0.2s.
[M::output_gtf] wrote gtf (1 refs, total 2 blocks)
[M::main_blockjoin] gtf written.
[M::main_blockjoin] writing vcf...
[T::load_intervals_from_file] used 0.0s
[M::output_modify_vcf] wrote vcf output, (24 ok + 0 dropped)/346 lines modified
[M::main_blockjoin] vcf written.

[M::main] CMD: /public/home/zenglingsen/04.software/48.Pomfret/down/pomfret/pomfret methphase --vcf variants.vcf -t 8 -o test.zls phased.bam
[M::main] used: 0.3s, peak RSS 0.0GiB

The below is the three results for your original phased block, your pomfret results, and my pomfret results:

your original phased block

Phasing statistics for sample hg002 from file variants.vcf
---------------- Chromosome chr6 ----------------
Variants in VCF: 79
Heterozygous: 31 ( 25 SNVs)
Phased: 25 ( 25 SNVs)
Heterozygous phased: 80.6% ( 100.0% SNVs)
Unphased: 6 (not considered below)
Singletons: 0 (not considered below)
Blocks: 2

Block sizes (no. of variants)
Sum of sizes: 25 variants
Median block size: 12.50 variants
Average block size: 12.50 variants
Largest block: 22 variants
Smallest block: 3 variants

Block lengths (basepairs)
Sum of lengths: 16206 bp
Median block length: 8103.00 bp
Average block length: 8103.00 bp
Longest block: 9691 bp
Shortest block: 6515 bp
Block NG50: 0 bp

your pomfret results

Phasing statistics for sample hg002 from file output.mp.vcf
---------------- Chromosome chr6 ----------------
Variants in VCF: 79
Heterozygous: 31 ( 25 SNVs)
Phased: 25 ( 25 SNVs)
Heterozygous phased: 80.6% ( 100.0% SNVs)
Unphased: 6 (not considered below)
Singletons: 0 (not considered below)
Blocks: 1

Block sizes (no. of variants)
Sum of sizes: 25 variants
Median block size: 25.00 variants
Average block size: 25.00 variants
Largest block: 25 variants
Smallest block: 25 variants

Block lengths (basepairs)
Sum of lengths: 71690 bp
Median block length: 71690.00 bp
Average block length: 71690.00 bp
Longest block: 71690 bp
Shortest block: 71690 bp
Block NG50: 0 bp

my pomfret results

Phasing statistics for sample hg002 from file test.zls.mp.vcf
---------------- Chromosome chr6 ----------------
Variants in VCF: 79
Heterozygous: 31 ( 25 SNVs)
Phased: 24 ( 24 SNVs)
Heterozygous phased: 77.4% ( 96.0% SNVs)
Unphased: 6 (not considered below)
Singletons: 1 (not considered below)
Blocks: 2

Block sizes (no. of variants)
Sum of sizes: 24 variants
Median block size: 12.00 variants
Average block size: 12.00 variants
Largest block: 21 variants
Smallest block: 3 variants

Block lengths (basepairs)
Sum of lengths: 69139 bp
Median block length: 34569.50 bp
Average block length: 34569.50 bp
Longest block: 61018 bp
Shortest block: 8121 bp
Block NG50: 0 bp

Would you mind have a look?

Thank you very much for your time and effort.

Best,

Lingsen

[xfengnefx]

For the example dataset you would need a -c 60 as mentioned in the readme. Sorry about the confusion. I will consider to update the example dataset for a smoother start.

I've also pushed fixups for a corner case and others, would be great if you could try out the current main branch. Thank you!

[Modernism-01]

Fine. it works. However, in my pigs data, I can not extend the length of phase block using methylation information. Is it relative to sequencing depth (~20X in my data) or quality contril? Do you have any comments or suggestions to improve it?

[xfengnefx]

Sorry for the late reply, I was away. Thanks for your patience. I think for your case, the sequencing depth might be too low. But just in case, if you still have the log (i.e. the STDERR output) from pomfret, could you paste it here so I can check if anything unexpected has happend for the pomfret run? Thank you!

[Modernism-01]

Thanks for your reply. Did you mean this log file below? I guess there are only one log file output.

19085706.err

[xfengnefx]

@Modernism-01 I've requested for viewer permission for this file (from xfengnefx at gmail), please give it a go when you are available. TIA!

[xfengnefx]

Thanks for your patience. From the log, it seems that the read depth is on the lower side, and for many phasing breakpoints, pomfret could not collect enough methylation info from reads to attempt to close the gaps. Relaxing the heuristics might let pomfret to make more attempts, but I believe that will risk making excessive switch errors.

I would suggest to either leave the phasing as-is, or try other tools or add other data types to improve phasing.