生信第二步:建立索引及比對

2021-02-15 小西瓜與大生信

xio'a做完質控與過濾之後，轉錄組分析上遊最核心部分就來了：比對。

那我們得到了測序文件，要幹嘛？為什麼要對比？如何比對有哪些軟體的選擇？別說你們我自己想想都頭大。

第一個問題:得到了乾淨測序文件，要幹嘛？


對於NGS流程來說肯定就是為了找到SNP和Indel，那Gwas流程就是找差異表達基因或尋找新的可變剪切。總之任何一個RNA-seq流程就一定會用到比對。這取決於你做什麼和下一步分析的必要性。

第二個問題:為什麼要比對？


在比對之前，我們得了解比對的目的是什麼？RNA-Seq數據比對和DNA-Seq數據比對有什麼差異？
RNA-Seq數據分析分為很多種，比如說找差異表達基因或尋找新的可變剪切。如果找差異表達基因單純只需要確定不同的read計數就行的話，我們可以用bowtie, bwa這類比對工具，或者是salmon這類align-free工具，並且後者的速度更快。
但是如果你需要找到新的isoform，或者RNA的可變剪切，看看外顯子使用差異的話，你就需要TopHat, HISAT2或者是STAR這類工具用於找到剪切位點。因為RNA-Seq不同於DNA-Seq，DNA在轉錄成mRNA的時候會把內含子部分去掉。所以mRNA反轉的cDNA如果比對不到參考序列，會被分開，重新比對一次，判斷中間是否有內含子。
連結：https://www.jianshu.com/p/681e02e7f9af

第三個問題：如何比對,有哪些軟體的選擇？

在2016年的一篇綜述A survey of best practices for RNA-seq data analysis，提到目前有三種RNA數據分析的策略。那個時候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已經被它的作者推薦改用HISAT進行替代。


最近的Nature Communication發表了一篇題為的Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis的文章--被稱之為史上最全RNA-Seq數據分析流程，也是我一直以來想做的事情，只不過他們做的超乎我的想像。文章中在基於參考基因組的轉錄本分析中所用的工具，是TopHat,HISAT2和STAR，結論就是HISAT2找到junction正確率最高，但是在總數上卻比TopHat和STAR少。從這裡可以看出HISAT2的二類錯誤(納偽）比較少，但是一類錯誤（棄真）就高起來。
就唯一比對而言，STAR是三者最佳的，主要是因為它不會像TopHat和HISAT2一樣在PE比對不上的情況還強行把SE也比對到基因組上。而且在處理較長的read和較短read的不同情況，STAR的穩定性也是最佳的。
就速度而言，HISAT2比STAR和TopHat2平均快上2.5~100倍。

如果學習RNA-Seq數據分析，上面提到的兩篇文獻是必須要看上3遍以上的，而且建議每隔一段時間回顧一下。但是如果就比對工具而言，基本上就是HISAT2和STAR選一個就行。

參考連結：https://www.jianshu.com/p/681e02e7f9af

本文我將著重降解HISAT2和bowtie、bowtie2的使用（主要就是我現在做的轉錄組和samllRNA會使用用到這兩款軟體相對來說更熟悉一些。）這麼多主流軟體各有優點，先給自己挖個坑，等過幾天時間充足了分別講一下各個軟體的優缺點，及該如何根據自己的測序文件選擇最佳比對軟體。

比對第一步：建立index


問題1：什麼是索引？

答：舉個小例子,在圖書館書籍有成千上萬本書，其中各種散文、傳記、文學、科普、小說、詩集等等，要是雜亂無章的去找那就真是太痛苦了。要是圖書管理員給你按照類別分好類。包括指定在那個書架，是不是就會方便很多？ 同理我們做比對的時候有成千上萬可能還上億，又來索引文件index是不是就更輕鬆？

問題2：如何獲得和建立index？

答：大多數模式物種都會有自己建立的index，一些非模式物種可能會有，但是一般來說都沒有，沒有現成的index，我們就需要自己重新構建索引；包括外顯子、剪切位點及SNP索引的建立。

因為我的課題是葫蘆科相關作物的育種和遺傳鑑定，典型的植物非模式物種，那我這裡就開始講一講自己該如何讓去構建index。

1.1用hisat2構建自己的index

# hisat2的安裝

# 源文件安裝

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

unzip hisat2-2.1.0-beta-Linux_x86_64.zip

export PATH=/home/xxx/xxx/bin/hisat2-2.1.0-beta:$PATH

source ~/.bashrc

# conda 安裝

conda install hisat2

建立索引

# 其實hisat2-buld在運行的時候也會自己尋找exons和splice_sites，但是先做的目的是為了提高運行效率

extract_exons.py watermelen.gtf > genome.exons

extract_splice_sites.py watermelen.gtf > genome.ss

#同時支持將snp文件加入

extract_snps.py watermelen.txt > genome.snp

# 建立index， 必須選項是基因組所在文件路徑和輸出的前綴
hisat2-build -p 16 --snp genomne.snp --ss genome.ss --exon genome.exons watermelen.fa xigua

參數簡單的說明

extract_exons.py watermelen.gtf > genome.exons&

將gtf文件中的exons（外顯子）

extract_splice_sites.py watermelen.gtf > genome.ss

將gtf文件中的剪切位點提取出來

extract_snps.py watermelen.txt > genome.snp

將SNP位置信息提取出來

hisat2-build -p 16 --snp genomne.snp --ss genome.ss --exon genome.exons watermelen.fa xigua

-p 指定線程數  --snp將snp信息加入 --ss將剪接位點信息加入 --exon將外顯子信息加入 watermlen.fa是參考基因組  xigua是index文件名

ps：說了半天第一步才搞完真的痛苦，今天同窗問我一個問題，我覺得很有意思可以拿出來簡單說一下。 他按照我給的步驟走了一遍說沒辦法提取exons，提取出來為0，當時我第一反應是不是參數錯了，從到到尾看了一遍，沒錯，很糾結，我記得我完整的做過這個物種的參考基因組和gtf文件，難受，懷疑我之前都做錯了。然後他告訴我說看了gtf裡面沒有exon文件，我開始還不信，後面自己看了一眼。果然沒有，在它較新的一個版本中有就有。後面百度了一些東西，發現是早期版本沒有這麼完善也就是沒有完整注釋，這種情況很正常。其實我們這樣寫入外顯子、剪切位點、snp是為了提高運行效率，hisat2自己就可以找到。所以你的參考基因組文件是早期版本也不用擔心哦。可以簡化參數為：

hisat2-build -p 16  watermelen.fa watermelen

$ hisat2-build -p 16  watermelen.fa xigua

比對第二步：將測序文件比對到參考基因組


老規矩hisat2 -h 查看用法說明

(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/data$ hisat2 -h
HISAT2 version 2.1.0 by Daehwan Kim (infphilo@gmail.com, www.ccb.jhu.edu/people/infphilo)
Usage:
hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-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).
<SRA accession number> Comma-separated list of SRA accession numbers, e.g. --sra-acc SRR353653,SRR353654.
<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'.

基本參數說明

-x <hisat2-idx>
參考基因組索引文件的前綴。

-1 <m1>
雙端測序結果的第一個文件。若有多組數據，使用逗號將文件分隔。Reads的長度可以不一致。

-2 <m2>
雙端測序結果的第二個文件。若有多組數據，使用逗號將文件分隔，並且文件順序要和-1參數對應。Reads的長度可以不一致。

-U <r>
單端數據文件。若有多組數據，使用逗號將文件分隔。可以和-1、-2參數同時使用。Reads的長度可以不一致。

–sra-acc <SRA accession number>
輸入SRA登錄號，比如SRR353653，SRR353654。多組數據之間使用逗號分隔。HISAT將自動下載並識別數據類型，進行比對。

-S <hit>
指定輸出的SAM文件。

參考連結：https://www.jianshu.com/p/479c7b576e6f

# hisat2的安裝

hisat2 -t -x index/xigua -1 watermelon/P1.R1.clean.fastq -2 watermelon_/P1.R2.clean.fastq -S sam/P1.sam &

參數說明（-t指定線程數， -x指定索引位置和參考序列，-1 指定雙端測序左鏈的位置，-2指定雙端測序右鏈的位置，-S輸出文件位置及輸出文件名稱）

得到了輸出結果類似下圖（圖是網上找的，因為我自己當時忘記保存圖片，重新開始又太耗費時間，理解一下）

輸出結果可以看出總體比對率，正向比對到了多少，反向比對到多少，沒有比對多少的有多少。。。等等重要的信息一目了然。

這個時候你就應該問了：那我跟你一樣我忘記保存了呢？我看別人統計都有比對率多少那我難道這個錯過了我就不知道了？（我不管，你就是有這個問題）

  這個時候就需要用到另一個神器了 那就是samtools軟體了。

  因為我得到的是sam文件，你自己做完之後發現動不動就20.30個G，無論是伺服器還是個人pc那都是及其恐怖的一個數據。 如果我們想查看數據，我怕等你下載完就過去了10天半個月。那有沒有吧這個sam文件給簡化呢？嘿嘿嘿，那就是將sam文件轉換成bam文件。（我知道你想問：你不是說看比對率嗎？怎麼就講啥sam、bam文件了，不要著急。我們慢慢來嘛）

先看看SAMTOOLS這個軟體能幹啥吧

  SAMtools的wiki介紹：SAMtools Wiki
  官方手冊：Manual page from samtools-1.9
  必看：詳細了解SAMtools的用法和來龍去脈

  而目前處理SAM格式的工具主要是SAMTools，這是Heng Li大神寫的.除了C語言版本，還有Java的Picard，Python的Pysam，Common lisp的cl-sam等其他版本。SAMTools的主要功能如下：

最常用的三板斧就是格式轉換，排序，索引。而進階教程就是看文檔提高。

參考連結：https://www.jianshu.com/p/681e02e7f9af
  我們主要講將sam文件轉換成bam文件，還有對其排序、建立索引和查看比對率。

基礎用法：

##第三步將sam文件轉換成bam文件並對bam文件進行排序。

samtools  view -bS sam/P1.sam > bam/P1.bam &(轉換成bam文件)

samtools view -bS sam/P2.sam >bam/P2.bam &(轉換成bam文件)

samtools sort sam/P1.bam -o bam/P1.sort.bam & （對P1按照染色體位置排序，-n就是按照read name排序，-r 按照染色體排序，不輸入則默認染色體排序）

samtools sort sam/P2.bam -o bam/P2.sort.bam & （對P2按照染色體位置排序，-n就是按照read name排序，-r 按照染色體排序，不輸入則默認染色體排序）

##第四步對參考基因組進行索引和排序後的bam文件進行索引（方便在IGV中查看，建立索引的目的就是為了讓文件快速找到，同時IGV中查看也必須有索引文件 不然無法打開）

samtools faidx  watermlen.fa（對參考基因組建立索引）

samtools index P1.sort.bam

samtools index P2.sort.bam(對bam文件建立索引)

  你發現bam文件比sam文件數據小了至少一半，sort之後bam文件又變小了

bam文件和sam文件之間的轉換可以多去了解我在上面放出了連結。簡單來說

為什麼 BAM 文件 sort 之後體積會變小因為 BAM 文件是壓縮的二進位文件，對文件內容排序之後相似的內容排在一起，使得文件壓縮比提高了，因此排序之後的 BAM 文件變小了，相對應的 SAM 文件就是純文本文件，對SAM 文件進行排序就不會改變文件大小。而且由於 RNA-seq 中由於基因表達量的關係，RNA-seq 的數據比對結果 BAM 文件使用 samtools 進行 sort 之後文件壓縮比例變化會比DNA-seq 更甚。另外，samtools 對 BAM 文件進行排序之後那些沒有比對上的 reads 會被放在文件的末尾。

回到最開始的問題我怎麼去查看比對率呢？

#查看比對情況

samtools flagstat XXX.bam  

比對結果

(rna-seq) [lmyang@compute01 bam3]$ samtools flagstat AB1.sort.bam 

64149401 + 0 in total (QC-passed reads + QC-failed reads)

3161137 + 0 secondary

0 + 0 supplementary

0 + 0 duplicates

60566539 + 0 mapped (94.41% : N/A)

60988264 + 0 paired in sequencing

30494132 + 0 read1

30494132 + 0 read2

56272878 + 0 properly paired (92.27% : N/A)

56357432 + 0 with itself and mate mapped

1047970 + 0 singletons (1.72% : N/A)

33272 + 0 with mate mapped to a different chr

26411 + 0 with mate mapped to a different chr (mapQ>=5)

  在這裡我提一個小問題：你們看看samtools得到的比對率和直接hisat2得到的對率有什麼不同？小驚喜等著你喲。（又是一個坑，可以自己網上查查有什麼解釋，人生處處是驚喜。）

  得到bam文件和它的索引文件之後我們就可以去IGV裡面看看情況如何，比如有沒有snp、indel、可變剪切之類的。沒錯下期預告來了，就是講解IGV的使用。

  講完這個我的狀態已經虛脫了，感覺還是有很多不足，希望大家批評指正！

聲明：本文章內容多參考《簡書》、《微信公眾號》等原創內容，個人對其優質內容進行整理總結，僅用於學習和交流，本人不負責其法律責任如侵權，聯繫本人刪除。

相關焦點

  • STAR:轉錄組數據比對工具簡介
    構建基因組索引運行比對前，首先需要對基因組建立索引，建立索引對應的runMode為genomeGenerate, 基本用法如下STAR --runMode genomeGenerate --runThreadN  20 --genomeFastaFiles hg19.fasta --genomeDir hg19
  • nanopore測序技術專題(二十):tablet可視化比對結果
    而tablet就是一款可視化高通量數據的一款軟體，我們可以通過該軟體直接看到每個位點的比對細節，例如該位點被測序了多少次，每個鹼基具體是什麼，記得還在深圳工作的時候，我們經常使用這款軟體一點點檢測利用二代測序做出的細菌基因組完成圖連接的是否正確，一看就是1個多小時，導致晚上做夢滿腦子都是花花綠綠的鹼基，這算是工傷。
  • 一文看懂聚集索引和非聚集索引的區別
    第二：為什麼聚集索引可以創建在任何一列上，如果此表沒有主鍵約束，即有可能存在重複行數據呢？粗一看，這還真是和聚集索引的約束相背，但實際情況真可以創建聚集索引。分析其原因是：如果未使用 UNIQUE 屬性創建聚集索引，資料庫引擎將向表自動添加一個四字節 uniqueifier 列。
  • 序列比對在biopython中的處理
    歡迎關注」生信修煉手冊」! 序列比對是生物信息學分析中的常見任務，包含局部比對和全局比對兩大算法，局部比對最經典的代表是blast, 全局比對則用於多序列比對。在biopython中，支持對序列比對的結果進行讀寫，解析，以及運行序列比對的程序。
  • 如何進行基因組序列比對?
    FM-index是基於Burrows-Wheeler transform進行全文壓縮和構建索引的算法。備註：大家看到下面的文件第二列是字符串，與下面的幾個截圖不一致，是因為在samtools view時添加了-X參數，可以將第二列的FLAG轉為字符形式，FLAG具體會在比對結果部分解說，在這主要是看FLAG這列最後的數字，「1」表示Paired end測序中的read1,「2」表示Paired end測序中的read2。
  • 手把手教你用DNAMAN進行多序列比對
    DNAMAN做多序列比對所用到的序列並不需要FASTA格式。但每一條序列都需要放在一個單獨的文件中。準備好序列之後，便可進行多序列比對了。打開DNAMAN軟體。點擊序列-比對-多序列比對。將要比對的序列都載入之後，點擊下一步-下一步-下一步-完成。稍等片刻便得出多序列比對結果了。
  • 河南省建立全省聯網離退休人員指紋比對資料庫
    原標題：河南省建立全省聯網離退休人員指紋比對資料庫 　　人民網鄭州2月26日專電 為確保養老保險基金安全，方便離退休人員領取，河南省建立了覆蓋全省企業離退休人員的指紋比對網絡管理系統。　　為此，從2011年6月開始，河南以省養老保險局資料庫為中心，以各省轄市、縣(市、區)、各大行業單位養老保險經辦機構為前臺，廣泛設立指紋採集比對服務網點，著手建立覆蓋全省企業離退休人員的指紋比對網絡管理系統。
  • 多序列比對分析-Dnaman很好用!
    序列載入打開Dnaman軟體，如下圖，第一欄為主菜單欄，有12個常用主菜單；第二欄為工具欄；再下面為瀏覽器欄多序列比對參數設置如下，可點擊channel選擇需要進行比對的序列，注意區分下方DNA/Protein選項，點擊OK進行下一步，其他參數見下圖，如無特殊需求可默認。
  • MySQL InnoDB 索引原理
    聚簇索引和二級索引3.1 聚簇索引每個InnoDB的表都擁有一個索引，稱之為聚簇索引，此索引中存儲著行記錄，一般來說，聚簇索引是根據主鍵生成的。為了能夠獲得高性能的查詢、插入和其他資料庫操作，理解InnoDB聚簇索引是很有必要的。
  • 宏基因組實戰7. bwa序列比對, samtools查看, bedtools豐度統計
    前情提要如果您在學習本教程中存在困難，可能因為缺少背景知識，建議先閱讀本系統前期文章測試數據劉博士幫助把測試數據建立了一個百度雲同步共享文件夾
  • 【科研工具】做序列比對,這個工具最好用!
    第二步，在方框中輸入序列並提交。如果比對人類基因，在方框上方的Genome 中下拉選擇Human。依具體研究需要，也可以選擇其它物種，甚至勾上Search all選中所有物種。第三步，預覽比對結果。Submit後幾乎立即就可得到BLAT比對結果，並以清晰明了的方式將該目標序列與人類基因組資料庫的匹配情況顯示出來。在本實例中，可以看到該序列BLAT產生了3個人類基因的匹配位置，按其分值大小上下排列。
  • SEO的索引和抓取是什麼意思,外貿自建站如何優化索引和抓取?
    抓取和索引這兩件事就是SEO領域中簡單而又重要的觀念，熟悉了解它們之後便可以優化搜尋引擎蜘蛛抓取、索引你的網站。 理解抓取 ( Crawl ) 、索引 ( Index ) 搜尋引擎運作原理我們可以簡單說為： 抓取 ( Crawl) – > 演算、建立索引到搜尋引擎上 ( Index ) – > 供查詢、使用 抓取 ( Crawl) 便是指搜尋引擎捕捉你網站上的資料的行為，包括網站的關鍵字、內容、反向連結等等
  • MySQL全文索引、聯合索引、like查詢、json查詢速度大比拼
    下面為該表建立一個聯合索引(本來想建一個type-del-is_leaf_outline的索引，但是outline欄位太長限制，所以只加type-del-is_leaf_的聯合索引ALTERTABLE tmp_test_course ADDKEY`type-del-is_leaf` (`type`,`del`
  • 【陪你學·生信】九、多序列比對-Multiple Sequence Alignment(MSA)
    ，這樣可以提供很多信息；（3）多序列比對選用10-15條序列開始比對（如果10條的結果不錯，又想再加別的序列進行分析也可以。首先進行序列兩兩比對，構建距離矩陣→基於兩兩比對距離矩陣，由關係近的序列逐漸加入關係遠的序列構建引導樹guide tree→進行多序列比對。由此可見，比對的準確性高度依賴於一開始的兩兩比對，比較適用於親緣關係較近的序列。Clustal Omega中改進的新兩兩比對和建guide tree算法使Omega在W的基礎上，速度、準確度和數據處理量上與所提升。
  • 為什麼 MongoDB 索引選擇B樹,而 Mysql 索引選擇B+樹?
    針對我們這個問題的最核心的特點如下：（1）多路，非二叉樹（2）每個節點既保存索引，又保存數據（3）搜索時相當於二分查找其他的基本上都是一些常見的數據結構，假定都已經了解了B樹相關的結構。（2）B+樹的數據保存在葉子節點，查詢時間複雜度固定是O(log(n))第二：區間查找（1）B樹每個節點 key 和 data 在一起，則無法區間查找。
  • Elasticsearch開始的第一步索引index
    讓我們建立一個員工目錄假設我們剛好在Megacorp工作，這時人力資源部門出於某種目的需要讓我們創建一個員工目錄，這個目錄用於促進人文關懷和用於實時協同工作，所以它有以下不同的需求：數據能夠包含多個值的標籤、數字和純文本。檢索任何員工的所有信息。支持結構化搜索，例如查找30歲以上的員工。
  • SEO如何診斷Google索引狀況?分析網站架構頁面健康與否
    為什麼要檢查Google索引的狀況 Google對於網站索引建立的狀況，可以很直接的反映出網站與搜尋引擎之間的關係是否健康，通常過多或過少的索引建立，代表你的網站SEO是有問題的。 假設你的網站有100 個頁面，但通過觀察，你發現在 Google 的搜尋引擎上，你只有被建立了 20頁的索引，那代表著有 80頁憑空不見了，這時候你的 SEO明顯是出了狀況的。 怎麼檢查網站被索引了幾頁?
  • 從原理到優化,深入淺出資料庫索引 - 計算機java編程
    使用複合索引時遵循最左前綴集合（4）唯一索引或者非唯一索引（5）空間索引空間索引是對空間數據類型的欄位建立的索引，MYSQL中的空間數據類型有4種，分別是GEOMETRY、POINT、LINESTRING、POLYGON。MYSQL使用SPATIAL關鍵字進行擴展，使得能夠用於創建正規索引類型的語法創建空間索引。