一?本系列的開啟原因

轉(zhuǎn)錄組分析已經(jīng)成為當(dāng)今生物信息分析領(lǐng)域最最最常見(jiàn)的一種組學(xué)分析了。并且近些年因?yàn)槠鋬r(jià)格的大幅度下降以及分析手段的日漸成熟，已經(jīng)成為打開每個(gè)生物學(xué)問(wèn)題的所必須經(jīng)歷的第一道門檻了。雖然這是一種常見(jiàn)的分析，但是很多人還是對(duì)這些問(wèn)題不夠明晰。正好，本公眾號(hào)最近需要重新分析一批轉(zhuǎn)錄組數(shù)據(jù)，這里就一邊分析一邊總結(jié)。希望這一個(gè)系列能夠給其他讀者展示我已經(jīng)踩過(guò)的一些坑，從而幫助讀者開始完成這一分析，進(jìn)入更深一步的功能驗(yàn)證環(huán)境。

那么，我們要做轉(zhuǎn)錄組分析，不可避免的就是要使用到hisat2、Stringtie和edgeR這三個(gè)軟件，并且這三個(gè)軟件也是目前轉(zhuǎn)錄組分析的主流流程。那么本系列前幾篇將先介紹一下這三個(gè)軟件的使用。考慮到原文檔太過(guò)繁瑣，我會(huì)盡量簡(jiǎn)化一下hisat的使用文檔。但是，如果讀者想要知道細(xì)節(jié)，或者覺(jué)得我介紹的不夠準(zhǔn)確，可以去官網(wǎng)查看。雖然盡可能的確保準(zhǔn)確，但本公眾號(hào)仍對(duì)介紹的準(zhǔn)確性不做任何保證。不過(guò)，如果有錯(cuò)誤，后續(xù)也會(huì)進(jìn)行相應(yīng)的更正。那么，接下來(lái)就讓我們先了解一下在使用hisat2軟件前需要做哪些準(zhǔn)備吧。

一?建立索引

現(xiàn)如今，隨著基因組的大量出現(xiàn)，轉(zhuǎn)錄組分析也更多的偏向于有參轉(zhuǎn)錄組，因此，hisat2首先提供了一個(gè)hisat2-build選項(xiàng)，用于對(duì)所使用的基因組進(jìn)行處理，即建立索引，方便分析。

1 主參數(shù)

<reference_in>:?以逗號(hào)分隔的FASTA文件列表，如每個(gè)染色體的fa文件。如果只有一個(gè)文件，那么不需要任何逗號(hào)。

<ht2_base>:?要寫入的索引文件的前綴名（會(huì)產(chǎn)生許多以此前綴名為開頭的文件）。

2 options選項(xiàng)

-f:?用于比對(duì)的輸入文件是FASTA文件。

-c:?在命令行中給出參考序列，即<reference_in>是一個(gè)用逗號(hào)分隔的序列列表，而不是FASTA文件的列表。

--large-index:?強(qiáng)制hisat2-build建立一個(gè)大的索引，即使參考序列的長(zhǎng)度小于~40億bp核苷酸長(zhǎng)度。

-a/--noauto:?禁用默認(rèn)行為，即hisat2-build根據(jù)可用內(nèi)存自動(dòng)選擇。

--bmax <int>:?一個(gè)區(qū)塊中允許的最大后綴數(shù)。

--bmaxdivn <int>:?一個(gè)區(qū)塊中允許的最大后綴數(shù)，以參考長(zhǎng)度的分?jǐn)?shù)表示。

--dcv <int>:?使用<int>作為差異覆蓋樣本的周期。

--nodc:?禁用差值覆蓋樣本的使用。

-r/--noref:?不建立索引的NAME.3.ht2和NAME.4.ht2部分，而這些部分包含參考序列的位包版本，可以用于成對(duì)端比對(duì)。

-3/--justref:?只建立索引中NAME.3.ht2和NAME.4.ht2的部分，其中包含參考序列的位包版本，用于配對(duì)端比對(duì)。

-o/--offrate <int>:?為了將比對(duì)結(jié)果映射到參考序列上的位置，有必要用基因組上的相應(yīng)位置來(lái)注釋（"標(biāo)記"）部分或全部Burrows-Wheeler行。

-t/--ftabchars<int ftab>:?是用于計(jì)算與查詢的第一個(gè)<int>字符有關(guān)的初始Burrows-Wheeler范圍的查找表。

--localoffrate<int>:?這個(gè)選項(xiàng)控制了在本地索引中標(biāo)記多少行：索引器將每2^<int>行標(biāo)記一次。

--localftabchars<int>:?本地ftab是本地索引中的查找表。

-p <int>:?構(gòu)建并行線程（默認(rèn)：1）。

--snp <path>:?提供SNPs的列表格式如下。

--haplotype <path>:?提供一個(gè)單倍型列表的文件，格式如下。

--ss <path>:?注意這個(gè)選項(xiàng)應(yīng)該和下面的--exon選項(xiàng)一起使用，提供一個(gè)剪接點(diǎn)的列表。

--exon <path>:?注意這個(gè)選項(xiàng)應(yīng)該和上面的--ss選項(xiàng)一起使用。提供一個(gè)外顯子列表，可以使用hisat2_extract_exons.py（HISAT2軟件包），從GTF文件中提取外顯子。

--seed <int>:?使用<int>作為偽隨機(jī)數(shù)發(fā)生器的種子。

--cutoff <int>:?只索引參考序列的前<int>個(gè)堿基（跨序列累積），而忽略其余部分。

-q/--quiet:?hisat2-build在默認(rèn)情況下是靜默的。使用這個(gè)選項(xiàng)，hisat2-build將只打印錯(cuò)誤信息。

-h/--help:?打印使用信息并退出。

--version:?打印版本信息并退出。

3 示例

上面的命令表明的含義有以下幾點(diǎn)：線程為2，基因組序列為genome.fa，生成索引文件的前綴為genome。

二?生成可變剪切的文件信息

注意：最好使用使用注釋的轉(zhuǎn)錄本建立的索引（如genome_tran或genome_snp_tran），這比使用這個(gè)選項(xiàng)效果更好。提供已經(jīng)包含在索引中的剪接位點(diǎn)沒(méi)有影響。

hisat2_extract_splice_sites.py:?包含在HISAT2軟件包中，

genes.gtf:?基因注釋文件

spliceites.txt:?提供給HISAT2的剪接點(diǎn)列表。

三?比對(duì)gap

本部分主要介紹一下比對(duì)過(guò)程中該如何根據(jù)測(cè)序長(zhǎng)度和讀長(zhǎng)來(lái)定義雙端測(cè)序中兩個(gè)read中所應(yīng)該產(chǎn)生的gap，方便在下一章中介紹如何設(shè)定參數(shù)進(jìn)行比對(duì)。

gap計(jì)算方法：如果指定了-I 60，并且一個(gè)配對(duì)端排列由兩個(gè)20bp的排列組成，在適當(dāng)?shù)姆较蛏?，它們之間有20bp的間隙，那么這個(gè)排列被認(rèn)為是有效的（只要-X也被滿足）。在這種情況下，19bp的間隙是無(wú)效的。如果同時(shí)使用了修剪選項(xiàng)-3或-5，-I約束將被應(yīng)用于未修剪的配體。

gap大小與比對(duì)性能：-I和-X之間的差異越大，HISAT2的運(yùn)行速度就越慢。這是因?yàn)?I和-X之間的差異越大，HISAT2就需要掃描一個(gè)更大的窗口來(lái)確定是否存在一致的排列。對(duì)于典型的片段長(zhǎng)度范圍（200到400個(gè)核苷酸），HISAT2是非常有效的。

六 慣例小結(jié)

通過(guò)構(gòu)建索引、生成可變剪切文件，我們就可以獲得除測(cè)序文件外，所有應(yīng)該被準(zhǔn)備好的數(shù)據(jù)了。那么，接下來(lái)，就輪到正式比對(duì)啦。

本公眾號(hào)開發(fā)的相關(guān)軟件，Multi-omics Hammer軟件和Multi-omics Visual軟件歡迎大家使用。文末是本公眾號(hào)在其他平臺(tái)的賬戶，也歡迎大家關(guān)注并多提意見(jiàn)。

簡(jiǎn)書：WJ的生信小院

公眾號(hào)：生信小院

博客園：生信小院

最后，也歡迎各位大佬能夠在本平臺(tái)上：1傳播和講解自己發(fā)表的論文；2：發(fā)表對(duì)某一科研領(lǐng)域的看法；3：想要達(dá)成的合作或者相應(yīng)的招聘信息；4：展示自己以尋找博后工作或者博士就讀的機(jī)會(huì);5:博導(dǎo)提供博后工作或者博士攻讀機(jī)會(huì)，都可以后臺(tái)給筆者留言。希望本平臺(tái)在進(jìn)行生信知識(shí)分享的同時(shí)，能夠成為生信分析者的交流平臺(tái)，能夠?qū)崿F(xiàn)相應(yīng)的利益互補(bǔ)和雙贏（不一定能實(shí)現(xiàn)，但是夢(mèng)想總得是有的吧）。

另外，怎么說(shuō)呢，投幣也可，不強(qiáng)求，但奢求。