【轉錄組入門】5：序列比對

作業要求：

比對軟件很多，首先大家去收集一下，因為我們是帶大家入門，請統一用hisat2，並且搞懂它的用法。
直接去hisat2的主頁下載index文件即可，然后把fastq格式的reads比對上去得到sam文件。
接着用samtools把它轉為bam文件，並且排序(注意N和P兩種排序區別)索引好，載入IGV，再截圖幾個基因看看！
順便對bam文件進行簡單QC，參考直播我的基因組系列。

 

 

【1】選擇比對工具

Hisat2 或 STAR（本次選擇Hisat2）

 

【2】下載hisat2的index文件（就是參考基因組）

1、進入hisat2官網，主頁面右側中間位置，

右擊“H.sapiens,UCSC hg19”下面的genome，“復制鏈接地址”

2、Ubuntu終端命令行下載

# 切換到下載目錄
$ cd ~/src
$ wget  ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz

# 解壓下載后的文件
$ tar -zxvf hg19.tar.gz

補充：

（1）可以在Windows系統中用迅雷下載，這樣下載比較快

（2）若沒有現成的index，可以通過HISAT2的方法創建

# 其實hisat2-buld在運行的時候也會自己尋找exons和splice_sites，但是先做的目的是為了提高運行效率
$ extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf &
# 建立index， 必須選項是基因組所在文件路徑和輸出的前綴
$ hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19

 

 【3】正式比對

# 編寫批量比對的腳本  hisat2_align.sh
   for i in `seq 56 58`
   do
         hisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 RNA-Seq/SRR35899${i}_2.fastq.gz  -S  RNA-Seq/aligned/SRR35899${i}.sam &
   done
# 運行腳本 hisat2_align.sh
$ bash  hisat2_align.sh

# hisat2用法
$ hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <hit>]
# -x index : 參考基因組
# 雙端測序：hisat2 -x hisat2_index -1 m1 -2 m2 -S name.sam
# 單端測序：hisat2 -x hisat2_index -U r1 -S name.sam

比對結果解釋：

1、全部數據都是100%的，6.37%的數據一次都沒有比對，85.71%的數據是唯一比對，7.92%是多個比對。

2、然后6.37%的一次都沒有比對的數據，如果不按照順序來，有4.94%的比對。

3、之后把剩下的部分用單端數據進行比對的話，58.21%數據沒比對上，34.94%的數據比對一次，6.86%比對超過一次。

4、零零總總的加起來是96.47%的比對！！

 

 【4】sam格式轉換成bam格式

# 編寫腳本
  for i in `seq 56 58`
  do
      samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam  #轉換成bam文件
      samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam  #采用默認方式排序
      samtools index SRR35899${i}_sorted.bam  #生成索引
  done

 

【5】比對結果QC

QC軟件有很多，這里我們使用RSeQC

• RSeQC——http://rseqc.sourceforge.net/

• Qualimap——http://qualimap.bioinfo.cipf.es/

• Picard——http://broadinstitute.github.io/picard/

使用RSeQC進行比對結果QC

 1、安裝RSeQC

# RSeQC的安裝，需要先安裝gcc；numpy；R；Python2.7，這里比較難裝的就是numpy——可以直接利用anaconda安裝（http://www.jianshu.com/p/14fd4de54402）
# 我的環境已經配置好了，所以直接可以用pip命令安裝
$ conda install RSeQC

# 進入存放bam文件的目錄
$ cd 
# 對bam文件進行質控，其余都同樣的進行
$ bam_stat.py  -i SRR3589956_sorted.bam

2、分析並查看結果

總read數是67327865，其中能夠比對上的有49466502，所以mapped的比率是73.5%，符合人類的70%~90%的結果

 3、比對到基因組各種原件上的情況：   （基因組覆蓋率的QC）

bed文件的下載：

hg19:https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/

mm10:https://sourceforge.net/projects/rseqc/files/BED/Mouse_Mus_musculus/

兩個都選擇RefSeq

bed文件還可以用gtf文件轉換，網上也有很多寫好的腳本可以用。

開始分析：

$ read_distribution.py -i RNA-Seq/aligned/SRR3589956_sorted.bam -r reference/hg19_RefSeq.bed

分析結果：

 

【6】IGV查看比對結果

載入參考基因組，基因組注釋文件，很bam文件，看一些基因。

 

 

理論知識

RNA-Seq數據分析分為很多種：說找差異表達基因、尋找新的可變剪切……

【】找差異表達基因：單純只需要確定不同的read計數就行的話，

可用工具：bowtie, bwa這類比對工具，或者是salmon這類align-free工具，並且后者的速度更快。

【】找到新的isoform，或者RNA的可變剪切，看看外顯子使用差異的話

可用工具：TopHat, HISAT2或者是STAR這類工具用於找到剪切位點。

因為RNA-Seq不同於DNA-Seq，DNA在轉錄成mRNA的時候會把內含子部分去掉。所以mRNA反轉的cDNA如果比對不到參考序列，會被分開，重新比對一次，判斷中間是否有內含子。

文章：Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis總結：

】二類錯誤（納偽）：hisat2最少

】一類錯誤（棄真）：hisat2較高

】唯一比對：STAR最佳

】穩定性：STAR最佳

】速度：hisat2最快（HISAT2比STAR和TopHat2平均快上2.5~100倍）

 

序列比對實質：把reads和index進行比較

 

SAM（sequence alignment/mapping）數據格式

是一種序列比對格式標准， 由Sanger制定，是以TAB為分割符的文本格式。主要應用於測序序列mapping到基因組上的結果表示，當然也可以表示任意的多重比對結果

是目前高通量測序中存放比對數據的標准格式，當然他可以用於存放未比對的數據。

sam文件由：頭文件和map結果組成

】頭文件：由一行行以@起始的注釋構成

】map結果：

每個read只占一行，只是它被tab分成了很多列，一共有12列，分別記錄了：

1.read名稱

2.SAM標記

3.chromosome號

4.5′端起始位置

5.MAPQ（mapping quality，描述比對的質量，數字越大，特異性越高）

6.CIGAR字串，記錄插入，刪除，錯配以及splice junctions（后剪切拼接的接頭）

7.mate名稱，記錄mate pair信息

8.mate的位置

9.模板的長度

10.read序列

11.read質量

12.程序用標記

 

處理SAM格式的工具主要是：SAMTools，SAMTools的主要功能如下：

    view: BAM-SAM/SAM-BAM 轉換和提取部分比對

    sort: 比對排序

    merge: 聚合多個排序比對

    index: 索引排序比對

    faidx: 建立FASTA索引，提取部分序列

    tview: 文本格式查看序列

    pileup: 產生基於位置的結果和 consensus/indel calling

最常用功能：格式轉換、排序、索引

# 格式轉換

$ samtools  view  -b  test.sam  >  test.bam

# 排序

$ samtools  sort  test.bam  -o  test_sorted.bam      # 采用默認方式排序

$ samtools  sort  -n  test.bam  -o  test_sorted_n.bam     # 根據reads名排序

# 索引

$ samtools  index  test_sorted.bam

 

常用的比對質控軟件有：

1.Picard https://broadinstitute.github.io/picard/

2.RSeQC http://rseqc.sourceforge.net/

3.Qualimap http://qualimap.bioinfo.cipf.es/

 

RNA-Seq數據分析分為很多種

找差異表達基因或尋找新的可變剪切……