Charpter_14 Sequence Alignment Maps(SAM)
SAM文件是由Tab(\t)分割，面向‘行’的文本格式文件，主要包括两部分：
- Header部分，包括
@HD,@SQ,@PG等开头的基本信息。 - Alignment部分，包括reads比对的所有信息。
sam是最常见的存放mapping数据的格式，各种组学分析只要存在mapping的步骤都会产生SAM格式文件用于后续进一步的下游分析。
比对+samtools sort一个快速示例
REF=ebola.fa
R1=SRR1972739_1.fq
R2=SRR1972739_2.fq
### bwa mem $REF $R1 $R2 > bwa.sam
### bowtie 比对输出直接sort排序产生bam文件
bowtie2 -x $REF -1 $R1 -2 $R2 |samtools sort -@ 10 >bowtie2_output.sorted.bam
SAM header
header部分由@起始，包括有：
-
HD：表示排序情况，unsorted或者已经按coordinate/names排过序 -
SQ：比对的参考序列的染色体信息，姓名/多长。 -
PG：运行命令的具体参数信息。 -
RG：每个read所在的group的信息。通常包含测序平台/测序文库/样本的ID等信息(header里最为重要的一个信息).ID(Read group identifier), LB(library), SM(Sample)
@RG信息：如果原来样本的测序深度比较深，一般会按照不同的lane分开比对，最后再合并到一起。那么这个时候会有多个RG，里面记录了不同的lane，甚至测序文库的信息，唯一不变的是SM的sample信息，这样合并后才能正确处理。
比对信息(Alignment部分)
image
image
Column 1. QNAME(query name)：fastq文件里的read ID
Column 2. FLAG信息
- 记录了许多有关read比对情况的信息，转换成二进制码后发现其包括12位，十进制取值范围就是0~2048，2的12次方。
| FLAG | 二进制 | 具体含义 |
|---|---|---|
| 1 | 000000000001 | 代表是PE双端测序，0代表SE单端测序 |
| 2 | 000000000010 | 代表read和参考序列完全匹配。如果是PE测序还代表没有插入确实 |
| 4 | 000000000100 | 该read未比对到参考序列上 |
| 8 | 000000001000 | PE测序的另一端R2未比对到参考序列上 |
| 16 | 000000010000 | 反向互补后比对到参考序列上，read比对到了反向互补链上 |
| 32 | 000000100000 | PE测序的另一条reads(R2)比对到了反向互补链上 |
| 64 | 000001000000 | PE测序read1 |
| 128 | 000010000000 | PE 测序的read2 |
| 256 | 000100000000 | 代表二次比对，该read在基因组上比对了多个位置，只有一个是首要比对位置，其它都是次要的，该位置是次要的，通常要过滤掉，但有些场景中是有用的信息 |
| 512 | 001000000000 | 代表此序列低于测序平台过滤阈值，QC失败无法过滤(不常用) |
| 1024 | 010000000000 | PCR重复序列或光学重复，可被Picard标记过滤掉（来自测序文库构建过程） |
| 2048 | 100000000000 | 该read可能存在嵌合，这个比对的部分只是来自其中的一部分序列(supplement alignment) |
- 对于未知的FLAG，可以
samtools flags 141 - 对于FLAG=77时，需要转换为二进制序列：77=000001001101=1+4+8+64即表示：PE read，比对不上参考序列，另一端也未比对上参考序列，它是PE测序read1
-
常见的FLAGS：
- 其中一条reads没有map上： 73, 133, 89 121, 165, 181, 101, 117, 153, 185, 59, 137
- 两条reads都未map上：77，141
- 比对上了，方向也对，也在插入大小(insert size)内： 99, 147, 83, 163
- 比对上了，也在插入大小(insert size)内， 但是反向不对:67, 131, 115, 179
- 单一配对，就是插入大小(insert size)不对: 81, 161, 97, 145, 65, 129, 113, 177
Column 3/4. RNAME(reference name)/POS(position)：
- RNAME:表示参考序列的名字
- POS：比对上的坐标位置，注意不管是正向还是反向，只记录比对上的最左边的坐标位置，坐标是以1为起始
Column 5/6. MAPQ(mapping quality)和CIGAR(compact idiosyncratic gapped alignment representation)
- MAPQ: 比对质量值，比对到参考序列上此位置的可靠程度，转化为-10logP。40/10=4 ——> 10^-4=1/10000，即比错的概率小于0.0001。
- CIGAR：reads比对到参考序列的具体细节情况：
-
M:(完全匹配或SNP单碱基错配) -
I:序列插入，包含潜在的insertion变异 -
D:序列删除，包含潜在的deletion变异 -
S:软跳过,跳过reads中部分序列，不改变read长度 -
H:硬跳过，切到reads中部分序列 -
N:跳过参考序列（常见于RNA-SEQ）
-
Column 7/8/9 RNEXT, PNEXT, TLEN(仅PE的数据才有)
- RNEXT:PE测序的配对read所比对到的染色体，
=代表两个比对到相同的染色体 - PNEXT:配对的read所比对到的位置。
- TLEN: 插入片段的长度，可根据POS和CIGAR信息得出
Column 10/11 SEQ 和 QUAL:query的序列以及fastq序列所对应的质量值
Column 12 and on
此列开始，不同测比对软件产出的数据会产生一定的差异。基本格式为TAG:TYPE:VALUE
SAM/BAM file Manipulation
- Select or Filter data from BAM files
### 查看特定的FLAG值
samtools flags 4 ### 0x4 unmap
### 查看所有含有4的比对, -c 计数
samtools view -f 4 bwa.bam
samtools view -f 4 -c bwa.bam
### 过滤所有含有4的比对，反向匹配
samtools view -F 4 bwa.bam
### -f -F参数的输出到新文件
samtools view -b -F 4 bwa.bam > bwa.aligned.bam
samtools index bwa.aligned.bam
- BAMw文件概况统计
samtools flagstat bwa.bam
samtools idxstats bwa.bam
bamtools stats -in bwa.bam
### 获取最佳比对的序列
samtools flags 2308
# unmap, secondary, supplement
### 此为proper_pair的最佳结果，其对于多次比对的reads只保留一次。
samtools view -f 2 -F 2308 bwa.bam
- 对于flags信息的运用
### 反链比对上的
samtools view -f 16 -c bwa.bam
# 6347
### 反链比对上的，且是非 未必对上的
samtools view -f 16 -F 4 -c bwa.bam
### 比对到正链上，且是非 未比对上的
samtools view -c -F 20 bwa.bam
-
过滤低质量比对, 5列MAPQ
- bwa比对的MAPQ值0代表未多次比对
-
q参数表示保留比对质量值大于等于此q值
samtools view -c -q 1 bwa.bam
###比对质量大于1，且比对到正链上
samtools view -q 1 -F 4 -F 16 -c bwa.bam
- Secondary alignment 二次比对:序列是多次比对，其中一个最好的比对为PRIMARY align，其余的都是二次比对，FLAG值256
samtools flags SECONDARY
# 0x100 256
samtools view -c -F 4 -f 256 bwa.bam
# 0
-
获得PRIMARY or Representative 比对
- 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'
## 在比对过程中添加@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 比对信息
BWA_TAG='@RG\tID:bwa\tSM:bwa\tLB:bwa'
samtools addreplacerg -m overwrite_all -r $BWA_TAG bwa.bam -O BAM > newbwa.bam
## bowtie2比对信息
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. **参数-l **
samtools view -c -l bwa all_merge.bam
samtools view -c -l bowtie all_merge.bam
- 指定比对到ref上的某一位置上查看信息
samtools view in.bam chr22:16050103-16050110
- 快速查看比对情况
- bam文件比较大，可以截取感兴趣的区段成小的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
参考文章：
- Biostar_handbook
- 解螺旋的矿工——从零开始完整学习全基因组测序数据分析：第5节 理解并操作BAM文件
- 徐州更——SAM格式及其相关工具
