轉(zhuǎn)錄組（5）：序列比對

參考http://www.itdecent.cn/p/681e02e7f9af
http://www.itdecent.cn/p/93f96e7538da
任務：

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

1. 比對軟件

2. HISAT2的使用

為什么要用index？官網(wǎng)有描述：為了用整個index代表整個基因組，HISAT2 用小的index覆蓋了整個基因組，每個index覆蓋了56 Kbp的范圍，覆蓋整個人類基因組需要55,000 indexes，這些index結(jié)合其他策略可以快速準確的比對序列。

#index 下載
nohup wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz &
tar -zxvf *.tar.gz #解壓
# 刪除壓縮包
rm -rf *.tar.gz
  • hisat2 -h查看幫助文件：
Usage: 
  hisat [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <sam>]

  <bt2-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'.
  • 參數(shù)說明：

-x 指定index文件名
-1 <m1> 雙端測序第一個文件
-2 <m2> 雙端測序第二個文件
-U 單端測序文件
--sra-acc 登錄號
-S 輸出文件為sam格式
-q 輸入文件為FASTQ .fq/.fastq格式

-比對到參考基因組,RNA-Seq數(shù)據(jù)是從 SRR3589957 ~ SRR3589962，6個樣本的PE數(shù)據(jù)，也就是有6次循環(huán)，可以寫腳本，也可以直接在命令行里運行

for i in `seq 56 62`
do
    hisat2 -t -x /opt/NfsDir/UserDir/qin/qin/Data/annotation/hg19/genome -1 /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/SRR35899${i}.sra_1.fastq.gz -2 /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/SRR35899${i}.sra_2.fastq.gz -S /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.sam &
done

運行時間如下：

Paste_Image.png
  • 很耗CPU啊！用的服務器！
Paste_Image.png
  • 結(jié)果：
Paste_Image.png

SAM（sequence Alignment/mapping)數(shù)據(jù)格式是目前高通量測序中存放比對數(shù)據(jù)的標準格式，當然他可以用于存放未比對的數(shù)據(jù)。目前處理SAM格式的工具主要是SAMTools。samtools功能眾多，在本次作業(yè)中，我們主要學會將sam文件轉(zhuǎn)換為bam文件，并對bam文件進行sorted（其中有兩種排序方式N和P），最后建立索引。

Program: samtools (Tools for alignments in the SAM format)
Version: 1.3.1-58-gcbee45e (using htslib 1.3.2-228-g0c32631)

Usage:   samtools <command> [options]

Commands:
  -- Indexing
     dict           create a sequence dictionary file
     faidx          index/extract FASTA
     index          index alignment

  -- Editing
     calmd          recalculate MD/NM tags and '=' bases
     fixmate        fix mate information
     reheader       replace BAM header
     rmdup          remove PCR duplicates #移除PCR產(chǎn)生的重復序列
     targetcut      cut fosmid regions (for fosmid pool only)
     addreplacerg   adds or replaces RG tags

  -- File operations
     collate        shuffle and group alignments by name
     cat            concatenate BAMs
     merge          merge sorted alignments
     mpileup        multi-way pileup
     sort           sort alignment file
     split          splits a file by read group
     quickcheck     quickly check if SAM/BAM/CRAM file appears intact
     fastq          converts a BAM to a FASTQ #格式轉(zhuǎn)換
     fasta          converts a BAM to a FASTA

  -- Statistics
     bedcov         read depth per BED region #bed文件的測序深度
     depth          compute the depth 
     flagstat       simple stats
     idxstats       BAM index stats
     phase          phase heterozygotes
     stats          generate stats (former bamcheck)

  -- Viewing
     flags          explain BAM flags
     tview          text alignment viewer
     view           SAM<->BAM<->CRAM conversion
     depad          convert padded BAM to unpadded BAM

主要功能：
view: BAM-SAM/SAM-BAM 轉(zhuǎn)換和提取部分比對
sort: 比對排序
merge: 聚合多個排序比對
index: 索引排序比對
faidx: 建立FASTA索引，提取部分序列
tview: 文本格式查看序列
pileup: 產(chǎn)生基于位置的結(jié)果和 consensus/indel calling

# 首先將比對后的sam文件轉(zhuǎn)換成bam文件
# 利用的是samtools的view選項，參數(shù)-S 輸入sam文件；參數(shù)-b 指定輸出的文件為bam；最后重定向?qū)懭隻am文件
$ for ((i=56;i<=62;i++));do samtools view -S /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.sam -b > /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.bam;done
# 將所有的bam文件進行排序
$ nohup for (( i=58;i<=62;i++ )); do samtools sort /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.bam -o /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.sorted.bam;done
# 將所有的排序文件建立索引，索引文件.bai后綴
$ for ((i=56;i<=62;i++));do samtools index /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR35899${i}.sorted.bam;done

合在一塊跑：
for i in `seq 56 58`
do
    samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam
    samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
    samtools index SRR35899${i}_sorted.bam
done
Paste_Image.png

比對質(zhì)控(QC)

常用工具有

java -jar /opt/NfsDir/BioDir/picard-tools-1.119/picard.jar CollectMultipleMetrics \
      I=/opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR3589956.sorted.bam \
      O=/opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/multiple_metrics \
      R=/opt/NfsDir/PublicDir/reference/ucsc.hg19.fasta

統(tǒng)計bam文件：

/opt/NfsDir/BioDir/RSeQC/RSeQC-2.6.4/scripts/bam_stat.py -i /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR3589956.sorted.bam

提示出錯：


Paste_Image.png

檢查發(fā)現(xiàn)是路徑和名稱寫錯了,修改后就可以了
對bam文件進行統(tǒng)計分析


Paste_Image.png
#下載hg19_RefSeq.bed文件
https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz/download
#查看基因組覆蓋率
/opt/NfsDir/BioDir/RSeQC/RSeQC-2.6.4/scripts/read_distribution.py -i /opt/NfsDir/UserDir/qin/qin/Data/RNAseq/align/SRR3589956.sorted.bam -r /opt/NfsDir/UserDir/qin/qin/Data/annotation/hg19/hg19_RefSeq.bed
Paste_Image.png
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容