Charpter_14 Sequence Alignment Maps(SAM)

SAM文件是由Tab(\t)分割，面向‘行’的文本格式文件惕橙，主要包括兩部分：

• Header部分，包括@HD,@SQ,@PG等開(kāi)頭的基本信息滨彻。
• Alignment部分觅丰，包括reads比對(duì)的所有信息。

sam是最常見(jiàn)的存放mapping數(shù)據(jù)的格式赫冬，各種組學(xué)分析只要存在mapping的步驟都會(huì)產(chǎn)生SAM格式文件用于后續(xù)進(jìn)一步的下游分析浓镜。

比對(duì)+samtools sort一個(gè)快速示例

REF=ebola.fa
R1=SRR1972739_1.fq
R2=SRR1972739_2.fq

### bwa mem $REF $R1 $R2 > bwa.sam

### bowtie 比對(duì)輸出直接sort排序產(chǎn)生bam文件
bowtie2 -x $REF -1 $R1 -2 $R2 |samtools sort -@ 10 >bowtie2_output.sorted.bam

SAM header

header部分由@起始，包括有：

• HD：表示排序情況劲厌，unsorted或者已經(jīng)按coordinate/names排過(guò)序
• SQ：比對(duì)的參考序列的染色體信息膛薛，姓名/多長(zhǎng)。
• PG：運(yùn)行命令的具體參數(shù)信息补鼻。
• RG：每個(gè)read所在的group的信息哄啄。通常包含測(cè)序平臺(tái)/測(cè)序文庫(kù)/樣本的ID等信息(header里最為重要的一個(gè)信息).ID(Read group identifier), LB(library), SM(Sample)

@RG信息：如果原來(lái)樣本的測(cè)序深度比較深，一般會(huì)按照不同的lane分開(kāi)比對(duì)风范，最后再合并到一起咨跌。那么這個(gè)時(shí)候會(huì)有多個(gè)RG，里面記錄了不同的lane硼婿，甚至測(cè)序文庫(kù)的信息锌半，唯一不變的是SM的sample信息，這樣合并后才能正確處理寇漫。

比對(duì)信息(Alignment部分)

image

image

Column 1. QNAME(query name)：fastq文件里的read ID

Column 2. FLAG信息

• 記錄了許多有關(guān)read比對(duì)情況的信息刊殉，轉(zhuǎn)換成二進(jìn)制碼后發(fā)現(xiàn)其包括12位，十進(jìn)制取值范圍就是0~2048州胳，2的12次方记焊。

• 對(duì)于未知的FLAG，可以samtools flags 141
• 對(duì)于FLAG=77時(shí)虎囚，需要轉(zhuǎn)換為二進(jìn)制序列：77=000001001101=1+4+8+64即表示：PE read角塑，比對(duì)不上參考序列，另一端也未比對(duì)上參考序列淘讥，它是PE測(cè)序read1
• 常見(jiàn)的FLAGS：
  • 其中一條reads沒(méi)有map上： 73, 133, 89 121, 165, 181, 101, 117, 153, 185, 59, 137
  • 兩條reads都未map上：77圃伶，141
  • 比對(duì)上了，方向也對(duì)，也在插入大小(insert size)內(nèi)： 99, 147, 83, 163
  • 比對(duì)上了窒朋，也在插入大小(insert size)內(nèi)搀罢， 但是反向不對(duì):67, 131, 115, 179
  • 單一配對(duì)，就是插入大小(insert size)不對(duì): 81, 161, 97, 145, 65, 129, 113, 177

Column 3/4. RNAME(reference name)/POS(position)：

• RNAME:表示參考序列的名字
• POS：比對(duì)上的坐標(biāo)位置侥猩，注意不管是正向還是反向榔至，只記錄比對(duì)上的最左邊的坐標(biāo)位置，坐標(biāo)是以1為起始

Column 5/6. MAPQ(mapping quality)和CIGAR(compact idiosyncratic gapped alignment representation)

• MAPQ: 比對(duì)質(zhì)量值欺劳，比對(duì)到參考序列上此位置的可靠程度唧取，轉(zhuǎn)化為-10logP。40/10=4 ——> 10^-4=1/10000划提，即比錯(cuò)的概率小于0.0001枫弟。
• CIGAR：reads比對(duì)到參考序列的具體細(xì)節(jié)情況：
  • M:(完全匹配或SNP單堿基錯(cuò)配)
  • I:序列插入，包含潛在的insertion變異
  • D:序列刪除鹏往，包含潛在的deletion變異
  • S:軟跳過(guò),跳過(guò)reads中部分序列媒区，不改變r(jià)ead長(zhǎng)度
  • H:硬跳過(guò)，切到reads中部分序列
  • N:跳過(guò)參考序列（常見(jiàn)于RNA-SEQ）

Column 7/8/9 RNEXT, PNEXT, TLEN(僅PE的數(shù)據(jù)才有)

• RNEXT:PE測(cè)序的配對(duì)read所比對(duì)到的染色體掸犬，=代表兩個(gè)比對(duì)到相同的染色體
• PNEXT:配對(duì)的read所比對(duì)到的位置。
• TLEN: 插入片段的長(zhǎng)度绪爸，可根據(jù)POS和CIGAR信息得出

Column 10/11 SEQ 和 QUAL:query的序列以及fastq序列所對(duì)應(yīng)的質(zhì)量值

Column 12 and on

此列開(kāi)始湾碎，不同測(cè)比對(duì)軟件產(chǎn)出的數(shù)據(jù)會(huì)產(chǎn)生一定的差異〉旎酰基本格式為TAG:TYPE:VALUE

SAM/BAM file Manipulation

• Select or Filter data from BAM files

### 查看特定的FLAG值
samtools flags 4 ### 0x4 unmap

### 查看所有含有4的比對(duì), -c 計(jì)數(shù)
samtools view -f 4 bwa.bam
samtools view -f 4 -c bwa.bam

### 過(guò)濾所有含有4的比對(duì)介褥，反向匹配
samtools view -F 4 bwa.bam

### -f -F參數(shù)的輸出到新文件
samtools view -b -F 4 bwa.bam > bwa.aligned.bam
samtools index bwa.aligned.bam

• BAMw文件概況統(tǒng)計(jì)

samtools flagstat bwa.bam
samtools idxstats bwa.bam
bamtools stats -in bwa.bam

### 獲取最佳比對(duì)的序列
samtools flags 2308
# unmap, secondary, supplement

### 此為proper_pair的最佳結(jié)果，其對(duì)于多次比對(duì)的reads只保留一次递惋。
samtools view -f 2 -F 2308 bwa.bam

• 對(duì)于flags信息的運(yùn)用

### 反鏈比對(duì)上的
samtools view -f 16 -c bwa.bam
# 6347

### 反鏈比對(duì)上的柔滔，且是非 未必對(duì)上的
samtools view -f 16 -F 4 -c bwa.bam

### 比對(duì)到正鏈上，且是非 未比對(duì)上的
samtools view -c -F 20 bwa.bam

• 過(guò)濾低質(zhì)量比對(duì), 5列MAPQ
  • bwa比對(duì)的MAPQ值0代表未多次比對(duì)
  • q參數(shù)表示保留比對(duì)質(zhì)量值大于等于此q值

samtools view -c -q 1 bwa.bam

###比對(duì)質(zhì)量大于1萍虽，且比對(duì)到正鏈上
samtools view -q 1 -F 4 -F 16 -c bwa.bam

• Secondary alignment 二次比對(duì):序列是多次比對(duì)睛廊，其中一個(gè)最好的比對(duì)為PRIMARY align，其余的都是二次比對(duì)杉编，F(xiàn)LAG值256

samtools flags SECONDARY
# 0x100 256

samtools view -c -F 4 -f 256 bwa.bam
# 0

• 獲得PRIMARY or Representative 比對(duì)
  • only have one primary alignment and other secondary and supplemental alignments.

samtools flags SUPPLEMENT,SECONDARY
# 0x900 2304

samtools view -c -F 4 -F 2304 bwa.bam

處理分析 SAM files

• 添加@RG分組信息

TAG='@RG\tID:xyz\tSM:Ebola\tLB:patient_100'

## 在比對(duì)過(guò)程中添加@RG信息tags
bwa mem -R $TAG $REF $R1 $R2 |samtools sort >bwa.bam
samtools index bwa.bam

### 改變分組@RG信息
NEWTAG='@RG\tID:abc\tSM:Ebola\tLB:patient_101'
samtools addreplacerg -m overwrite_all -r $NEWTAG bwa.bam -O BAM > newbam.bam

• 合并BAM文件

## bwa 比對(duì)信息
BWA_TAG='@RG\tID:bwa\tSM:bwa\tLB:bwa'
samtools addreplacerg -m overwrite_all -r $BWA_TAG bwa.bam -O BAM > newbwa.bam

## bowtie2比對(duì)信息
BOWTIE_TAG='@RG\tID:bowtie\tSM:bowtie\tLB:bowtie'
samtools addreplacerg -m overwrite_all -r $BOWTIE_TAG bowtie.bam -O BAM > newbowtie.bam

## 合并成all_merge.bam
samtools merge all_merge.bam newbwa.bam newbowtie.bam 

### 去合并之后分組信息為bwa的bam,bowtie.  **參數(shù)-l **
samtools view -c -l bwa all_merge.bam
samtools view -c -l bowtie all_merge.bam

• 指定比對(duì)到ref上的某一位置上查看信息

samtools view in.bam chr22:16050103-16050110

• 快速查看比對(duì)情況
  • bam文件比較大超全，可以截取感興趣的區(qū)段成小的bam文件再查看
  • 直接samtools tview --reference hg38.fa in.bam

###截成小的bam文件
samtools view -h bwa.bam AF086833.2:1000-1100 |samtools view -Sb - >small.bam

### 直接tview 查看
samtools tview --reference ref/ebola.fa bwa.bam

參考文章：

1. Biostar_handbook
2. 解螺旋的礦工——從零開(kāi)始完整學(xué)習(xí)全基因組測(cè)序數(shù)據(jù)分析：第5節(jié) 理解并操作BAM文件
3. 徐州更——SAM格式及其相關(guān)工具