如何解析對比基因組工具hisat2

這篇文章將為大家詳細講解有關如何解析對比基因組工具hisat2，文章內容質量較高，因此小編分享給大家做個參考，希望大家閱讀完這篇文章后對相關知識有一定的了解。

由于測序儀機器讀長的限制，在構建文庫的過程中首先需要將DNA片段化，測序得到的序列只是基因組上的部分序列。為了確定測序reads在基因組上的位置，需要將reads比對回參考基因組上，這個步驟叫做mapping。

在進行mapping時，需要考慮以下幾個因素

1. 硬件資源的消耗

通常來說，基因組越大，占用的內存越大。對于大型基因組，比如人類基因組而言，優化內存消耗是很關鍵的一點。

2. 運行速度

隨著測序價格的下降和數據深入挖掘的需求，測序量越來越大，海量測序reads的比對，要求速度上必須夠快。

3. 準確性

SNP/indel, 測序錯誤率等因素都使得測序的reads和基因組上的原始序列會存在幾個bp的誤差，所以mapping的算法必須支持堿基的錯配，或者是gap的存在。同時由于測序的短序列可能和基因組多個位置存在同源，一條reads會比對到基因組上多個位置。雙端測序技術在一定程度上能夠校正多個位置，因為雙端reads 來自同一個DNA片段，二者在基因組上的位置不會相距太遠，但是僅靠這一點并不能解決所有的同源比對，這就要求比對算法對多個位置進行判斷和打分，給出比對結果的可靠性。

4. RNA

對于轉錄組數據， 真核生物可變剪切的存在，導致cDNA片段在基因組上的位置并不是連續的，中間可能存在內含子。在比對轉錄組數據時，就需要考慮跳過剪切位點。

目前mapping的工具有很多，比如bwa, hisat, star等。hisat 是其中速度最快的，是tophat軟件的升級版本。采用了改進的FM index 算法，對于人類基因組，只需要4.3GB左右的內存。同時支持DNA和RNA數據的比對，軟件官網如下

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

目前最新版為為hisat2. 安裝過程如下

wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip
unzip hisat2-2.1.0-Linux_x86_64.zip

下載解壓縮即可。

在進行比對前，首先需要對參考基因組建立索引, 基本用法如下

hisat2-build -p 20   hg19.fa hg19

對于轉錄組數據，在構建索引時，可以通過gtf文件，得到剪切位點和exon的信息，用法如下

hisat2_extract_splice_sites.py hg19.gtf > hg19.ss
hisat2_extract_exons.py hg19.gtf > hg19.exon
hisat2-build -p 20  --ss hg19.ss --exon hg19.exon  hg19.fa hg19

hisat2 支持多種格式的輸入文件，常見格式有以下兩種

1. fasta

2. fastq

   
   

-f參數表示輸入問下格式為fasta, -q參數表示輸入文件格式為fastq。輸入文件可以是經過gzip壓縮之后的文件，默認輸入文件是fastq格式。

對于單端數據，采用-U指定輸入文件；對于雙端數據，采用-1和-2分別指定R1端和R2端的輸入文件。

reads比對到基因組上的一個位置，我們稱之為一個alignment。 軟件會對所有的alignments 進行打分和判斷，能夠符合過濾條件的alignment 稱之為valid  alignment， 只有valid alignments , 才會輸出。

和blast類似，每個alignment也有對應的打分機制。hisat 從以下幾個方面對alignment 進行打分

1. 錯配堿基罰分

錯配堿基的罰分通過--mp參數指定，其值為逗號分隔的兩個數字，第一個數字為最大的罰分，第二個數字為最小的罰分

2. reads上的gap罰分

gap的罰分通過分成兩個部分，第一次出現gap的罰分和gap延伸的罰分，reads上的gap罰分通過--rdg參數指定，其值為逗號分隔的兩個數字，第一個數字為gap第一個位置的罰分，第二個數字為gap延伸的罰分。

3. reference上的gap罰分

reference上的gap罰分通過--rdg參數指定，其值為逗號分隔的兩個數字，第一個數字為gap第一個位置的罰分，第二個數字為gap延伸的罰分。

經過一系列的罰分機制，每個alignment會有一個對應的得分，然后會根據一個閾值，來判斷這個得分是否滿足valid  alignment的要求。

hisat通過--score--min參數指定該閾值，指定方式是一個和reads程度相關的函數，默認值為L,0,-0.2, 對應函數為

f(x) = 0 - 0.2 * x

根據reads長度，可以計算出得分的閾值，大于該閾值的alignment 被認為是valid alignment , 才可能被輸出。L代表線性函數，此外，也支持其他類型的函數，比如常量，自然對數等，更多選擇請參考官方文檔。

一條reads可能會擁有多個valid  alignments, 在輸出時，并不會輸出所有的alignments, 而是只輸出-k參數指定的N個alignments，-k參數的默認值為5。
輸出結果以SAM格式保存，默認輸出到屏幕上，可以通過-S參數指定輸出文件。

通常情況下，默認參數就能夠滿足我們的需求了。單端數據比對的用法如下

hisat -x hg19 -p 20 -U reads.fq -S align.sam

雙端數據用法如下

hisat -x hg19 -p 20 -1 R1.fq -2 R2.fq -S align.sam

關于如何解析對比基因組工具hisat2就分享到這里了，希望以上內容可以對大家有一定的幫助，可以學到更多知識。如果覺得文章不錯，可以把它分享出去讓更多的人看到。