README.md
1# Graph-Based Genome Alignment and Genotyping with HISAT2 and HISAT-genotype
2
3## Contact
4
5[Daehwan Kim](https://kim-lab.org) (infphilo@gmail.com) and [Chanhee Park](https://www.linkedin.com/in/chanhee-park-97677297/) (parkchanhee@gmail.com)
6
7## Abstract
8
9Rapid advances in next-generation sequencing technologies have dramatically changed our ability to perform genome-scale analyses. The human reference genome used for most genomic analyses represents only a small number of individuals, limiting its usefulness for genotyping. We designed a novel method, HISAT2, for representing and searching an expanded model of the human reference genome, in which a large catalogue of known genomic variants and haplotypes is incorporated into the data structure used for searching and alignment. This strategy for representing a population of genomes, along with a fast and memory-efficient search algorithm, enables more detailed and accurate variant analyses than previous methods. We demonstrate two initial applications of HISAT2: HLA typing, a critical need in human organ transplantation, and DNA fingerprinting, widely used in forensics. These applications are part of HISAT-genotype, with performance not only surpassing earlier computational methods, but matching or exceeding the accuracy of laboratory-based assays.
10
11
12
13For more information, see the following websites:
14* [HISAT2 website](http://ccb.jhu.edu/software/hisat2)
15* [HISAT-genotype website](http://ccb.jhu.edu/software/hisat-genotype)
16
17## HISAT2
18HISAT2 is a fast and sensitive alignment program for mapping next-generation sequencing reads (whole-genome, transcriptome, and exome sequencing data) to a population of human genomes (as well as to a single reference genome). Based on an extension of BWT for a graph [1], we designed and implemented a graph FM index (GFM), an original approach and its first implementation to the best of our knowledge. In addition to using one global GFM index that represents general population, HISAT2 uses a large set of small GFM indexes that collectively cover the whole genome (each index representing a genomic region of 56 Kbp, with 55,000 indexes needed to cover human population). These small indexes (called local indexes) combined with several alignment strategies enable effective alignment of sequencing reads. This new indexing scheme is called Hierarchical Graph FM index (HGFM). We have developed HISAT2 based on the HISAT [2] and Bowtie 2 [3] implementations. See the [HISAT2 website](http://ccb.jhu.edu/software/hisat2/index.shtml) for
19more information.
20
21A few notes:
22
231) HISAT2's index (HGFM) size for the human reference genome and 12.3 million common SNPs is 6.2GB. The SNPs consist of 11 million single nucleotide polymorphisms, 728,000 deletions, and 555,000 insertions. Insertions and deletions used in this index are small (usually <20bp). We plan to incorporate structural variations (SV) into this index.
24
252) The memory footprint of HISAT2 is relatively low, 6.7GB.
26
273) The runtime of HISAT2 is estimated to be slightly slower than HISAT (30–100% slower for some data sets).
28
294) HISAT2 provides greater accuracy for alignment of reads containing SNPs.
30
315) We released a first (beta) version of HISAT2 in September 8, 2015.
32
33## License
34
35[GPL-3.0](LICENSE)
36
37# For reviwers, follow the instructions below to reproduce some of the results in the manuscript.
38
39## Code
40
41A specific version of [HISAT2 and HISAT-genotype](http://github.com/infphilo/hisat2) at GitHub is used (a branch name: hisat2_v2.2.0_beta).
42
43## Initial setup
44
45HISAT-genotype requires a 64-bit computer running either Linux or Mac OS X and at least 8 GB of RAM (16 GB of RAM is preferred). All the commands used should be run from the Unix shell prompt within a terminal window and are prefixed with a '$' character.
46
47We refer to <b>hisat-genotype-top</b> as our top directory where all of our programs are located. <b>hisat-genotype-top</b> is a place holder that can be changed to another name according to user preference.
48Run the following commands to install HISAT2 and HISAT-genotype.
49
50 $ git clone https://github.com/infphilo/hisat2 hisat-genotype-top
51 $ cd hisat-genotype-top
52 hisat-genotype-top$ git checkout hisat2_v2.2.0_beta
53 hisat-genotype-top$ make hisat2-align-s hisat2-build-s hisat2-inspect-s
54
55To make the binaries built above and other python scripts available everywhere, add the hisat-genotype-top directory to the PATH environment variable (e.g. ~/.bashrc)
56
57 export PATH=hisat-genotype-top:hisat-genotype-top/hisatgenotype_scripts:$PATH
58 export PYTHONPATH=hisat-genotype-top/hisatgenotype_modules:$PYTHONPATH
59
60To reflect the change, run the following command:
61
62 $ source ~/.bashrc
63
64Download real reads, simulated reads, and HISAT2 indexes, then move them into appropriate directories:
65
66 hisat-genotype-top$ cd evaluation
67 hisat-genotype-top/evaluation$ wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hisat2_20181025.tar.gz
68 hisat-genotype-top/evaluation$ tar xvzf hisat2_20181025.tar.gz
69 hisat-genotype-top/evaluation$ mkdir aligners aligners/bin; cd aligners/bin; ln -s ../../../hisat2* .; cd ../..
70 hisat-genotype-top/evaluation$ mv hisat2/* .
71 hisat-genotype-top/evaluation$ cd simulation; ./init.py; cd ../real; ./init.py; cd ..
72
73## Run HISAT2 on the following simulated and real data sets.
74### 10 million simulated read pairs with SNPs and with sequencing errors
75
76 hisat-genotype-top/evaluation$ cd simulation/10M_DNA_mismatch_snp_reads_genome
77 hisat-genotype-top/evaluation/simulation/10M_DNA_mismatch_snp_reads_genome$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
78
79### 10 million simulated read pairs with SNPs and without sequencing errors
80
81 hisat-genotype-top/evaluation$ cd simulation/10M_DNA_snp_reads_genome
82 hisat-genotype-top/evaluation/simulation/10M_DNA_snp_reads_genome$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
83
84### 10 million simulated read pairs without SNPs and with sequencing errors
85 hisat-genotype-top/evaluation$ cd simulation/10M_DNA_mismatch_reads_genome
86 hisat-genotype-top/evaluation/simulation/10M_DNA_mismatch_reads_genome$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
87
88### 10 million simulated read pairs without SNPs and without sequencing errors
89 hisat-genotype-top/evaluation$ cd simulation/10M_DNA_reads_genome
90 hisat-genotype-top/evaluation/simulation/10M_DNA_reads_genome$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
91
92### 10 million real read pairs
93 hisat-genotype-top/evaluation$ cd real/DNA/10M
94 hisat-genotype-top/evaluation/real/DNA/10M$ ./calculate_read_cost.py --aligner-list hisat2 --paired-end --fresh
95
96### Interpreting output
97 Example alignment output for simulated reads
98 aligned: 1000000, multi aligned: 2654390
99 correctly mapped: 999963 (100.00%)
100 uniquely and correctly mapped: 967631 (96.76%)
101 54694 reads per sec (all)
102 Memory Usage: 86MB
103
104The above lines show that 1,000,000 read pairs are aligned and the total number of alignments is 2,654,390. 999,963 pairs (100.00%) are correctly aligned (e.g. one of the alignments is correct). 967,631 (96.76%) pairs are uniquely and correctly aligned. HISAT2 aligns 54,594 reads with a peak memory usage of 86 MB of RAM.
105
106Each run is expected to take up to several hours mostly due to the comparison of HISAT2’s reported alignments and true alignments and the expansion of repeat alignments.
107
108## Details on HISAT-genotype run for HLA typing and assembly
109
110To create a directory where we perform our analysis for HLA typing and assembly, which here is referred to as hla-analysis but can be changed by the user, execute the following command.
111
112 hisat-genotype-top/evaluation$ mkdir hla-analysis
113
114The current directory can be changed to hla-analysis as follows:
115
116 hisat-genotype-top/evaluation$ cd hla-analysis
117
118Additional program requirements: SAMtools (version 1.3 or later)
119
120### Downloading a Graph Reference and Index
121The graph reference we are going to build incorporates variants of numerous HLA alleles into the linear reference using a graph. The graph reference also includes some known variants of other regions of the genome (e.g. common small variants). To copy the graph reference, type:
122
123 hisat-genotype-top/evaluation/hla-analysis$ mv ../hisat2-genotype/* .
124
125### Typing and Assembly
126Since whole genome sequencing (WGS) data includes reads that are from the whole genome, the first step is to extract the reads that belong to the HLA genes by aligning them to the graph reference with HISAT2. We provide these extracted reads in hisat-genotype-top/evaluation/hla-analysis/ILMN_20181025.
127
128HISAT-genotype performs both HLA typing and assembly as follows.
129You can perform HLA typing and assembly for HLA-A gene on sequencing reads from the genome NA12892 (Illumina's HiSeq 2000 platform).
130
131 hisat-genotype-top/evaluation/hla-analysis$ hisatgenotype_locus.py --base hla --locus-list A --assembly -1 ILMN_20181025/NA12892.hla.extracted.1.fq.gz -2 ILMN_20181025/NA12892.hla.extracted.2.fq.gz
132
133### DNA Fingerprinting
134This function can be performed with the same commands used for “Typing and Assembly” and just replacing --base hla with --base codis.
135
136### Interpreting Output
137 Typing Output
138 Number of reads aligned: 1507
139 1 A*02:01:01:02L (count: 571)
140 2 A*02:01:31 (count: 557)
141 3 A*02:20:02 (count: 557)
142 4 A*02:29 (count: 557)
143 5 A*02:321N (count: 556)
144 6 A*02:372 (count: 556)
145 7 A*02:610:02 (count: 556)
146 8 A*02:249 (count: 555)
147 9 A*02:479 (count: 555)
148 10 A*02:11:01 (count: 554)
149
150The above lines show the top ten alleles that the most number of reads are mapped to or compatible with. For example, the allele first ranked, A\*02:01:01:02L, is compatible with 571 reads. This raw estimate based on the number of reads should not be used to determine the two true alleles because the alleles that resemble both but are not true alleles often tend to be compatible with more reads than either of the true alleles. Thus, we apply a statistical model to identify the two true alleles as described in the main text.
151
152 Abundance of alleles
153 1 ranked A*02:01:01:01 (abundance: 54.32%)
154 2 ranked A*11:01:01:01 (abundance: 45.20%)
155 3 ranked A*24:33 (abundance: 0.48%)
156
157The above rankings show the top three alleles that are most abundant in the sample. Normally, the top two alleles in this estimate (e.g. A\*02:01:01:01 and A\*11:01:01:01) are considered as the two alleles that best match a given sequencing data.
158
159Additional tutorials and details are available at the HISAT-genotype website: https://ccb.jhu.edu/hisat-genotype
160
161
162## Data
163
164The Data directory (`/data`) contains all input files for reproducing some of our results such as from the
165evaluation of HISAT2 and other programs using both simulated and real reads, from typing and assembling
166HLA genes of Illumina Platinum Genomes using HISAT-genotype, and from building a HISAT2 graph index.
167
168* **Simulated read pairs**
169
170| Type | Number of pairs | Path |
171| - | - | - |
172| SNPs and sequencing errors included | 10,000,000 | hisat-genotype-top/evaluation/reads/simulation/10M_DNA_mismatch_snp_reads_genome |
173| SNPs included | 10,000,000 | hisat-genotype-top/evaluation/reads/simulation/10M_DNA_snp_reads_genome |
174| Sequencing errors included | 10,000,000 | hisat-genotype-top/evaluation/reads/simulation/10M_DNA_mismatch_reads_genome |
175| No SNPs nor sequencing errors included | 10,000,000 | hisat-genotype-top/evaluation/reads/simulation/10M_DNA_reads_genome |
176
177Each directory comes with a true alignment file in SAM format so that users know where the reads were generated in the human reference genome.
178
179* **Real read pairs**
180
181| Number of read pairs | Path |
182| - | - |
183| 10,000,000 | hisat-genotype-top/evaluation/reads/real/DNA |
184
185* **Human reference genome, SNPs, haplotypes, and HISAT2's indexes**
186
187| Type | Path |
188| - | - |
189| GRCh38 reference | hisat-genotype-top/evaluation/data/genome.fa |
190| SNPs | hisat-genotype-top/evaluation/data/genome.snp |
191| Haplotypes | hisat-genotype-top/evaluation/data/genome.haplotype |
192| HISAT2's prebuilt graph index for comparison with other aligners | hisat-genotype-top/evaluation/indexes/HISAT2/genome.[1-8].ht2 |
193| HISAT2's prebuilt graph index for genotyping | hisat-genotype-top/evaluation/hla-analysis/genotype_genome.[1-8].ht2 |
194
195
196## References
197
198[1] Sirén J, Välimäki N, Mäkinen V (2014) Indexing graphs for path queries with applications in genome research. IEEE/ACM Transactions on Computational Biology and Bioinformatics 11: 375–388. doi: 10.1109/tcbb.2013.2297101
199
200[2] Kim D, Langmead B, and Salzberg SL HISAT: a fast spliced aligner with low memory requirements, Nature methods, 2015
201
202[3] Langmead B, Salzberg SL: Fast gapped-read alignment with Bowtie 2. Nat Methods 2012, 9:357-359
203
