-
Notifications
You must be signed in to change notification settings - Fork 0
Description
I am testing whether I can get (partial) haplotype reconstructed from Illumina short reads derived from a diploid heterozygous tree genome.
I generated a vcf file with freeBayes after mapping the reads to a reference genome for just one scaffold, 5.5 Mb long. pHapCompass can generate the .frag file within a few seconds, but when it comes to the actual haplotyping I get an error message.
2026-01-25 21:33:27,904 | INFO | Inferring ploidy from VCF: /home/dfbs_user/Results/Tg2_scaffold48_FreeBayes.vcf
2026-01-25 21:33:27,904 | INFO | Inferred ploidy: 2
2026-01-25 21:33:27,912 | INFO | No fragment file provided. Will generate using extractHAIRS at: /home/dfbs_user/Results/Tg2_vs_TgS01_scaffold48_sorted.frag
Using extractHAIRS binary: /home/dfbs_user/pHapCompass/src/phapcompass/bin/extractHAIRS
Running: /home/dfbs_user/pHapCompass/src/phapcompass/bin/extractHAIRS --bam /home/dfbs_user/Results/Tg2_vs_TgS01_scaffold48_sorted.bam --VCF /home/dfbs_user/Results/Tg2_scaffold48_FreeBayes.vcf --out /home/dfbs_user/Results/Tg2_vs_TgS01_scaffold48_sorted.frag --mbq 13
Extracting haplotype informative reads from bamfiles /home/dfbs_user/Results/Tg2_vs_TgS01_scaffold48_sorted.bam minQV 13 minMQ 20 maxIS 1000
VCF file /home/dfbs_user/Results/Tg2_scaffold48_FreeBayes.vcf has 1686 variants
adding chrom scaffold48_cov80 to index
vcffile /home/dfbs_user/Results/Tg2_scaffold48_FreeBayes.vcf chromosomes 1 hetvariants 1082 1686
reading sorted bamfile /home/dfbs_user/Results/Tg2_vs_TgS01_scaffold48_sorted.bam
processing reads mapped to chrom "scaffold48_cov80"
final cleanup of fragment list: 11618 current chrom (null) 0 prev 0
✓ Fragment file created successfully: /home/dfbs_user/Results/Tg2_vs_TgS01_scaffold48_sorted.frag
2026-01-25 21:33:29,183 | INFO | Data type: short
2026-01-25 21:33:29,183 | INFO | Result path: /home/dfbs_user/Results/Tg2_scaffold48_haplotypes
2026-01-25 21:33:29,183 | INFO | Running short-read model with mw=10.0000, lw=1.0000, sw=1.0000
Traceback (most recent call last):
File "/home/dfbs_user/conda/bin/phapcompass", line 7, in
sys.exit(main())
^^^^^^
File "/home/dfbs_user/pHapCompass/src/phapcompass/cli.py", line 126, in main
run_pHapCompass_short(args)
File "/home/dfbs_user/pHapCompass/src/phapcompass/run_pHapCompass_short.py", line 54, in run_pHapCompass_short
gi = g_per_snp[col_to_snp[nodes_cols[:, 0]]] # g at SNP position for col_i
~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IndexError: index 1681 is out of bounds for axis 0 with size 1681
OS Linux Ubuntu 24.04.3 LTS
Python 3.12
wc -l ./Results/Tg2_vs_TgS01_scaffold48_sorted.frag
3763 ./Results/Tg2_vs_TgS01_scaffold48_sorted.frag