macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

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

2019 6/18 コマンド追記

2019 6/26 インストール追記

2019 6/28 samtoolsコマンドエラー修正

 

 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を紹介する。

 

公式ページ

f:id:kazumaxneo:20180714114102p:plain

マニュアル

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

すでに色々な人が日本語でも記事にしてますが、このブログでも依存関係で取り上げる機会が多いのでまとめておきます。

 

インストール

brewやcondaで導入できる。

#bioconda(linux)
conda install -c bioconda -y hisat2

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

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

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

> hisat2 -h

r$ hisat2 -h

HISAT2 version 2.0.5 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

 

 Alignment:

  --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)

 

 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

 

 Scoring:

  --ma <int>         match bonus (0 for --end-to-end, 2 for --local) 

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

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

  --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> (default: 5) report up to <int> alns per read

 

 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.)

  --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          supppress header lines, i.e. lines starting with @

  --no-sq            supppress @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

 

ラン

まずリファレンスのindexを作成する。

hisat2-build input_genome.fa genome_index

ペアエンドのマッピングを実行する。(シングルエンドfastqは"-U input.fastq")

hisat2 -x genome_index -1 pair1.fastq -2 pair2.fastq -k 3 -p 24 -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).

豊富なパラメータがあり、インサートサイズの最大値など様々設定可能です。詳細はマニュアルから確認してください。

 

samをexportせず、そのままcoordinate sortされたbamとbam.baiを作成。スレッドは12指定。

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

 

追記

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

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

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

 

 

引用

HISAT: a fast spliced aligner with low memory requirements.

Kim D, Langmead B, Salzberg SL

Nat Methods. 2015 Apr;12(4):357-60.

 

参考ページ

二階堂さんのマニュアル

RNA-seqのリードをHISAT2でゲノムへ高速にマッピングする | Tips for NGS Data Analysis

アメリエフさんのブログ

http://staffblog.amelieff.jp/entry/2016/07/01/144419

ノート

http://mctakashima.blog.jp/archives/2587379.html

 

HISAT2 for aligning reads against human population

http://nextgenseek.com/2015/11/hisat2-for-aligning-reads-against-human-population/