HipSTR: genotyping forensic STRs from a BAM

Before running anything, look at what you have. Here, our example: the reads (.bam) with their index (.bai), the list of STRs (.bed) with its table of allele numbers (.tsv), the reference (.fa.gz, with its indexes .fai and .gzi), and the two folders from the setup: HipSTR/ and hipstr-env/. With your data you need the same: your BAM or CRAM files with their indexes, your reference with its index, and your BED.

HipSTR takes the sample name from the read group of the BAM (the @RG line): SM is the sample and LB the library. If your BAM has no read group, HipSTR does not start.

Each line is one STR. The columns, as the HipSTR README defines them: chromosome, start, end, motif length (4 = tetranucleotide), number of copies of the motif in the reference allele, and (optional) the locus name, which HipSTR writes in the ID column of the VCF. The BED has no header line.

Our BED uses the ISFG minimum ranges (STRidER FSSG v6-1, Recommendation 3 of Gettings et al. 2024): the stretch of each STR that forensic nomenclature is reported over. They are a little wider than the bare repeat, because they include the bases that the forensic allele number counts. That is what lets you turn HipSTR's results into forensic allele numbers in step 8. For your forensic STRs you can use the same BED: it is in strhub-demo-data.

This is the basic command of the HipSTR README (--bams, --fasta, --regions, --str-vcf) plus --viz-out, which the official HipSTR tutorial also uses and which you need for VizAln (step 12). The log goes to the screen and tee keeps a copy in NA12878.log.

The summary at the end says that only 14 of the 22 loci were genotyped, and nothing else on screen warns you:

• 2 loci skipped for length: D2S1338 (101 bp in its ISFG range) and D21S11 (138 bp) are longer than --max-str-len (default 100 bp).
• 6 loci skipped for too few reads to learn the stutter model: D1S1656, D12S391, D13S317, D19S433, Penta D and FGA (in the log: TOTAL=65 to 97, MIN=100). HipSTR learns a stutter model per locus from the reads of all samples; with one sample this may not be enough. The ISFG ranges are wider than the bare repeat, so fewer reads cross them completely.

The log suggests --min-reads, but for very few samples the HipSTR README recommends --def-stutter-model.

The table below says, for each option, what you put there with your own data.

Three options are added at the end, the rest stays the same. The summary now says 22 of 22 (careful: in step 10 you will see that one of them, D21S11, does not pass the filter).

HipSTR writes its result to NA12878.vcf.gz: a VCF (Variant Call Format) file compressed with bgzip. gunzip -k decompresses it and keeps the .gz, which the next steps read (-f overwrites NA12878.vcf if you run it again). Without -k, the .gz disappears and the next commands (bcftools and the filter) stop working.

The lines that start with ## are the header, and they define every field (try grep '^##FORMAT' NA12878.vcf). Then comes one line per STR, with tab-separated columns: CHROM, POS, ID (the locus name), REF and ALT (the allele sequences), QUAL, FILTER, INFO (data for the locus), FORMAT (the names of the per-sample fields, separated by :) and one column per sample, with the values in the same order.

To open it in a spreadsheet (Excel, LibreOffice Calc, Google Sheets), open it from inside the program and choose tab as the separator. Do not double-click it: .vcf is also the extension of contact cards (vCard), and your computer may open your contacts app instead. To see each field next to its value, split the FORMAT and NA12878 columns at : (menu Data, Text to Columns).

The fields below are the ones in the FORMAT column. The descriptions come from the VCF header and from the HipSTR README.

bcftools query (bcftools is not part of HipSTR: it comes from the samtools team) prints only the fields you ask for:

• GT: which sequences were called (indexes into REF/ALT, not alleles).
• GB: difference in bp of each allele from the reference. 0|4 at D13S317 = one allele as long as the reference, the other 4 bp longer. When GB is not a multiple of the motif, it is a microvariant: TH01 (motif of 4) shows 0|11 = 2 repeats plus 3 bases.
• DP: valid reads HipSTR used (the total coverage). Q: posterior probability of the genotype.
• FILTER: PASS, or why the call was filtered (thanks to --output-filters).

All 22 loci have a genotype. To turn GB into a forensic allele number you need the allele number of the reference at each locus: that is step 8. Careful with D21S11 (0|0): HipSTR gives it a genotype, but in step 10 the filter removes it.

HipSTR reports each allele as a length difference from the reference (GB). Forensic laboratories report allele numbers (the CE designation). Because our BED uses the ISFG minimum ranges, the conversion is a sum: allele = ce_ref + GB / motif length. ce_ref is the allele number of the GRCh38 reference over that range: column 6 of isfg_ranges_grch38.tsv, from the STRidER FSSG v6-1 (Gettings et al. 2024, FSI: Genetics 68:102946).

D13S317: ce_ref 11, GB 0|4, motif of 4 bp: alleles 11 and 12. When the difference is not a whole number of repeats, the remainder goes after a point: TH01 (ce_ref 7, motif 4) with GB 11 is 7 + 2 repeats + 3 bases = 9.3.

With other coordinates (for example the repeat alone, as in HipSTR's own BED files) the same sum gives numbers that do not match the forensic designation in several loci. That is why the BED matters.

MALLREADS (from the HipSTR VCF) gives, for each length difference in bp, how many reads cross the repeat with at least 5 bp on each side. At D13S317: 47 reads at 0 bp and 26 at +4 bp, about 64/36: a balanced heterozygote. The VCF field AB tests this same balance (here -0.42, close to 0 = balanced). There are fewer reads than DP (413) because the ISFG range is wider than the bare repeat, and fewer reads cross all of it with 5 bp to spare.

HipSTR comes with scripts/filter_vcf.py, and its README suggests these thresholds. The filtered VCF keeps every locus but removes, sample by sample, the calls that fail, and writes the reason in FILTER. The second command lists every call that is not PASS. The script runs with the Python of hipstr-env, which has PyVCF3.

Here D21S11 is removed: 45 % of its reads have a stutter-like indel (the limit is 15 %). Its ISFG range is 138 bp and the reads are 148 bp, so almost no read crosses it with room to spare: HipSTR still writes a genotype (0|0), but it is not reliable. A genotype in the VCF is not always a reliable genotype. Every other call passes (Q = 1, stutter ≤ 2.5 %, allele bias ≥ -1.44): this sample has ~300x.

Only to see what a removed call looks like: with a 2 % stutter threshold, D2S1338 (2.5 % of its reads with stutter) is removed too; D21S11 was already out. The genotype becomes ./. and FILTER says STUTTER_FRACTION. Use the README thresholds for real data.

--viz-out (step 4) saved the alignments. The first command indexes them with tabix -p bed, as the HipSTR README says. The second activates the environment of the setup, because VizAln calls python. The third runs VizAln with the file, the chromosome and start of the locus (columns 1 and 2 of its line in your BED: for D13S317, chr13 82148020) and the sample.

VizAln writes an HTML file in /tmp and tries to open it in your browser. On macOS it usually opens; on WSL or a server there is no browser, so open the file the message names: the reference on top, the call in red (NA12878: 0|4) and one line per read. On macOS, VizAln also needs GNU cut (see the errors below).