Issue with haplotagging and ASJ analysis

[Bioman233]

I first ran the following command to haplotag the Iso-Seq alignment BAM:

longcallR -b input.sorted.bam -f Ref.genome.fa -p hifi-masseq -o haplotagged --input-vcf input.dip.vcf.gz --direct-haplotag -t 5

I observed many messages like the following in the log file:
Error: unknown genotype 4 at "chr06":51620118

However, the BAM file was still successfully generated. I then ran:

longcallR asj -a Ref.genome.gff3 -b haplotagged.phased.bam -f Ref.genome.fa -o test.ASJ -t 1

But the output was:
Total genes to be analyzed: 0
Total junctions: 0
number of junctions with at least 10 reads: 0
number of genes with allele-specific junctions: 0

Could you please provide some guidance on:

1. What does the “unknown genotype 4” error indicate?
2. Could this issue affect the haplotagging results?
3. Why does the ASJ analysis return zero genes and junctions?

Any suggestions would be greatly appreciated. Thank you!

[huangnengCSU]

Sorry for late reply.

Regarding "unknown genotype 4":

This is not a critical error, and I will update it to be reported as a warning in future versions. It typically indicates that the GT field contains missing, partial, or otherwise non-standard genotype encodings (e.g., ./., 0/., ./1). These variants are automatically filtered out during processing.

In general, this should not affect the haplotagging results unless all of SNPs have unexpected genotype formats.

Regarding the ASJ analysis returning zero genes and junctions:

There are a few common causes you may need to check:

• Ensure that chromosome names are consistent across all input files (BAM, reference FASTA, and annotation). For example, mismatches like "chr1" vs "1" can prevent proper matching.
• By default, only protein_coding genes and lncRNA are considered. If your annotation file does not include these biotype annotations, no genes will be selected, resulting in empty output. If your annotation file is from Gencode or NCBI, they will have these category tags.
• Please check whether reads in haplotagged.phased.bam contain the HP and PS tags. While some reads may remain unphased (e.g., due to lack of SNP coverage), at least a subset of reads should carry both tags.

[bikekillsme]

I encountered almost the same issue. After checking all the points mentioned above, I found that running the longcallR-ase.py script directly produced normal output, whereas using the software’s built-in function did not.

[huangnengCSU]

Hi @bikekillsme

I’ve updated the code to support both standard and GENCODE/Ensembl GFF3/GTF formats. Could you please compile the latest codes and check again? I know @Bioman233 's issue has already been resolved with the latest codes.

Thanks
Neng

[bikekillsme]

I installed longcallR using mamba, and this issue still persists.

[huangnengCSU]

@bikekillsme

I haven't uploaded to conda and crates.io yet, please clone the source code and compile.

git clone https://github.com/huangnengCSU/longcallR.git
cd longcallR
cargo build --release # The executable will be located at: target/release/longcallR