Skip to content
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

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

Open
jan-glx opened this issue Oct 26, 2023 · 1 comment

Comments

@jan-glx
Copy link

jan-glx commented Oct 26, 2023

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 0eb7075 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
Copy link
Owner

vibansal commented Nov 8, 2023

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants