tags： HISAT2 RNA-seq
HISAT2 發(fā)表的文章中強調了它的速度很快晤愧，我就測試了一下這個工具。
HISAT2 建立索引：
hisat2-build -p 4 rRNA.fa rRNA.fa.tran
然而沒多久就看到這樣的警告：
Reading reference sizes
Warning: Encountered reference sequence with only gaps
Warning: Encountered reference sequence with only gaps
Warning: Encountered reference sequence with only gaps
Warning: Encountered reference sequence with only gaps
Warning: Encountered reference sequence with only gaps
只是警告，并沒有報錯。
HISAT2 建參考索引很慢，等 HISAT2 建完索引（建索引花了 16 個小時），然后用 HISAT2 比對 RNA-seq 數據測試嫂拴。
hisat2 -S /dev/null -p 4 -x rRNA.fa.tran -1 A_1.fq.gz -2 A_2.fq.gz --un-conc-gz A_filter_rRNA_%.fq.gz
幾分鐘之后報錯了：
(ERR): hisat2-align died with signal 11 (SEGV) (core dumped)
github 上有人在 HISAT2 項目中報告過這個錯誤，雖然沒有最終討論出解決辦法贮喧，但是都覺得跟建索引不完整有關筒狠，或許與建索引時候的警告有關。我查了一些資料箱沦，綜合 biostar 和 SEQanswer 中的討論辩恼，建立索引時遇到的警告是由于參考序列中存在大段的 n 堿基導致的，例如其中一條 fasta 中 n （我遇到的是小寫的 n）太多谓形。解決辦法也很簡單灶伊，過濾掉參考序列中長度小于 50bp 的 contig 和序列中連續(xù) n 堿基超過 40bp 的contig 。然后重新建索引寒跳，就沒有任何警告了聘萨，但是會損失部分參考序列。這樣處理之后建的索引再用 HISAT2 比對 RNA-seq 數據童太，就沒有問題了米辐。
#!/usr/bin/env python3
from Bio import SeqIO
long_seq = []
for record in SeqIO.parse("rRNA.fa","fasta"):
if len(record.seq) > 50 and not 'n'*40 in record.seq:
long_seq.append(record)
SeqIO.write(long_seq, "filtered.rRNA.fa", "fasta")
HISAT2 比對速度確實很快碾牌，一個樣本轉錄組數據比對 2G 的核糖體 RNA 參考基因組約 25 分鐘，bowtie2 需要 170 分鐘（線程數都是 4）儡循。 bowtie2 的 local 模式比對出 21% 的 rRNA 污染，而 HISAT2 比對 2% 的 rRNA 污染征冷，差異也挺大的……
但是择膝，HISAT2 的線程數不能設置太高，用戶手冊中建議對人類參考基因組進行比對時检激，線程數設置在 1-8 之間肴捉。我用文獻 Transcript-level expression analysis of RNA-seq
experiments with HISAT, StringTie and Ballgown 中的數據測試也發(fā)現當線程數超過 12 時整個比對步驟消耗的時間反而會增加。
