2018-03-21

高速なRNA seqのマッピングツール HISAT2

2019 6/18 コマンド追記、6/26 インストール追記、6/28　samtoolsコマンドエラー修正

2020 3/22 help更新、4/16 multiqcとの連携例、4/29 誤解のある表現を修正、8/28 index追記, help更新

2021 1/21、8/26 インストール追記

2024/02/16 インストール手順（pythonバージョン）とhelp更新

　RNA-seqは、2008年に導入されて以来、遺伝子発現、転写体構造、長い非コード化RNAと融合転写物の同定のためのツールとして普及してきた（論文より ref.2-5） RNA-seq解析は、リードを参照ゲノムに対してアライメントさせ、リードの起点となる場所を決定するが、時間のかかるボトルネックステップになっている。例えば、TopHat2（ref.9）およびGSNAP(ref.10)のような広く使用されているアラインメントプログラムは、単一のRNA-seq実験を処理するのに数日かかることがある。最近導入されたSTARプログラム(ref.11)は、サフィックス配列を使用して、TopHat2を含む他のほとんどの方法よりも大幅に高速な処理を提供する（STAR紹介）。しかしながら、STARのアルゴリズムはBurrows-Wheeler変換を用いる方法と比較して、非常に大きなメモリ要件（ヒトゲノムについては28ギガバイト（GB））を有する。

　著者らは、適度な量のランダムアクセスメモリ（RAM）を使用する高速スプライスアライナを作成するために、Burrows-Wheelerトランスフォーム（ref.12）とFMインデックス（ref.13）に基づく新しいインデックス作成戦略を使用してHISATを設計した。HISAT（hierarchical indexing for spliced alignment of transcripts）は、Burrows-WheelerトランスフォームとFerragina-Manzini（FM）インデックスに基づいたインデックス作成スキームを使用している。ヒトゲノムのHISATの階層インデックスは、48,000個のローカルFMインデックスを含み、それぞれが約64,000bpのゲノム領域を表す。実際のデータセットとシミュレートされたデータセットのテストでは、HISATが現在利用可能な最速のシステムであることが示されている（論文執筆時点）。 HISATはインデックス数が多いにもかかわらず、わずか4.3ギガバイトのメモリしか必要としない。 HISATは、40億塩基以上のゲノムを含め、あらゆるサイズのゲノムをサポートしている。HISATのメモリ要件が大幅に削減された結果、ユーザはこれらの計算を専用サーバから従来の単一のデスクトップコンピュータで実行できる。複数のサンプルを同時に実行することも可能とされる（TopHatの開発者として、HISATはそのプログラムの次のメジャーバージョンであるTopHat3のコアアラインメントエンジンにする予定と記載されている）。

マニュアルにはHISAT2はHISATとtophat2の後継とある（HISAT2 is a successor to both HISAT and TopHat2. We recommend that the HISAT and TopHat2 users switch to HISAT2.）。ここではHISAT2を紹介する。

公式ページ

マニュアル

http://ccb.jhu.edu/software/hisat2/manual.shtml

インストール

Github

brewやcondaで導入できる。

#bioconda(linux) python環境を指定して仮想環境に入れる
mamba create -n hisat2 -y python=3 -y
conda activate hisat2
mamba install -c bioconda hisat2 -y

#homebrew
#tapしてなければtapしておく
brew tap brewsci/bio
brew tap brewsci/science

brew install hisat2 bwa #brewで両方入れる。

Bioconda環境下ならcondaでも導入できます。

> hisat2-build #v2.2.1

$ hisat2-build

No input sequence or sequence file specified!

HISAT2 version 2.2.1 by Daehwan Kim (infphilo@gmail.com, http://www.ccb.jhu.edu/people/infphilo)

Usage: hisat2-build [options]* <reference_in> <ht2_index_base>

reference_in comma-separated list of files with ref sequences

hisat2_index_base write ht2 data to files with this dir/basename

Options:

-c reference sequences given on cmd line (as

<reference_in>)

--large-index force generated index to be 'large', even if ref

has fewer than 4 billion nucleotides

-a/--noauto disable automatic -p/--bmax/--dcv memory-fitting

-p <int> number of threads

--bmax <int> max bucket sz for blockwise suffix-array builder

--bmaxdivn <int> max bucket sz as divisor of ref len (default: 4)

--dcv <int> diff-cover period for blockwise (default: 1024)

--nodc disable diff-cover (algorithm becomes quadratic)

-r/--noref don't build .3/.4.ht2 (packed reference) portion

-3/--justref just build .3/.4.ht2 (packed reference) portion

-o/--offrate <int> SA is sampled every 2^offRate BWT chars (default: 5)

-t/--ftabchars <int> # of chars consumed in initial lookup (default: 10)

--localoffrate <int> SA (local) is sampled every 2^offRate BWT chars (default: 3)

--localftabchars <int> # of chars consumed in initial lookup in a local index (default: 6)

--snp <path> SNP file name

--haplotype <path> haplotype file name

--ss <path> Splice site file name

--exon <path> Exon file name

--repeat-ref <path> Repeat reference file name

--repeat-info <path> Repeat information file name

--repeat-snp <path> Repeat snp file name

--repeat-haplotype <path> Repeat haplotype file name

--seed <int> seed for random number generator

-q/--quiet disable verbose output (for debugging)

-h/--help print detailed description of tool and its options

--usage print this usage message

--version print version information and quit

> hisat2 -h #v2.2.1

$ hisat2

No index, query, or output file specified!

HISAT2 version 2.2.1 by Daehwan Kim (infphilo@gmail.com, www.ccb.jhu.edu/people/infphilo)

Usage:

hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]

<ht2-idx> Index filename prefix (minus trailing .X.ht2).

<m1> Files with #1 mates, paired with files in <m2>.

Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).

<m2> Files with #2 mates, paired with files in <m1>.

Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).

<r> Files with unpaired reads.

Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).

<sam> File for SAM output (default: stdout)

<m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be

specified many times. E.g. '-U file1.fq,file2.fq -U file3.fq'.

Options (defaults in parentheses):

Input:

-q query input files are FASTQ .fq/.fastq (default)

--qseq query input files are in Illumina's qseq format

-f query input files are (multi-)FASTA .fa/.mfa

-r query input files are raw one-sequence-per-line

-c <m1>, <m2>, <r> are sequences themselves, not files

-s/--skip <int> skip the first <int> reads/pairs in the input (none)

-u/--upto <int> stop after first <int> reads/pairs (no limit)

-5/--trim5 <int> trim <int> bases from 5'/left end of reads (0)

-3/--trim3 <int> trim <int> bases from 3'/right end of reads (0)

--phred33 qualities are Phred+33 (default)

--phred64 qualities are Phred+64

--int-quals qualities encoded as space-delimited integers

Presets: Same as:

--fast --no-repeat-index

--sensitive --bowtie2-dp 1 -k 30 --score-min L,0,-0.5

--very-sensitive --bowtie2-dp 2 -k 50 --score-min L,0,-1

Alignment:

--bowtie2-dp <int> use Bowtie2's dynamic programming alignment algorithm (0) - 0: no dynamic programming, 1: conditional dynamic programming, and 2: unconditional dynamic programming (slowest)

--n-ceil <func> func for max # non-A/C/G/Ts permitted in aln (L,0,0.15)

--ignore-quals treat all quality values as 30 on Phred scale (off)

--nofw do not align forward (original) version of read (off)

--norc do not align reverse-complement version of read (off)

--no-repeat-index do not use repeat index

Spliced Alignment:

--pen-cansplice <int> penalty for a canonical splice site (0)

--pen-noncansplice <int> penalty for a non-canonical splice site (12)

--pen-canintronlen <func> penalty for long introns (G,-8,1) with canonical splice sites

--pen-noncanintronlen <func> penalty for long introns (G,-8,1) with noncanonical splice sites

--min-intronlen <int> minimum intron length (20)

--max-intronlen <int> maximum intron length (500000)

--known-splicesite-infile <path> provide a list of known splice sites

--novel-splicesite-outfile <path> report a list of splice sites

--novel-splicesite-infile <path> provide a list of novel splice sites

--no-temp-splicesite disable the use of splice sites found

--no-spliced-alignment disable spliced alignment

--rna-strandness <string> specify strand-specific information (unstranded)

--tmo reports only those alignments within known transcriptome

--dta reports alignments tailored for transcript assemblers

--dta-cufflinks reports alignments tailored specifically for cufflinks

--avoid-pseudogene tries to avoid aligning reads to pseudogenes (experimental option)

--no-templatelen-adjustment disables template length adjustment for RNA-seq reads

Scoring:

--mp <int>,<int> max and min penalties for mismatch; lower qual = lower penalty <6,2>

--sp <int>,<int> max and min penalties for soft-clipping; lower qual = lower penalty <2,1>

--no-softclip no soft-clipping

--np <int> penalty for non-A/C/G/Ts in read/ref (1)

--rdg <int>,<int> read gap open, extend penalties (5,3)

--rfg <int>,<int> reference gap open, extend penalties (5,3)

--score-min <func> min acceptable alignment score w/r/t read length

(L,0.0,-0.2)

Reporting:

-k <int> It searches for at most <int> distinct, primary alignments for each read. Primary alignments mean

alignments whose alignment score is equal to or higher than any other alignments. The search terminates

when it cannot find more distinct valid alignments, or when it finds <int>, whichever happens first.

The alignment score for a paired-end alignment equals the sum of the alignment scores of

the individual mates. Each reported read or pair alignment beyond the first has the SAM ‘secondary’ bit

(which equals 256) set in its FLAGS field. For reads that have more than <int> distinct,

valid alignments, hisat2 does not guarantee that the <int> alignments reported are the best possible

in terms of alignment score. Default: 5 (linear index) or 10 (graph index).

Note: HISAT2 is not designed with large values for -k in mind, and when aligning reads to long,

repetitive genomes, large -k could make alignment much slower.

--max-seeds <int> HISAT2, like other aligners, uses seed-and-extend approaches. HISAT2 tries to extend seeds to

full-length alignments. In HISAT2, --max-seeds is used to control the maximum number of seeds that

will be extended. For DNA-read alignment (--no-spliced-alignment), HISAT2 extends up to these many seeds

and skips the rest of the seeds. For RNA-read alignment, HISAT2 skips extending seeds and reports

no alignments if the number of seeds is larger than the number specified with the option,

to be compatible with previous versions of HISAT2. Large values for --max-seeds may improve alignment

sensitivity, but HISAT2 is not designed with large values for --max-seeds in mind, and when aligning

reads to long, repetitive genomes, large --max-seeds could make alignment much slower.

The default value is the maximum of 5 and the value that comes with -k times 2.

-a/--all HISAT2 reports all alignments it can find. Using the option is equivalent to using both --max-seeds

and -k with the maximum value that a 64-bit signed integer can represent (9,223,372,036,854,775,807).

--repeat report alignments to repeat sequences directly

Paired-end:

-I/--minins <int> minimum fragment length (0), only valid with --no-spliced-alignment

-X/--maxins <int> maximum fragment length (500), only valid with --no-spliced-alignment

--fr/--rf/--ff -1, -2 mates align fw/rev, rev/fw, fw/fw (--fr)

--no-mixed suppress unpaired alignments for paired reads

--no-discordant suppress discordant alignments for paired reads

Output:

-t/--time print wall-clock time taken by search phases

--un <path> write unpaired reads that didn't align to <path>

--al <path> write unpaired reads that aligned at least once to <path>

--un-conc <path> write pairs that didn't align concordantly to <path>

--al-conc <path> write pairs that aligned concordantly at least once to <path>

(Note: for --un, --al, --un-conc, or --al-conc, add '-gz' to the option name, e.g.

--un-gz <path>, to gzip compress output, or add '-bz2' to bzip2 compress output.)

--summary-file <path> print alignment summary to this file.

--new-summary print alignment summary in a new style, which is more machine-friendly.

--quiet print nothing to stderr except serious errors

--met-file <path> send metrics to file at <path> (off)

--met-stderr send metrics to stderr (off)

--met <int> report internal counters & metrics every <int> secs (1)

--no-head suppress header lines, i.e. lines starting with @

--no-sq suppress @SQ header lines

--rg-id <text> set read group id, reflected in @RG line and RG:Z: opt field

--rg <text> add <text> ("lab:value") to @RG line of SAM header.

Note: @RG line only printed when --rg-id is set.

--omit-sec-seq put '*' in SEQ and QUAL fields for secondary alignments.

Performance:

-o/--offrate <int> override offrate of index; must be >= index's offrate

-p/--threads <int> number of alignment threads to launch (1)

--reorder force SAM output order to match order of input reads

--mm use memory-mapped I/O for index; many 'hisat2's can share

Other:

--qc-filter filter out reads that are bad according to QSEQ filter

--seed <int> seed for random number generator (0)

--non-deterministic seed rand. gen. arbitrarily instead of using read attributes

--remove-chrname remove 'chr' from reference names in alignment

--add-chrname add 'chr' to reference names in alignment

--version print version information and quit

-h/--help print this usage message

(ERR): hisat2-align exited with value 1

> hisat2_extract_snps_haplotypes_UCSC.py

$ hisat2_extract_snps_haplotypes_UCSC.py

usage: hisat2_extract_snps_haplotypes_UCSC.py [-h] [--inter-gap INTER_GAP] [--intra-gap INTRA_GAP] [-v] [--testset] [genome_file] [snp_fname] [base_fname]

Extract SNPs and haplotypes from a SNP file downloaded from UCSC (e.g. http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/snp144.txt.gz)

positional arguments:

genome_file input genome file (e.g. genome.fa)

snp_fname input snp file downloaded from UCSC (plain text or gzipped file is accepted: snp144Common.txt or snp144Common.txt.gz)

base_fname base filename for SNPs and haplotypes

optional arguments:

-h, --help show this help message and exit

--inter-gap INTER_GAP

Maximum distance for variants to be in the same haplotype

--intra-gap INTRA_GAP

Break a haplotype into several haplotypes

-v, --verbose also print some statistics to stderr

--testset print test reads

> hisat2_extract_exons.py

$ hisat2_extract_exons.py

usage: hisat2_extract_exons.py [-h] [-v] [gtf_file]

Extract exons from a GTF file

positional arguments:

gtf_file input GTF file (use "-" for stdin)

optional arguments:

-h, --help show this help message and exit

-v, --verbose also print some statistics to stderr

> hisat2_extract_splice_sites.py

$ hisat2_extract_splice_sites.py

usage: hisat2_extract_splice_sites.py [-h] [-v] [gtf_file]

Extract splice junctions from a GTF file

positional arguments:

gtf_file input GTF file (use "-" for stdin)

optional arguments:

-h, --help show this help message and exit

-v, --verbose also print some statistics to stderr

ラン
1、indexing

まずリファレンスゲノム（遺伝子上にイントロンがある真核生物ゲノム）のindexを作成する。

hisat2-build -p 12 input_genome.fa genome_index

-p number of threads
--ss Splice site file name
--exon Exon file name
--snp SNP file name
--bmax max bucket sz for blockwise suffix-array builder
--large-index force generated index to be 'large', even if ref
has fewer than 4 billion nucleotides

- hisat2-build は、サフィックスが .1.ht2, .2.ht2, .3.ht2, .4.ht2, .5.ht2, .6.ht2, .7.ht2, そして .8.ht2 の 6 つのファイルのセットを出力する。ラージゲノムの大きなインデックスの場合、これらのサフィックスは ht2l で終了する。これらのファイルが一緒になってインデックスを構成する。インデックスが構築されると、元のシーケンス FASTA ファイルは HISAT2 によって使用されない。

- メモリが足りない場合、スレッドを減らすか”--bmax”を調整する（マニュアル参照）。

- リファレンスゲノムが40億文字（32bit）を超えなくても大きなインデックスが望ましい場合、--large-indexを指定して、hisat2-buildが大きなインデックスを生成するように強制することができる。

- 複数の FASTA ファイルにインデックスを作成する場合は、ファイル名をカンマで区切ってすべてのファイルを指定する。

その他、indexに関する補足

--snpで指定する既知SNPsのファイルはdnSNPなどからhisat2_extract_snps_haplotypes_UCSC.pyコマンドで作成することができる。
--ssと--exonの両方のフラグを立てることで、既知スプライスサイト情報をindex内に含め、ヒトなどのモデル生物のマッピング精度を上げることができる。この--exonで指定するテキストはGTFファイルからhisat2_extract_exons.pyコマンドで作成することができる。
--ssで指定するsplice siteのテキストはhisat2_extract_splice_sites.pyコマンドでGTFファイルから作成することができる。このコマンドを使わずにindexファイルを作成している場合、マッピング時に--known-splicesite-infile splicesites.txt をつけることで既知スプライスサイト情報を考慮したマッピングが可能だが、マニュアルにはindex作成時に含めることが望ましいとの記載がマニュアルにある。
splicesites.txt はGTFファイルからhisat2_extract_splice_sites.py hg38.gtf > splicesites.txtのようにして作成することができる。

２，Mapping

ペアエンドのマッピング（シングルエンドは"-U input.fastq"）

hisat2 -x genome_index -1 pair1.fastq -2 pair2.fastq -k 3 -p 12 -S output.sam

-1　Comma-separated list of files containing mate 1s (filename usually includes _1), e.g. -1 flyA_1.fq,flyB_1.fq. Sequences specified with this option must correspond file-for-file and read-for-read with those specified in <m2>. Reads may be a mix of different lengths. If - is specified, hisat2 will read the mate 1s from the "standard in" or "stdin" filehandle.
-2　Comma-separated list of files containing mate 2s (filename usually includes _2), e.g. -2 flyA_2.fq,flyB_2.fq. Sequences specified with this option must correspond file-for-file and read-for-read with those specified in <m1>. Reads may be a mix of different lengths. If - is specified, hisat2 will read the mate 2s from the "standard in" or "stdin" filehandle.
-k 　it searches for at most <int> distinct, primary alignments for each read. Primary alignments mean alignments whose alignment score is equal or higher than any other alignments. The search terminates when it can't find more distinct valid alignments, or when it finds <int>, whichever happens first. The alignment score for a paired-end alignment equals the sum of the alignment scores of the individual mates. Each reported read or pair alignment beyond the first has the SAM 'secondary' bit (which equals 256) set in its FLAGS field. For reads that have more than <int> distinct, valid alignments, hisat2 does not guarantee that the <int> alignments reported are the best possible in terms of alignment score. Default: 5 (HFM) or 10 (HGFM)
--phred33 　Input qualities are ASCII chars equal to the Phred quality plus 33. This is also called the "Phred+33" encoding, which is used by the very latest Illumina pipelines.
--phred64　 Input qualities are ASCII chars equal to the Phred quality plus 64. This is also called the "Phred+64" encoding.
-p　Launch NTHREADS parallel search threads (default: 1). Threads will run on separate processors/cores and synchronize when parsing reads and outputting alignments. Searching for alignments is highly parallel, and speedup is close to linear. Increasing -p increases HISAT2's memory footprint. E.g. when aligning to a human genome index, increasing -p from 1 to 8 increases the memory footprint by a few hundred megabytes. This option is only available if bowtie is linked with the pthreads library (i.e. if BOWTIE_PTHREADS=0 is not specified at build time).
-S　File to write SAM alignments to. By default, alignments are written to the "standard out" or "stdout" filehandle (i.e. the console).

豊富なパラメータがあり、インサートサイズの最大値など様々設定可能です。例えば、--max-intronlenではインサートサイズの最大値がデフォルトでは500000となっていますが、これはおそらくヒトゲノム向けの設定です。長いイントロンがない生物では下げたほうが良いと思います。いくつかの論文では3000から5000程度の値が使われています（数値に根拠があるのか不明）。

詳細はマニュアルから確認してください。

３、sam形式 => bam形式変換

samtoolsを使います。samtools viewでbamに変換後にcoordinateソートしても良いですし、sortコマンドでsamからcoordinateソートされたbamに１コマンドで変換する事も出来ます。

samtools sort -@ 6 -O BAM hisat2.sam > output.bam
#indexxing
samtools index -@ 6 output.bam

samをexportせず、そのままcoordinate sortされたbamとbam.baiを作成するならhisat2の結果をsamtoolsにpipingします。スレッドは12指定。samtools sortはI/O集約型なのでスレッド数は控えめにします。

hisat2 -x genome_index -1 pair1.fastq -2 pair2.fastq -k 3 -p 12 \
 |samtools sort -@ 6 -O BAM - > sorted.bam && \
 samtools index -@ 6 sorted.bam

その他

１、中川先生がRNA seq解析のフローを説明した記事を書かれています。AsperaによるSRAのシーケンスデータの高速ダウンロードから、マッピング、定量、検定、Rを使った分析まで丁寧に説明されています。

https://ncrna.jp/blog/item/385-2018-aspera-connect-hisat2

リンクありがとうございます。

２、multiqcを使って統合レポートを出力するにはhisat2実行時に --summary-file を使ってレポートファイルを出力しておく必要がある。

hisat2 --summary-file rep1_hisat2 --new-summary -x genome_index \
 -1 sample1_QT1.fq.gz -2 sample1_QT2.fq.gz -k 3 -p 20 \
 |samtools sort -@ 20 -O BAM - > cont_rep1.bam && \
samtools index -@ 12 cont_rep1.bam

hisat2 --summary-file rep2_hisat2 --new-summary -x genome_index \
 -1 sample2_QT1.fq.gz -2 sample2_QT2.fq.gz -k 3 -p 20 \
 |samtools sort -@ 20 -O BAM - > cont_rep2.bam && \
samtools index -@ 12 cont_rep2.bam

hisat2 --summary-file rep3_hisat2 --new-summary -x genome_index \
 -1 sample3_QT1.fq.gz -2 sample3_QT2.fq.gz -k 3 -p 20 \
 |samtools sort -@ 20 -O BAM - > cont_rep3.bam && \
samtools index -@ 12 cont_rep3.bam

hisat2 --summary-file tr1_hisat2 --new-summary -x genome_index \
 -1 sample4_QT1.fq.gz -2 sample4_QT2.fq.gz -k 3 -p 20 \
 |samtools sort -@ 20 -O BAM - > tre_rep1.bam && \
samtools index -@ 12 tre_rep1.bam

hisat2 --summary-file tr2_hisat2 --new-summary -x genome_index \
 -1 sample5_QT1.fq.gz -2 sample5_QT2.fq.gz -k 3 -p 20 \
 |samtools sort -@ 20 -O BAM - > tre_rep2.bam && \
samtools index -@ 12 tre_rep2.bam

hisat2 --summary-file tr3_hisat2 --new-summary -x genome_index \
 -1 sample6_QT1.fq.gz -2 sample6_QT2.fq.gz -k 3 -p 20 \
 |samtools sort -@ 20 -O BAM - > tre_rep3.bam && \
samtools index -@ 12 tre_rep3.bam

multiqc .