轉錄組(5):序列比對

參考http://www.reibang.com/p/681e02e7f9af
http://www.reibang.com/p/93f96e7538da
任務:

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

1. 比對軟件

2. HISAT2的使用

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

#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
  • 結果:
Paste_Image.png

SAM(sequence Alignment/mapping)數(shù)據(jù)格式是目前高通量測序中存放比對數(shù)據(jù)的標準格式吨枉,當然他可以用于存放未比對的數(shù)據(jù)仇箱。目前處理SAM格式的工具主要是SAMTools。samtools功能眾多东羹,在本次作業(yè)中剂桥,我們主要學會將sam文件轉換為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 #格式轉換
     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 轉換和提取部分比對
sort: 比對排序
merge: 聚合多個排序比對
index: 索引排序比對
faidx: 建立FASTA索引,提取部分序列
tview: 文本格式查看序列
pileup: 產(chǎn)生基于位置的結果和 consensus/indel calling

# 首先將比對后的sam文件轉換成bam文件
# 利用的是samtools的view選項冤议,參數(shù)-S 輸入sam文件斟薇;參數(shù)-b 指定輸出的文件為bam;最后重定向寫入bam文件
$ 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

比對質控(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
最后編輯于
?著作權歸作者所有,轉載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末堪滨,一起剝皮案震驚了整個濱河市,隨后出現(xiàn)的幾起案子蕊温,更是在濱河造成了極大的恐慌袱箱,老刑警劉巖遏乔,帶你破解...
    沈念sama閱讀 216,843評論 6 502
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異发笔,居然都是意外死亡盟萨,警方通過查閱死者的電腦和手機,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,538評論 3 392
  • 文/潘曉璐 我一進店門了讨,熙熙樓的掌柜王于貴愁眉苦臉地迎上來捻激,“玉大人,你說我怎么就攤上這事前计“罚” “怎么了?”我有些...
    開封第一講書人閱讀 163,187評論 0 353
  • 文/不壞的土叔 我叫張陵男杈,是天一觀的道長丈屹。 經(jīng)常有香客問我,道長势就,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 58,264評論 1 292
  • 正文 為了忘掉前任脉漏,我火速辦了婚禮苞冯,結果婚禮上,老公的妹妹穿的比我還像新娘侧巨。我一直安慰自己舅锄,他們只是感情好,可當我...
    茶點故事閱讀 67,289評論 6 390
  • 文/花漫 我一把揭開白布司忱。 她就那樣靜靜地躺著皇忿,像睡著了一般。 火紅的嫁衣襯著肌膚如雪坦仍。 梳的紋絲不亂的頭發(fā)上鳍烁,一...
    開封第一講書人閱讀 51,231評論 1 299
  • 那天,我揣著相機與錄音繁扎,去河邊找鬼幔荒。 笑死,一個胖子當著我的面吹牛梳玫,可吹牛的內(nèi)容都是我干的爹梁。 我是一名探鬼主播,決...
    沈念sama閱讀 40,116評論 3 418
  • 文/蒼蘭香墨 我猛地睜開眼提澎,長吁一口氣:“原來是場噩夢啊……” “哼姚垃!你這毒婦竟也來了?” 一聲冷哼從身側響起盼忌,我...
    開封第一講書人閱讀 38,945評論 0 275
  • 序言:老撾萬榮一對情侶失蹤积糯,失蹤者是張志新(化名)和其女友劉穎掂墓,沒想到半個月后,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體絮宁,經(jīng)...
    沈念sama閱讀 45,367評論 1 313
  • 正文 獨居荒郊野嶺守林人離奇死亡梆暮,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 37,581評論 2 333
  • 正文 我和宋清朗相戀三年,在試婚紗的時候發(fā)現(xiàn)自己被綠了绍昂。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片啦粹。...
    茶點故事閱讀 39,754評論 1 348
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖窘游,靈堂內(nèi)的尸體忽然破棺而出唠椭,到底是詐尸還是另有隱情,我是刑警寧澤忍饰,帶...
    沈念sama閱讀 35,458評論 5 344
  • 正文 年R本政府宣布贪嫂,位于F島的核電站,受9級特大地震影響艾蓝,放射性物質發(fā)生泄漏力崇。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點故事閱讀 41,068評論 3 327
  • 文/蒙蒙 一赢织、第九天 我趴在偏房一處隱蔽的房頂上張望亮靴。 院中可真熱鬧,春花似錦于置、人聲如沸茧吊。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,692評論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽搓侄。三九已至,卻和暖如春话速,著一層夾襖步出監(jiān)牢的瞬間讶踪,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 32,842評論 1 269
  • 我被黑心中介騙來泰國打工泊交, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留俊柔,地道東北人。 一個月前我還...
    沈念sama閱讀 47,797評論 2 369
  • 正文 我出身青樓活合,卻偏偏與公主長得像雏婶,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子白指,可洞房花燭夜當晚...
    茶點故事閱讀 44,654評論 2 354

推薦閱讀更多精彩內(nèi)容