處理SAM、BAM你需要Samtools

. 關於SAMSam: I Want You!

sam是短序列比對默認的標準格式，由sanger制定，是以TAB為分割符的文本格式。主要應用於測序序列mapping到基因組上的結果表示，另外也可以表示其他的多重比對結果。一般把測序reads比對到參考基因組以後，通常得到的就是sam文件，全稱是sequence alignment/map format，BAM就是SAM的二進位文件(B即：Binary)
官方參考文檔http://samtools.github.io/hts-specs/SAMv1.pdf

需要知道的名詞

template：DNA/RNA序列的一部分，用它進行測序或者由原始數據拼接得到它，作為研究的模版；
【參考序列和比對上的序列共同組成的序列為template】

read： 測序儀下機得到的原始數據，這些數據依照測序時間的先後被打上不同標籤；雙端測序存在read1，read2

segment：比對到read的時候的一段連續序列或者被分開幾條子序列。一條read可以包含多條segment；

linear alignment: 一條read比對到參考序列上，可以存在插入(insert)、缺失(delete)、跳躍(skip)、剪切(clip)，但是方向不變（不能是一部分和正鏈匹配，另一部分又和負鏈匹配），sam文件中只佔用一行記錄；

chimeric alignment: 「嵌合比對」 的形成是由於一條測序read比對到基因組上時分別比對到兩個不同的區域，而這兩個區域基本沒有overlap。因此它在sam文件中需要佔用多行記錄顯示。只有第一個記錄被稱作"representative",其他的都是"supplementary"【Chimeric reads are also called split reads】

嵌合 reads 與嵌合/融合基因（重組由來源與功能不同的基因序列剪接而形成的雜合基因）並不完全一樣，RNA-seq中的chimeric read或許可以說明有融合基因存在，但在基因組中一般作為結構變異的證據

【下面👇的圖1是原始數據(其中Coor是坐標的簡寫；ref是參考序列)；圖2是使用bwa、bowtie2、hisat2等軟體比對的結果sam文件】
【其中r001/1與r001/2是paired end read； r003是嵌合read，有兩個記錄】

read alignment： 不管黑貓(linear alignment)、白貓(chimeric alignment)，抓住老鼠(能夠完整表示一個read)就是好貓(read alignment)

multiple mapping:  假如有一個短序列，他在比對的時候看到哪哪都有他的影子，這種就是受重複區域影響比較大，所以read越長出現這種的可能性越低。一般指定首先匹配上的為最優匹配結果primary。

坐標系統 ：

1-based coordinate system ：序列的第一個鹼基設為數字1，如SAM,VCF,GFF,wiggle格式

0-based coordinate system ：序列的第一個鹼基設為數字0，如BAM, BCFv2, BED, PSL格式

SAM要處理好的問題：

非常多的序列（read)，mapping到多個參考基因組（reference）上；

同一條序列，分多段（segment）比對到參考基因組上；

各種結構化信息表示，包括錯配、刪除、插入等比對信息；

具體的Sam格式：

花花之前有介紹，並且在Linux考試題中也有所提及

sam格式的頭文件可有可無，主要有@HD，說明符合標準的版本、對比序列的排列順序；@SQ，參考序列說明；@RG，比對上的序列（read）說明；@PG，使用的程序說明；@CO，任意的說明信息

. 關於Samtools

一般測序文件都比較大，生成的單個sam文件從幾百M到十幾G不等，為了學習samtools的使用首先要有sam文件

#
ascp -v -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -k 1 -T -l500m [e-mail]:/sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR3589956/SRR3589956.sra 
./raw_data
#
fastq-dump --gzip --split-3 *.sra
#
mkdir -p genome/hg19  && cd genome/hg19
for i in $(seq 1 22) X Y M; #指定染色體號下載
do echo $i;
# 一般下載UCSC的數據
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
done
gunzip *.gz
for i in $(seq 1 22) X Y M; # 指定染色體號拼成整體
do cat chr${i}.fa >> hg19.fasta
done
rm -fr chr*.fasta
#
mkdir -p ~/reference/genome/hg19 && cd ~/reference/genome/hg19
nohup time bwa index -a bwtsw -p hg19 ~/reference/genome/hg19/hg19.fa 1>hg19.bwa_index.log 2>&1 &
#
bwa mem -t 10 -M ~/reference/index/bwa/hg19/hg19 SRR3589956_1.fastq.gz SRR3589956_2.fastq.gz 1>SRR3589956.sam 2>SRR3589956.bwa.align.log
done

1. View

作用：轉換sam與bam，對bam進行排序（sort）和提取的操作

sam轉bam，-S指輸入文件格式（不加-S默認輸入是bam），-b指定輸出文件（默認輸出sam）

samtools view -Sb SRR3589956.sam >SRR3589956.bam

【如果要bam轉sam，-h設置輸出sam時帶上頭注釋信息】

samtools view -h SRR3589956.bam > SRR3589956.sam

過濾功能-F：後接flag數字，常用有4（表示序列沒比對上）、8（配對的另一條序列，即mate序列沒比對上）以及12（兩條序列都沒比對上）。加上-F就表示過濾掉這些情況

#
samtools view -bF4 SRR3589956.bam>SRR3589956.F4.bam
#
samtools view -bF12 SRR3589956.bam>SRR3589956.F12.bam

如果是-f，就是提取指定flag的序列

#
samtools view -bf4 SRR3589956.bam>SRR3589956.f4.bam

轉換後的bam文件大小是2.8G，可以看到，相比於sam文件，節省了四分之三的空間；

另外得到的SRR3589956.F4.bam是2.7G，說明本來bam文件中比對上的reads就很多了，才過濾掉了不到四十分之一的reads

2. Sort

作用：對bam進行排序，一些軟體需要使用排序後的bam，拿到bam先按照比對位置的順序排一下，百利無一害

samtools sort -@ 10 -o SRR3589956.sorted.bam SRR3589956.bam

sort後的文件SRR3589956.sorted.bam大小是1.5G

3. merge

作用：將兩個及以上的sort過的bam文件融合成一個bam文件

這裡就一個示例文件就不演示了

samtools merge xxx.merged.bam xxx.sorted.bam

4. index

作用：對bam文件構建索引，產生.bai文件，方便以後的快速處理
bam文件進行排序sort後，才能進行index，否則報錯；
要顯示比對結果時，比如用IGV導入bam，就需要有.bai的存在

samtools index SRR3589956.sorted.bam

插播一條

當有多個bam文件時，一般思路就是對每一個bam進行sort、index後，再merge成一個整體merged.bam，然後對merged.bam再進行sort、index，才算能用了，得到最終結果應該是是sorted.merge.bam

5. faidx

作用：對fasta文件建立提取索引，索引文件後綴是.fai。利用索引文件可以快速提取fasta文件中的某些序列

samtools faidx hg19.fasta

6. tview

作用：直觀顯示reads比對到基因組的情況，與IGV類似
需要先sort和index

samtools tview SRR3589956.sorted.bam hg19.fasta  

7. flagstat

作用：統計bam文件的比對結果

samtools flagstat SRR3589956.sorted.bam > SRR3589956.sorted.flagstat.txt

8. depth

作用：統計每個鹼基位點的測序深度
需要使用重定向定義輸出文件；要是用構建過索引的bam

-r：（region）加染色體號；
-q：要求測序鹼基質量最低值；
-Q：要求比對的質量最低值

samtools depth SRR3589956.sorted.bam >SRR3589956.depth

9. mpileup

作用：用於生成bcf文件，之後結合bcftools進行SNP與InDel的分析
安裝samtools時，包含了bcftools

-f：輸入有索引的參考基因組fasta；
-g：輸出到二進位的bcf格式【不使用-g，就不生成bcf格式，而是一個文本文件，統計了參考序列中每個鹼基位點的比對情況；每一行代表參考序列中某一個鹼基位點的比對結果】

samtools mpileup -f hg19.fa SRR3589956.sorted.bam > SRR3589956.mpileup.txt

結果文件大小 7.1G Aug  3 21:48 SRR3589956.mpileup.txt

結果包含6列：參考序列名、匹配位置、參考鹼基、比對上的reads數、比對的情況、比對的鹼基質量

在第5列比對具體情況中:
. 表示與參考序列正鏈匹配；
, 表示與參考序列負鏈匹配；
ATCGN 表示在正鏈不匹配；
atcgn 表示在負鏈不匹配；
* 模糊鹼基；
^ 匹配的鹼基是一個read的開始，後面的ASCII碼-33表示比對質量，再向後修飾的(.,ATCGNatcgn) 表示該read的第一個鹼基；
$ 表示一個read結束，修飾前面鹼基;
正則表達式+[0-9][ATCGNatcgn] 表示在該位點後面插入的鹼基；
正則表達式-[0-9][ATCGNatcgn] 表示該位點後面缺失的鹼基

10. bam轉fastq

作用：方便提取出一段比對到參考序列的reads進行分析

利用軟體：http://www.hudsonalpha.org/gsl/information/software/bam2fastq

11. rmdup

作用：將測序數據中由於PCR duplicate得到的reads去掉，只保留比對質量最高的reads

samtools rmdup 默認只對PE數據進行處理

12. idxstats

作用：輸出一個表格，包含「序列名、序列長度、比對上的reads數、沒有比對上的reads數」；其中第四列指PE reads中的一條read能匹配到參考基因組的染色體A，另一條read不能匹配到A上

samtools idxstats

13. reheader

作用：替換bam文件的頭文件

samtools reheader