search

vibansal / HapCUT2

software tools for haplotype assembly from sequence data
BSD 2-Clause "Simplified" License
207 stars 36 forks source link

HapCUT2 does not phase this variant directly adjacent to another and overlapping another other. #132

Open jan-glx opened 1 year ago

jan-glx commented 1 year ago

For below data example with two haplotypes, each with a different overlapping indel variant, that are represented in the unphased input vcf file as three variants (a one 1 bp deletion immediately followed by a 8 bp deletion for haplotype "earlyIndel" and a 5 bp deletion for haplotype "lateDeletion") HapCUT2 (with --indel 1) leaves the the second variant unphased. Is this expected?

image

Fully reproducible example (using current master v1.3.1 0eb7075424afbb40259ced6b5b1cbd95d0cedf55 and current develop of htslib: 1.18 99415e2a)

wget https://hgdownload.cse.ucsc.edu/goldenpath/mm10/chromosomes/chr19.fa.gz && gunzip chr19.fa.gz
samtools faidx chr19.fa
bind 'set disable-completion on' # temporarily disable auto-completion to allow pasting tabs.
cat <<'EOF' > unphased.vcf
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele not already represented at this location by REF and ALT">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##contig=<ID=chr19,length=61431566>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  CGAAGAGGTAGGTGCGAG-1
chr19   55910646    .   AC  A   7447.85 PASS    AC=1    GT:AD:DP    0/1:47,35:82
chr19   55910648    .   AAATCCCCC   A   12412.8 PASS    AC=1    GT:AD:DP    0/1:47,35:82
chr19   55910653    .   CCCCAT  C   13372.2 PASS    AC=1    GT:AD:DP    0/1:0,47:82
EOF

cat <<'EOF' > reads.sam
@HD VN:1.4  SO:coordinate
@SQ SN:chr19    LN:61431566
@RG ID:CGAAGAGGTAGGTGCGAG-1 SM:CGAAGAGGTAGGTGCGAG-1
@PG ID:bwa  VN:0.7.12-r1039 CL:bwa mem -C -M -t 32 -H .sam.header .fa .fastq.gz .fastq.gz   PN:bwa
A01382:376:H5HVGDRX3:1:2131:30807:22842 147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGAAATCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTCTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFF    MC:Z:109M   MD:Z:65^CAAATCCCC0C54T39    RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:11 MQ:i:60 AS:i:135    XS:i:19
A01382:376:H5HVGDRX3:1:2252:12048:7200  147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGAAATCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF    MC:Z:109M   MD:Z:65^CAAATCCCC0C94   RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:10 MQ:i:60 AS:i:140    XS:i:20
A01382:376:H5HVGDRX3:1:2264:23258:10457 147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGAAATCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFF,FFFFFFFFFFFFFF:FFFFFFF:FFFF:FFF:,FFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFF:FF,FFFFFFFFFFF,FFF:F:FF:FFFFFF:FFFFF:F:,:FF,FFFF:FFFFFFF:FFFFFF:FFFF,FFFFFFFF:FFF,FFFF    MC:Z:109M   MD:Z:65^CAAATCCCC0C94   RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:10 MQ:i:60 AS:i:140    XS:i:20
A01382:376:H5HVGDRX3:2:2105:19253:5024  147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGAAATCCGCTAGGATGGTTAGCACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFF:F:FFF::FFFFFFFFFFFFFFF:FFF:FFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFF    MC:Z:109M   MD:Z:65^CAAATCCCC0C18T75    RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:11 MQ:i:60 AS:i:135    XS:i:19
A01382:376:H5HVGDRX3:2:2105:19271:5431  147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGAAATCCGCTAGGATGGTTAGCACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFF:FFFFFFFF    MC:Z:109M   MD:Z:65^CAAATCCCC0C18T75    RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:11 MQ:i:60 AS:i:135    XS:i:19
A01382:376:H5HVGDRX3:2:2167:16938:3302  147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGACATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTGGGAAATCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FF:FFFFF,F:FFFFFFFFFFF::FF:FFFFFFF,F,,FFFFFFFFFFFFFFFF,FF:FFFF,FFFFF:FFFFFFFFFFF:FFF,FF:FFFFFFFFFF:FFFFFFFF:::FFFFFFFF:FFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFF    MC:Z:109M   MD:Z:25T35A3^CAAATCCCC0C94  RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:12 MQ:i:60 AS:i:130    XS:i:20
A01382:376:H5HVGDRX3:2:2231:9245:11741  147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGAAATCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFF:FFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFF    MC:Z:109M   MD:Z:65^CAAATCCCC0C94   RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:10 MQ:i:60 AS:i:140    XS:i:20
A01382:376:H5HVGDRX3:2:2236:29125:20431 147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGAAATCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FF:FFFFFFFFFFFFFFFFFFFFFFF:FFFFF,F,F:,FFFFFFF:FFFFFFF:F:FFFFFFFFF:F:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF,:,:FFFFFFFFF:FFFFFFFFFFFFFFF:FFFF::FFFFFFFFF    MC:Z:109M   MD:Z:65^CAAATCCCC0C94   RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:10 MQ:i:60 AS:i:140    XS:i:20
A01382:376:H5HVGDRX3:2:2259:3025:9549   147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGAAATCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF,FFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,:FFFFFFFF    MC:Z:109M   MD:Z:65^CAAATCCCC0C94   RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:10 MQ:i:60 AS:i:140    XS:i:20
A01382:376:H5HVGDRX3:2:2274:11460:6934  147 chr19   55910582    60  65M9D95M    =   55910535    -216    TCCCAAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGAAATCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFF    MC:Z:109M   MD:Z:65^CAAATCCCC0C94   RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:10 MQ:i:60 AS:i:140    XS:i:20
A01382:376:H5HVGDRX3:1:2247:22309:12164 147 chr19   55910586    60  68M5D92M    =   55910535    -216    AAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGACAAATCCCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    :FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:,FFFFFFFFFFFFFFFFFFFFFFFFFFF:F:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFF    MC:Z:109M   MD:Z:68^CCCAT92 RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:5  MQ:i:60 AS:i:149    XS:i:20
A01382:376:H5HVGDRX3:2:2137:18249:27336 147 chr19   55910586    60  68M5D92M    =   55910535    -216    AAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGACAAATCCCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    :FFFFFFFFFFFFFFFFFFFFFF,FFFFFFFF:F:,FFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFF,F:FFFFFFFF    MC:Z:109M   MD:Z:68^CCCAT92 RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:5  MQ:i:60 AS:i:149    XS:i:20
A01382:376:H5HVGDRX3:2:2206:5638:6026   147 chr19   55910586    60  68M5D92M    =   55910535    -216    AAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGACAAATCCCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFF:,,F:FF:F,FFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,:FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF    MC:Z:109M   MD:Z:68^CCCAT92 RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:5  MQ:i:60 AS:i:149    XS:i:20
A01382:376:H5HVGDRX3:2:2227:24740:17628 147 chr19   55910586    60  68M5D92M    =   55910535    -216    AAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGACAAATCCCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFFFFFFFFFFFFFFFFFFFFFFF:FFFFF::F,FFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:F,FFFFFFFF    MC:Z:109M   MD:Z:68^CCCAT92 RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:5  MQ:i:60 AS:i:149    XS:i:20
A01382:376:H5HVGDRX3:2:2252:15700:25144 147 chr19   55910586    60  68M5D92M    =   55910535    -216    AAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGACAAATCCCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFF    MC:Z:109M   MD:Z:68^CCCAT92 RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:5  MQ:i:60 AS:i:149    XS:i:20
A01382:376:H5HVGDRX3:2:2255:32018:34898 147 chr19   55910586    60  68M5D92M    =   55910535    -216    AAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGACAAACCCCCGCTAGGATGGTTGGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    ,:FFFFFFFFFFFFFFF:FFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF:FFFFFFFFFFFFFFFFF:F,:FFFF:FF    MC:Z:109M   MD:Z:65T2^CCCAT14A77    RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:7  MQ:i:60 AS:i:139    XS:i:19
A01382:376:H5HVGDRX3:2:2270:19180:5181  147 chr19   55910586    60  68M5D92M    =   55910535    -216    AAGGCCTCCGCACCCTCCAGATATCTCTCCATATTACCCGCTGTCGCCCGGCACCGTAGGACAAATCCCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTCTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    FF,FFFFF:FFFFFFFFFFFF,FF:FFFFFFF:F,FFFFFFFFFFFFFFF:FFFFFFFFFFFF:FF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFF:FFFFFFFFFFFFFFFF:FF,F:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF    MC:Z:109M   MD:Z:68^CCCAT56T35  RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:6  MQ:i:60 AS:i:144    XS:i:19
A01382:376:H5HVGDRX3:1:2145:13051:6527  147 chr19   55910620    60  34S34M5D92M =   55910535    -216    CCCGCCCCCCCACCCACCAGATCTCTCGCCACACTACCCGCCGTCGCCCGGCACCGTAGGACAACTCCCCGCTAGGATGGTTAGTACCACAGTAAGCAATTCCAGTTTTTAATTCCCTTTTTGTTTCTTGCATGAGCATGCTTATTAGTTTACGTGTACG    :,,,FF,FF,F,FFF,FFF,:,,,F,F,FF,,,,,FFFFFF,F:F:FFFFFF:FFF:FF:FF:F,FFFFF,FFFF,:FF:F:FF,FFFF:FF,F:F:FF:FFFFFFFF:F,::FFF::FF:FFF,FFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFF    MC:Z:109M   MD:Z:7T22A3^CCCAT92 RG:Z:CGAAGAGGTAGGTGCGAG-1   NM:i:7  MQ:i:60 AS:i:105    XS:i:20
EOF

bind 'set disable-completion off'

extractHAIRS --bam reads.sam --VCF unphased.vcf --indels 1 --ref chr19.fa --out fragment_file --verbose 1 && HAPCUT2 --fragments fragment_file --VCF unphased.vcf --verbose 1 --out haplotype_output_file && cat haplotype_output_file.phased.VCF
Extracting haplotype informative reads from bamfiles reads.bam minQV 13 minMQ 20 maxIS 1000

VCF file unphased.vcf has 3 variants
adding chrom chr19 to index
vcffile unphased.vcf chromosomes 1 hetvariants 3 variants 3
detected 0 variants with two non-reference alleles, these variants will not be phased
reading fasta index file chr19.fa.fai ... fasta file chr19.fa has 1 chromosomes/contigs

found match for reference contig chr19 in VCF file index
contig chr19 length 61431566
reading reference sequence file chr19.fa with 1 contigs
read reference sequence file in 0.37 sec
reading sorted bamfile reads.bam
processing reads mapped to chrom "chr19"
final cleanup of fragment list: 8 current chrom 0 prev 0

[2023:10:26 11:04:40] input fragment file: fragment_file
[2023:10:26 11:04:40] input variantfile (VCF format):unphased.vcf
[2023:10:26 11:04:40] haplotypes will be output to file: haplotype_output_file
[2023:10:26 11:04:40] solution convergence cutoff: 5
[2023:10:26 11:04:40] read 3 variants from unphased.vcf file
[2023:10:26 11:04:40] read fragment file and variant file: fragments 7 variants 3
mean number of variants per read is 2.00
[2023:10:26 11:04:40] building read-variant graph for phasing
[2023:10:26 11:04:40] no of non-trivial connected components 1 max-Degree 7 connected variants 2 coverage-per-variant 7.000000
[2023:10:26 11:04:40] comp 0 first 0 last 2 phased 2 fragments 7
[2023:10:26 11:04:40] fragments 7 snps 3 component(blocks) 1
[2023:10:26 11:04:40] starting Max-Likelihood-Cut based haplotype assembly algorithm
[2023:10:26 11:04:40] HIC ITER 0
[2023:10:26 11:04:40] PHASING ITER 0
[2023:10:26 11:04:40] component 0 length 3 phased 2 0...2
[2023:10:26 11:04:40] component 0 offset 0 length 3 phased 2 LL -0.0 cutvalue 13.909938 bestLL -0.03
[2023:10:26 11:04:40] PHASING ITER 1
[2023:10:26 11:04:40] component 0 offset 0 length 3 phased 2 LL -0.0 cutvalue 13.909938 bestLL -0.03
[2023:10:26 11:04:40] PHASING ITER 2
[2023:10:26 11:04:40] component 0 offset 0 length 3 phased 2 LL -0.0 cutvalue 13.909938 bestLL -0.03
[2023:10:26 11:04:40] PHASING ITER 3
[2023:10:26 11:04:40] component 0 offset 0 length 3 phased 2 LL -0.0 cutvalue 13.909938 bestLL -0.03
[2023:10:26 11:04:40] PHASING ITER 4
[2023:10:26 11:04:40] component 0 offset 0 length 3 phased 2 LL -0.0 cutvalue 13.909938 bestLL -0.03
[2023:10:26 11:04:40] PHASING ITER 5
[2023:10:26 11:04:40] component 0 offset 0 length 3 phased 2 LL -0.0 cutvalue 13.909938 bestLL -0.03
[2023:10:26 11:04:40] PHASING ITER 6
[2023:10:26 11:04:40] starting to post-process phased haplotypes to further improve accuracy
[2023:10:26 11:04:40] starting to output phased haplotypes
[2023:10:26 11:04:40] OUTPUTTING PRUNED HAPLOTYPE ASSEMBLY TO FILE haplotype_output_file
[2023:10:26 11:04:40] N50 haplotype length is 0.01 kilobases
[2023:10:26 11:04:40] OUTPUTTING PHASED VCF TO FILE haplotype_output_file.phased.VCF
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele not already represented at this location by REF and ALT">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##INFO=<ID=hapcut2,Number=1,Type=Integer,Description="phased by HapCUT2 or not">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##contig=<ID=chr19,length=61431566>
##FORMAT=<ID=PS,Number=1,Type=Integer,Description="ID of Phase Set for Variant">
##FORMAT=<ID=PQ,Number=1,Type=Integer,Description="Phred QV indicating probability that this variant is incorrectly phased relative to the haplotype">
##FORMAT=<ID=PD,Number=1,Type=Integer,Description="phased Read Depth">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  CGAAGAGGTAGGTGCGAG-1
chr19   55910646        .       AC      A       7447.85 PASS    AC=1    GT:AD:DP:PQ:PD:PS       1|0:47,35:82:100:7:55910646
chr19   55910648        .       AAATCCCCC       A       12412.8 PASS    AC=1    GT:AD:DP:PS     0/1:47,35:82:.
chr19   55910653        .       CCCCAT  C       13372.2 PASS    AC=1    GT:AD:DP:PQ:PD:PS       0|1:0,47:82:100:7:55910646

If I add a second * allele to the third variant to denote the overlapping deletion (as suggested by VCFv4.3) , and add --triallelic 1 I get the same results. Nor do --new_format 1 and --realign_variants 1 change the results.

vibansal commented 1 year ago

Thanks for reporting this. I have not been able to figure out the problem in the first pass. I will have to look deeper.