RNAseq005 轉(zhuǎn)錄組入門(5)：序列比對

1 比對的目的

1.1 差異表達基因
只需要確定不同的read計數(shù)，可以用bowtie, bwa這類比對工具抱怔，或者是salmon這類align-free工具，并且后者的速度更快。
1.2 尋找新剪切位點
但是如果需要找新的isoform，或RNA的可變剪切之景，看外顯子使用差異的話，需要TopHat, HISAT2或者是STAR這類工具用于找到剪切位點膏潮。

2 RNA-Seq數(shù)據(jù)比對和DNA-Seq數(shù)據(jù)比對有什么差異

RNA-Seq不同于DNA-Seq锻狗，DNA在轉(zhuǎn)錄成mRNA的時候會把內(nèi)含子部分去掉，所以mRNA反轉(zhuǎn)的cDNA如果比對不到參考序列戏罢，會被分開屋谭，重新比對一次脚囊，判斷中間是否有內(nèi)含子龟糕。

3 工具抉擇

在2016年的一篇綜述A survey of best practices for RNA-seq data analysis，提到目前有三種RNA數(shù)據(jù)分析的策略悔耘。那個時候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已經(jīng)被它的作者推薦改用HISAT進行替代讲岁。
最近的Nature Communication發(fā)表了一篇題為的Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis的文章—被稱之為史上最全RNA-Seq數(shù)據(jù)分析流程，文章中在基于參考基因組的轉(zhuǎn)錄本分析中所用的工具衬以，是TopHat,HISAT2和STAR缓艳，結(jié)論就是HISAT2找到junction正確率最高，但是在總數(shù)上卻比TopHat和STAR少看峻。從這里可以看出HISAT2的二類錯誤(納偽）比較少阶淘，但是一類錯誤（棄真）就高起來。
就唯一比對而言互妓，STAR是三者最佳的溪窒，主要是因為它不會像TopHat和HISAT2一樣在PE比對不上的情況還強行把SE也比對到基因組上坤塞。而且在處理較長的read和較短read的不同情況，STAR的穩(wěn)定性也是最佳的澈蚌。
就速度而言摹芙，HISAT2比STAR和TopHat2平均快上2.5~100倍。

3 index下載

3.1 hisat2主頁下載

http://daehwankimlab.github.io/hisat2/

image.png

http://daehwankimlab.github.io/hisat2/download/

image.png

3.2 http://ccb.jhu.edu/software.shtml

http://ccb.jhu.edu/software/hisat2/manual.shtml
ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/

image.png

cd referece && mkdir index && cd index
# hg19 index
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
# hg38 index
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg38.tar.gz

將上述鏈接復制放入迅雷宛瞄，十分鐘搞定下載

4.比對alignment

4.1 說明書：https://daehwankimlab.github.io/hisat2/manual/

4.2 用法

hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | –sra-acc <SRA accession number>} [-S <hit>]

主要參數(shù)：
-x <hisat2-idx>: 參考基因組索引文件的前綴
如下圖浮禾，hisat2 hg19index解壓后如下，-x參數(shù)只需要index的前綴份汗，也就是genome盈电，所以直接在路徑后加入genome即可，如果直接使用路徑裸影，或者使用正則表達式genome.*,或是cat合并成一個genome都是不可取的

hisat2 hg19index

-1 <m1>: 雙端測序結(jié)果的第一個文件挣轨。若有多組數(shù)據(jù)，使用逗號將文件分隔轩猩。Reads的長度可以不一致卷扮。
-2 <m2>: 雙端測序結(jié)果的第二個文件。若有多組數(shù)據(jù)均践，使用逗號將文件分隔晤锹，并且文件順序要和-1參數(shù)對應。Reads的長度可以不一致彤委。
-U <r>: 單端數(shù)據(jù)文件鞭铆。若有多組數(shù)據(jù)，使用逗號將文件分隔焦影〕邓欤可以和-1、-2參數(shù)同時使用斯辰。Reads的長度可以不一致舶担。
–sra-acc <SRA accession number>: 輸入SRA登錄號，比如SRR353653彬呻，SRR353654衣陶。多組數(shù)據(jù)之間使用逗號分隔。HISAT將自動下載并識別數(shù)據(jù)類型闸氮，進行比對剪况。
-S <hit> : 指定輸出的SAM文件。

輸入選項：
-q:輸入文件為FASTQ格式蒲跨。FASTQ格式為默認參數(shù)译断。
-qseq:輸入文件為QSEQ格式。
-f:輸入文件為FASTA格式或悲。
-r:輸入文件中孙咪，每一行代表一條序列藏姐，沒有序列名和測序質(zhì)量等。選擇此項時该贾，–ignore-quals參數(shù)也會被選擇羔杨。
-c:此參數(shù)后是直接比對的序列，而不是包含序列的文件名杨蛋。序列間用逗號隔開兜材。選擇此項時，–ignore-quals參數(shù)也會被選擇逞力。
-s/–skip <int>:跳過輸入文件中前條序列進行比對曙寡。
-u/–qupto <int>:只使用輸入文件中前條序列進行比對，默認是沒有限制寇荧。
-5/–trim5 <int>:比對前去除每條序列5’端個堿基
-3/–trim3 <int>:比對前去除每條序列3’端個堿基
–phred33:輸入的FASTQ文件堿基質(zhì)量值編碼標準為phred33举庶，phred33為默認參數(shù)。

文章中提到：Samples 9-15 are mRNA-seq to determine effect of AKAP95 knockdown in human 293 cells (9-11) or mouse ES cells (12-15).

# -t參數(shù)可以顯示時間
# -x參數(shù)指定index位置：/mnt/e/bioinfo/reference/gtf/hisat2_index/hg19/揩抡，需要在路徑后添加index的前綴genome户侥，否則會報錯
# -1，-2（雙向測序）參數(shù)指定數(shù)據(jù)存放位置：/mnt/e/bioinfo/rna_seq/ENA_data/
# -S指定輸出sam的路徑峦嗤，僅有三個是人基因組數(shù)據(jù)
for i in {56..58};do hisat2 -t -x /mnt/e/bioinfo/reference/gtf/hisat2_index/hg19/genome  -1 /mnt/e/bioinfo/rna_seq/ENA_data/SRR35899${i}_1.fastq.gz  -2  /mnt/e/bioinfo/rna_seq/ENA_data/SRR35899${i}_2.fastq.gz   -S  /mnt/e/bioinfo/rna_seq/rna_seq/aligned/SRR35899$i.sam;done

我的筆記本是8G內(nèi)存蕊唐，差不多是半小時跑完一個數(shù)據(jù)，以下為比對運行結(jié)果

image.png

Hisat2 bowtie2比對結(jié)果解讀（Hisat2 Alignment summary）

第一部分描述的是pair-end模式下的一致比對結(jié)果
1820962 (6.31%) aligned concordantly 0 times 是什么意思烁设？aligned concordantly就是read1和read2同時合理的比對到了基因組/轉(zhuǎn)錄組上替梨。
24723672 (85.68%) aligned concordantly exactly 1 time，exactly 1 time 就是只有一種比對結(jié)果装黑。>1 times 就是read1和read2可以同時比對到多個地方副瀑。
第二部分，pair-end模式下不一致的比對結(jié)果
1820962 pairs aligned concordantly 0 times; of these: 90371 (4.96%) aligned discordantly 1 time恋谭，concordantly 0 times就是read1和read2不能同時合理的同時比對到基因組/轉(zhuǎn)錄組上
aligned discordantly 1 time 最難理解糠睡，discordantly比對就是read1和read2都能比對上，但是不合理箕别，
1.比對方向不對铜幽，pair-end測序的方向是固定的滞谢；
2.read1和read2的插入片段長度是有限的
第三部分就是對剩余reads（既不能concordantly串稀，也不能discordantly 1 time）的單端模式的比對，沒什么好說的狮杨。

5 SAMtools

Samtools是一個用于操作sam和bam格式文件的應用程序集合母截，具有眾多的功能。它從SAM（序列比對/映射）格式導入和導出橄教，進行排序清寇，合并和索引喘漏，并允許快速檢索任何區(qū)域中的讀數(shù)。SAM和BAM格式的比對文件主要由bwa华烟、bowtie2翩迈、tophat和hisat2等序列比對工具產(chǎn)生，用于記錄測序reads在參考基因組上的映射信息盔夜。其中负饲，BAM格式文件是SAM文件的的二進制格式，占據(jù)內(nèi)存較小且運算速度更快喂链。

5.1 語法

Program: samtools (Tools for alignments in the SAM format)
Version: 1.10 (using htslib 1.10)
Usage:   samtools <command> [options]
Commands:
  -- Indexing 索引
     dict           創(chuàng)建一個序列字典文件
     faidx          建立Fasta索引文件返十，以.fai后綴結(jié)尾
     fqidx          建立Fastq索引文件
     index          建立sam索引文件，需對bam文件進行排序后才能構(gòu)建索引
  -- Editing 編輯
     calmd          重新計算MD/NM標簽和'='基數(shù)
     fixmate        為以名稱排序定位的alignment填入配對坐標
     reheader       替換BAM頭注釋
     targetcut      切割fosmid區(qū)域（僅適用于fosmid池）
     addreplacerg   添加或替換RG標簽
     markdup        重復標記
  -- File operations 文件操作
     collate        根據(jù)read name信息將bam文件進行隨機打散和分組
     cat            將多個bam文件合并為一個bam文件椭微，無需對sam文件進行排序
     merge          將多個已經(jīng)sort的bam文件進行合并
     mpileup        對參考基因組每個位點做堿基堆積洞坑，用于call SNP和INDEL。主要是生成BCF蝇率、VCF文件或者pileup一個或多個bam文件
     sort           對比對后的bam文件進行排序迟杂，默認按堿基位置坐標進行排序
     split          按讀取組分割文件
     quickcheck     快速檢查SAM/BAM/CRAM是否完整
     fastq          將BAM轉(zhuǎn)換為FASTQ
     fasta          將BAM轉(zhuǎn)換為FASTA
  -- Statistics 統(tǒng)計
     bedcov         統(tǒng)計給定region區(qū)間內(nèi)reads的深度，每個堿基的深度疊加在一起
     coverage       統(tǒng)計每個堿基位點的測序深度及百分比
     depth          統(tǒng)計每個堿基位點的測序深度本慕，注意計算前必須對bam文件排序并構(gòu)建索引
     flagstat       統(tǒng)計bam文件中read的比對情況
     idxstats       統(tǒng)計bam索引文件里的比對信息
     phase          雜合子項
     stats          對bam文件做詳細統(tǒng)計, 統(tǒng)計結(jié)果可用mics/plot-bamstats作圖
  -- Viewing 查看
     flags          查看不同flag值的含義
     tview          直觀地顯示reads比對到參考基因組的情況逢慌，類似于基因組瀏覽器
     view           SAM<->BAM<->CRAM轉(zhuǎn)換
     depad          將填充的BAM轉(zhuǎn)換為未填充的BAM

5.2 Samtools常用命令

最常用的三板斧就是格式轉(zhuǎn)換，排序间狂，索引

轉(zhuǎn)換：samtools view -S SRR3589912.sam -b > SRR3589912.bam（-S是最新版的samtools為了兼容以前的版本寫的）
排序：samtools sort SRR3589912.bam -o SRR3589912_sorted.bam
索引：samtools index SRR3589912_sorted.bam
比對后的sam文件應該先進行格式轉(zhuǎn)換攻泼，接著是排序，最后根據(jù)排序的文件建立索引文件

5.2.1 格式轉(zhuǎn)換 samtools view

5.2.1.1 用法

Usage: samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]

Options:
  -b             輸出BAM
  -C             輸出CRAM (requires -T)
  -1             使用快速BAM壓縮(implies -b)
  -u             未壓縮的BAM輸出(implies -b)
  -h             在SAM輸出中包括頭注釋
  -H             僅打印SAM標頭（無比對信息）
  -c             僅打印匹配記錄的數(shù)量
  -o FILE        輸出文件名[stdout]
  -U FILE        輸出未過濾的reads到文件 [null]
  -t FILE        列出參考序列的名稱和長度[null]
  -X             包含通用的索引文件
  -L FILE        僅包括與該BED FILE重疊的Reads [null]
  -r STR         僅包含在STR組中的Reads [null]
  -R FILE        僅包含F(xiàn)ILE [null]中列出的讀取組的Reads
  -d STR:STR     僅包含帶有標簽STR和關聯(lián)值STR的Reads [null]
  -D STR:FILE    僅包含帶有標簽STR以及FILE [null]中列出的相關值的Reads
  -q INT         僅包含映射質(zhì)量> = INT [0]的Reads
  -l STR         僅包括對庫STR的Reads[null]
  -m INT         僅包含消耗CIGAR操作數(shù)的Reads
  -f INT         僅包括存在INT中所有FLAG的Read[0]
  -F INT         僅包括在INT中不存在任何FLAGS的Read[0] [0]
  -G INT         僅排除當前FLAGs為INT的Read [0]
  -s FLOAT       子樣本讀燃蟆（給定INT.FRAC選項值忙菠，0。FRAC是要保留的模板/讀取對的分數(shù)纺弊； INT部分設置種子）
  -M             使用多區(qū)域迭代器（提高速度牛欢，刪除重復項并按文件中的順序輸出讀取結(jié)果）
  -x STR         讀取標記以剝離[null]
  -B             折疊CIGAR操作
  -?             幫助，包括有關區(qū)域規(guī)范的注釋
  -S             忽略（自動檢測輸入格式）
  --no-PG  do not add a PG line
      --input-fmt-option OPT[=VAL]                      以OPTION或OPTION = VALUE的形式指定一個輸入文件格式選項
  -O, --output-fmt FORMAT[,OPT[=VAL]]...                指定輸出格式（SAM淆游，BAM傍睹，CRAM）
      --output-fmt-option OPT[=VAL]                     以OPTION或OPTION = VALUE的形式指定單個輸出文件格式
  -T, --reference FILE                參考序列FASTA FILE 
  -@, --threads INT                   要使用的其他線程數(shù)[0]
      --write-index                   自動索引輸出文件[關閉]
      --verbosity INT                 設置詳細程度

5.2.1.2 samtools的view不僅可以進行格式轉(zhuǎn)換，還可以進行數(shù)據(jù)的提取

# 提取1號染色體上1234~123456區(qū)域的以對read
samtools view SRR3589957_sorted.bam chr1:1234-123456 | head

5.2.1.3 SAM文件中FLAG值的理解

FLAG列在SAM文件的第二列犹菱，這是一個很重要的列拾稳，包含了很多mapping過程中的有用信息

image.png

SAM文件的flag標簽包含0x4

### 小寫的f是提取，大寫的F是過濾
samtools view -f4 sample.bam sample.unmapped.sam

對于PE數(shù)據(jù)腊脱，在未比對成功的測序數(shù)據(jù)可以分成3類：

僅reads1沒有比對成功访得，該提取條件包括：
該read是read1，對應于二進制FLAG的第7位嘁灯，該位取1涯捻，其十進制值為64；
該read未成功比對到參考基因組烟勋，對應于二進制FLAG的第3位搜骡，該位取1拂盯，其十進制值為4；
另一配對read成功比對到參考基因組记靡，對應于二進制FLAG的第4位磕仅，該位取0，其十進制值為8簸呈；

 # 對于取1的位點采用提取的策略榕订，用-f參數(shù)，值設為64+4=68
 # 對于取0的位點采取過濾的策略蜕便，用-F參數(shù)劫恒，值設為8
 samtools view -u -f 68 -F 8 alignments.bam >read1_unmap.bam

僅reads2沒有比對成功，該提取條件包括：
該read是read2轿腺，對應于二進制FLAG的第8位两嘴，該位取1，其十進制值為128族壳；
該read未成功比對到參考基因組憔辫，對應于二進制FLAG的第3位，該位取1仿荆，其十進制值為4贰您；
另一配對read成功比對到參考基因組，對應于二進制FLAG的第4位拢操，該位取0锦亦，其十進制值為8；

# 對于取1的位點采用提取的策略令境，用-f參數(shù)杠园，值設為128+4=132
# 對于取0的位點采取過濾的策略，用-F參數(shù)舔庶，值設為8
samtools view -u -f 132 -F 8 alignments.bam >read2_unmap.bam

兩端reads都沒有比對成功抛蚁，該提取條件包括：
該read未成功比對到參考基因組，對應于二進制FLAG的第3位惕橙，該位取1瞧甩，其十進制值為4；
另一配對read未成功比對到參考基因組吕漂，對應于二進制FLAG的第4位亲配，該位取1尘应，其十進制值為8惶凝；

# 對于取1的位點采用提取的策略吼虎，用-f參數(shù)，值設為4+8=12
$ samtools view -u -f 12 alignments.bam >pairs_unmap.bam

5.2.2 排序 samtools sort

Usage: samtools sort [options...] [in.bam]
Options:
  -l INT     設置輸出文件壓縮等級苍鲜。0-9思灰，0是不壓縮，9是壓縮等級最高混滔。不設置此參數(shù)時洒疚，使用默認壓縮等級
  -m INT     設置每個線程運行時的內(nèi)存大小，可以使用K坯屿，M和G表示內(nèi)存大小
  -n         設置按照read名稱進行排序
  -t TAG     Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
  -o FILE    設置最終排序后的輸出文件名
  -T PREFIX  設置臨時文件的前綴 PREFIX.nnnn.bam
  --no-PG不添加PG行
      --input-fmt-option OPT [= VAL]                      以OPTION或OPTION = VALUE的形式指定一個輸入文件格式選項
  -O油湖，--output-fmt FORMAT [，OPT [= VAL]] ...             指定輸出格式（SAM领跛，BAM乏德，CRAM）
      --output-fmt-option OPT [= VAL]                     以OPTION或OPTION = VALUE的形式指定單個輸出文件格式選項
      --reference FILE                                    參考序列FASTA FILE [null]
  -@, --threads INT               設置排序和壓縮是的線程數(shù)量，默認是單線程
      --verbosity INT               Set level of verbosity

對bam文件進行排序吠昭，不能對sam文件進行排序喊括。以leftmost coordinates的方式對比對結(jié)果進行排序，或者使用-n參數(shù)以read名稱進行排序矢棚。將會添加適當?shù)腀HD-SO排序順序標頭標簽或者如果有必要的話郑什，將會更新現(xiàn)存的一個排序順序標頭標簽。sort命令的輸出默認是標準輸出寫入蒲肋，或者使用-o參數(shù)時蘑拯，指定bam文件輸出名。sort命令還會在內(nèi)存不足時創(chuàng)建臨時文件tmpprefix.%d.bam兜粘。

samtools的排序方式有兩種（常用）

# 默認方式强胰，按照染色體的位置進行排序
 for i in {56..58};do samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam;done
# 參數(shù)-n則是根據(jù)read名進行排序。
 for i in {56..58};do samtools sort -n SRR35899${i}.bam -o SRR35899${i}_nsorted.bam;done

image.png

bam文件是Sam 文件的二進制壓縮格式妹沙，保留了與sam 完成相同的內(nèi)容信息偶洋。SAM/BAM 文件可以是未排序的，但是按照坐標(coodinate)排序可以線性的監(jiān)控數(shù)據(jù)處理過程距糖。samtools可以用來轉(zhuǎn)化bam/sam文件玄窝，可以merg，sort aligment,可以去除duplicate,可以call snp及indels.
BAM 文件是壓縮的二進制文件悍引，對文件內(nèi)容排序之后相似的內(nèi)容排在一起恩脂，使得文件壓縮比提高了，因此排序之后的 BAM 文件變小了趣斤，相對應的 SAM 文件就是純文本文件俩块，對 SAM 文件進行排序就不會改變文件大小。而且由于 RNA-seq 中由于基因表達量的關系，RNA-seq 的數(shù)據(jù)比對結(jié)果 BAM 文件使用 samtools 進行 sort 之后文件壓縮比例變化會比 DNA-seq 更甚玉凯。

5.2.3 索引 samtools index

Usage: samtools index [-bc] [-m INT] <in.bam> [out.index]
Options:
  -b       創(chuàng)建bai索引文件势腮，未指定輸出格式時，此參數(shù)為默認參數(shù)
  -c       創(chuàng)建csi索引文件漫仆，默認情況下捎拯，索引的最小間隔值為2^14，與bai格式一致
  -m INT   創(chuàng)建csi索引文件盲厌，最小間隔值2^INT
  -@ INT   設置線程數(shù)[無]

為了能夠快速訪問bam文件署照，可以為已經(jīng)基于坐標排序后bam或者cram的文件創(chuàng)建索引，生成以.bai或者.crai為后綴的索引文件吗浩。,必須使用排序后的文件建芙。另外，不能對sam文件使用此命令懂扼。如果想對sam文件建立索引岁钓，那么可以使用tabix命令創(chuàng)建。在使用與region參數(shù)相關的其它samtools命令時微王，需要創(chuàng)建索引

# -i 參數(shù)直接用 {}就能代替管道之前的標準輸出的內(nèi)容
ls *sorted.bam | xargs -i samtools index {}
# 或是ls *.bam | xargs -n1  samtools index

xargs命令可參考：http://www.reibang.com/writer#/notebooks/45449543/notes/76201605/preview

image.png

報錯：默認排序可以建立index屡限，但是按reads名排序（sort -n）無法建立index

Reference

http://www.reibang.com/p/681e02e7f9af
http://www.reibang.com/p/68f6e35fa4a2
https://ming-lian.github.io/2019/02/07/Advanced-knowledge-of-SAM/