hiphap

KO | EN

gitlite — search

#javascript #hacktoberfest #ai #llm #deep-learning #docker #nodejs #java #awesome #cli #security #pytorch #openai #open-source

(★ 10)

rust tool to choose the best read alignment between haploids to a diploid reference genome

hiphap Latest Version Download

Download Latest Version (.zip)

// repository documentation

Was this content helpful?

★ 0 (0 ratings)

# HipHap (previously Diplinator)

Diploid genome assemblies are now routinely available, but most read aligners were designed for haploid references. When reads are aligned to a diploid assembly, the aligner sees two nearly identical alignments to either haplotype, and thus reduces the mapping quality (MapQ) score to reflect this ambiguity. This can cause downstream tools to discard reads from easily mappable regions.
HipHap resolves this issue by aligning reads to each haplotype assembly separately and assigning each read to its best-supported haplotype. We also introduce a haplotype assignment quality score (HapQ) in HipHap to quantify confidence in the haplotype of origin of a read.

HipHap is implemented in Rust, and supports SAM, BAM, CRAM, and PAF formats

## Installation

Recommended: Download precompiled binary:
```bash
# Todo
```

Build from source:
```bash
git clone https://github.com/jheinz27/hiphap.git
cd hiphap
cargo build --release
./target/release/hiphap
```

## Usage

```
HipHap: Choose the best alignment to each haploid of a diploid assembly

Usage: hiphap [OPTIONS] <ASM1> <ASM2>

Arguments:
  <ASM1>  asm1 alignment file (sam/bam/cram/paf)
  <ASM2>  asm2 alignment file (sam/bam/cram/paf)

Options:
  -1, --s1 <NAME>          label for asm1 sample (used in output file names and summary) [default: asm1]
  -2, --s2 <NAME>          label for asm2 sample (used in output file names and summary) [default: asm2]
      --paf                input files are PAF
      --ms                 use ms:i: tag rather than AS:i: for alignment score
  -b, --both               write reads with equal alignment scores to both output files
  -m, --merge              write a single merged output file (hiphap_{s1}_{s2}_merged.*) instead of one file per haplotype
      --ref-merged <FILE>  combined reference FASTA for merged CRAM output (must contain all contigs of both inputs); required with --merge on CRAM input
  -u, --unmapped <DEST>    where to write reads unmapped in both assemblies: asm1, asm2, or discard [default: asm1] [possible values: asm1, asm2, discard]
      --ref1 <FILE>        reference FASTA for cram file (asm1)
      --ref2 <FILE>        reference FASTA for cram file (asm2)
  -A, --match-sc <FLOAT>   per-base match score from aligner scoring scheme (auto-estimated from ms:i tags if omitted)
      --no-hapq            skip HAPQ score calculation and hq tag output (e.g. for comparing GRCh38 vs CHM13)
      --no-span-chrom      disable writing hiphap_{s1}_{s2}_span_chrom.fastq
  -t, --threads <INT>      Total thread pool size (min 4). Multiples of 8 recommended for optimal read/write balance. [default: 8]
  -h, --help               Print help
  -V, --version            Print version
```

Each output record is annotated with an `hq:i:` tag carrying the HapQ score, unless `--no-hapq` is set.

The per-base match score used by the HAPQ calculation is auto-estimated from the `ms:i:` tags of the input files. Pass `-A`/`--match-sc <FLOAT>` to set it explicitly (e.g. `2.0` for minimap2 long-read defaults).

## Example Workflow

HipHap only works on name-sorted files, which is the default [minimap2](https://github.com/lh3/minimap2) output. Therefore, coordinate-sorted files need to be name-sorted first, for example:

```bash
samtools sort -n -o name_sort.bam index_sort.bam
```

### Diploid assembly alignment

```bash
# If needed, split diploid genome assembly FASTA into respective haplotypes
separate_haps_fasta -1 MATERNAL -2 PATERNAL hg002v1.1.fa
# writes hg002v1.1.MATERNAL.fa and hg002v1.1.PATERNAL.fa

# Align reads to each haplotype
minimap2 -ax map-hifi -o asm1_alignments.sam hg002v1.1.MATERNAL.fa reads.fastq
minimap2 -ax map-hifi -o asm2_alignments.sam hg002v1.1.PATERNAL.fa reads.fastq

# Run hiphap
hiphap -1 mat -2 pat asm1_alignments.sam asm2_alignments.sam
# Output: hiphap_mat.sam  hiphap_pat.sam  hiphap_mat_pat_span_chrom.fastq

# Merge best alignments into one file (if desired)
samtools merge -@ 12 merged.sam hiphap_mat.sam hiphap_pat.sam

# Save as sorted BAM
samtools sort -@ 12 -o hiphap_mat_pat_merged.bam hiphap_mat_pat_merged.sam
```

By default, HipHap also writes `hiphap_{s1}_{s2}_span_chrom.fastq`, a FASTQ of reads whose alignments span more than one chromosome. These reads are emitted for easy realignment. Pass `--no-span-chrom` to skip writing this file. If inputs are PAF files, a tsv of the chromosome spamnign reads rather than a fastq is output, as the read sequences are not stored in the input PAF files.

### Merged output (`--merge`)

Instead of one file per haplotype, `--merge` writes a single combined output file with a merged header, so the separate `samtools merge` step is no longer needed:

```bash
hiphap --merge -1 mat -2 pat asm1_alignments.sam asm2_alignments.sam
# Output: hiphap_mat_pat_merged.sam  hiphap_mat_pat_span_chrom.fastq

# Save as sorted BAM
samtools sort -@ 12 -o merged.bam hiphap_mat_pat_merged.sam
```

#### Notes:
- In merge mode the **two assemblies must have unique contig names**. The merged header concatenates the two inputs `@SQ` lists, so any contig name shared between them is ambiguous 
- In merge mode HipHap writes only the merged file; the per-haplotype `hiphap_{s1}.*` / `hiphap_{s2}.*` files are not produced.
- `--merge` cannot be combined with `--both`.
- For **CRAM** inputs, merged output is also CRAM and requires a single combined reference (`--ref-merged <FILE>`) containing all contigs of both haplotypes, typically the original diploid assembly:

```bash
  hiphap --merge --ref1 hg002v1.1.MATERNAL.fa --ref2 hg002v1.1.PATERNAL.fa \
         --ref-merged hg002v1.1.fa mat_alignments.cram pat_alignments.cram
  ```

### Comparing different reference genomes

HipHap can also be used to select best alignment of a read between different reference genomes (e.g. GRCh38 and CHM13). For this use case, the HAPQ score is generally not meaningful, so pass `--no-hapq` to skip its calculation:

```bash
hiphap --no-hapq -1 grch38 -2 chm13 grch38_alignments.sam chm13_alignments.sam
# Output: hiphap_grch38.sam  hiphap_chm13.sam
```

### CRAM input files

If input files are CRAM format, the original reference genomes must be provided, for example:

```bash
hiphap --ref1 asm1_hap.fasta --ref2 asm2_hap.fasta asm1_alignments.cram asm2_alignments.cram
```

## Example PAF Usage

#### Notes: 
- It is important to use the `--paf-no-hit` flags when aligning with minimap2 to output unmapped reads to the file
- If a SAM file is converted to a PAF file with `paftools.js sam2paf`, it will **NOT** have the required AS:i: tag and HipHap will fail to run

```bash
minimap2 -cx map-hifi --paf-no-hit -o asm1_alignments.paf hg002v1.1.MATERNAL.fa reads.fastq
minimap2 -cx map-hifi --paf-no-hit -o asm2_alignments.paf hg002v1.1.PATERNAL.fa reads.fastq

hiphap --paf out_asm1.paf out_asm2.paf
# Output: hiphap_asm1.paf  hiphap_asm2.paf  hiphap_asm1_asm2_span_chrom.txt
```

`--merge` also works with PAF

```bash
hiphap --paf --merge out_asm1.paf out_asm2.paf
# Output: hiphap_asm1_asm2_merged.paf
```

## Weighted Alignment Scoring Mechanism

For each read, HipHap computes a single weighted alignment score per assembly using all primary and supplementary alignments (secondary alignments are passed through to the output but ignored when scoring).

Let $n$ be the number of primary and supplementary alignments for a given read to that reference genome.

Let $L$ be the full read length (sum of query-consuming CIGAR operations on a non-secondary record, or the `qlen` field of the PAF record).

Let $I_i = [r_i^{\mathrm{start}}, r_i^{\mathrm{end}})$ denote the interval on the read covered by alignment $i$.

Let

$$
B = \left| \bigcup_{i=1}^{n} I_i \right|
$$

be the number of **read bases** covered by at least one alignment.

Let $a_i$ be the alignment score for alignment $i$ (`AS:i:` by default, or `ms:i:` if `--ms` is set).

Let $l_i$ be the alignment length in read coordinates for alignment $i$.

$$
S = \frac{\sum_{i=1}^{n} a_i}{\sum_{i=1}^{n} l_i} \cdot B \cdot \frac{B}{L}
$$

The first factor is the average alignment score per aligned base. It is multiplied by the number of unique read bases covered, $B$, and then scaled by the read coverage fraction $B/L$, so that reads which align over a large fraction of their length are weighted more heavily than reads which align only over a small portion.

For each read, the assembly with the higher $S$ wins; its full alignment cluster (including secondary alignments) is written to the corresponding output file. If $S$ is equal in both assemblies, the "better" assignment is determined by a hash of the read name, or the read is written to both output files when `--both` is used.

## Haplotype Assignment Quality (HapQ)

For each read assigned to a winning haplotype, HipHap reports a HapQ score in the `hq:i:` tag of the output record. HapQ is a Phred-like confidence [0-60] that the read was assigned to the correct haplotype. The calculation is modeled on BWA-MEM's [`mem_approx_mapq_se`](https://github.com/lh3/bwa/blob/b92993c1161e73167181558856567ef2f367e3f0/bwamem.c#L982-L1006).

Let $S_w$ and $S_l$ be the weighted alignment scores of the winning and losing assemblies, $m$ the per-base match score of the alignment software used (`-A`/`--match-sc`), and $k$ the number of non-secondary alignments (splits) on the winning side.

HapQ is the product of:

d, the approximate difference, in matching bases, between the winning and losing alignments between haplotypes.
```math
d = \frac{S_w - S_l}{m}
```

The split penalty, $\rho$ , which penalizes reads with more than three split alignments, which often fall in complex/repetitive regions where haplotype assignment is less reliable.
```math
\rho = \begin{cases} 1 & k \le 3 \\ \frac{3}{k} & k > 3 \end{cases}
```

The final HapQ calculation is then: 
```math
\text{HapQ}\;=\;\lfloor\, 6.02 \cdot \rho \cdot d \,\rfloor
```

Clamped to the range [0,60].

Special cases:
- Read mapped in only one assembly: HapQ = 60.
- Read tied between assemblies (winner = `Both`): HapQ = 0.
- Read unmapped in both assemblies: no `hq` tag is written.

If `--no-hapq` is set, HAPQ is not computed and no `hq:i:` tag is added (recommended when the two inputs are not haplotypes of the same sample, e.g. GRCh38 vs CHM13).

## Citation
If HipHap has helped you in your research, please cite our preprint at: TODO

hiphap

hiphap Latest Version Download

Today's picks

Was this content helpful?

Recent Feedback

hiphap

hiphap Latest Version Download

Today's picks

Was this content helpful?

Recent Feedback

Download README