再次理解SAM/BAM操作
刘小泽写于2020.3.14
要说生信的重难点，不断转换数据格式一定是排在前三位的
上次做了一次外显子分析，再一次认识了sam/bam格式，整理一些内容和你分享
以前写过一篇：处理SAM、BAM你需要Samtools ，其中介绍了bam/sam格式
首先是格式问题
-
PE reads在bam中默认是两条记录，除非设置参数将多比对结果也输出
-
一般是11列，其中1、10、11合起来就是fq格式
-
2列：二进制flag（查看：https://www.samformat.info/sam-format-flag）
-
3、7列：判断某条reads是否比对成功到了基因组的染色体，左右两条reads是否比对到同一条染色体 [左右端测序数据比对到不同染色体的情况，比较有意义，可能是融合基因，也可能是基因之间本来就相似性很大]
-
4列：染色体位置
-
5列：比对结果的质量值MAPQ（http://biofinysics.blogspot.com/2014/05/how-does-bowtie2-assign-mapq-scores.html），结果越高越可信。为0表示在参考基因组有多种定位的可能性
-
6列：CIGAR （其中soft-clipping是指一条reads未匹配上当前基因组位置的部分，如果有多个reads在这种情况并且这些reads的soft-clipping碱基都能够比对在基因组另一位置，那么就可能存在SV）
-
9列是我们建库的时候打断的片段长度（PE 150测序，一般是350左右）
-
11列以后附带的（各种tag）：
- RG代表着你的sam文件比对来自于哪个样本的fastq结果；
- NM是编辑距离（reads如果想转变成参考基因组，需要改变多少个碱基），编辑距离是0才说明150bp长度的序列跟参考基因组一致
- MD是序列跟参考基因组不同在哪里
关于BAM Flag
来自：https://davetang.org/muse/2014/03/06/understanding-bam-flags/
获得flag信息
#take the second column of the BAM file
#and output all the unique entries
#the second column in the BAM flag column
samtools view test.bam | cut -f2 | sort | uniq -c
- 只有一条reads没有比对上: 73, 133, 89, 121, 165, 181, 101, 117, 153, 185, 69, 137
- 两条reads都没有比对上: 77、141
- 比对上了，方向也对，并且在插入片段大小范围（insert size）内: 99, 147, 83, 163
- 比对上了，也在插入片段大小范围内， 但是方向不对:67, 131, 115, 179
- 唯一配对，就是插入片段大小范围不对: 81, 161, 97, 145, 65, 129, 113, 177
将bam拆分成小bam
方法一 使用工具
bamtools split \
-in lung-1.sort.bam \
-reference # 指定按照reference来分离bam文件
# 另外，可以指定-mapped来分离比对成功与否的bam文件
# 可以指定 -tag RG 来把这个bam文件按照原来的测序上样品的lane给分离开
方法二：使用脚本
for chrom in `seq 1 22` X Y MT do (samtools view -bh lung-1.sort.bam\
chr${chrom} | samtools sort - chr${chrom} && samtools index chr${chrom}.bam) done
提取未比对成功的bam
注意：未比对成功的测序数据可以分成3类（仅reads1，仅reads2和两端reads都没有比对成功）
方法一：直接使用samtools
samtools view -f4 sample.bam > sample.unmapped.sam
# 小写的f是提取，大写的F是过滤
方法二：根据bam格式
第3列表示read1，第7列表示read2，因此提取它们就是：要么第3列是*，要么第7列是*，或者都是*
samtools view sample.bam |perl -alne '{print if $F[2] eq "*" or $F[6] eq "*" }' > sample.unmapped.sam
方法三：bamtools
bamtools -split -in my.bam -mapped
bam过滤低质量
仅仅过滤MAPQ=0
samtools view -q 1 tmp.bam
# 过滤掉MAPQ值低于1的情况
去除掉multiple mapping、PCR duplication及低质量比对
samtools view -h -F4 -q 5 test.bam |samtools view -bS |samtools rmdup -o test.filter.rmdup.bam
samtools index test.filter.rmdup.bam
bam去PCR重复
单端数据
samtools rmdup -s tmp.sorted.bam tmp.rmdup.bam
# 比对flag情况只有0,4,16，只需要去掉比对起始、终止坐标一致的reads
双端数据
最好用picard的MarkDuplicates
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" MarkDuplicates \
-I test.bam \
--REMOVE_DUPLICATES=true \
-O test.marked.bam \
-M test.metrics
bam转为fq
# https://bedtools.readthedocs.io/en/latest/content/tools/bamtofastq.html
# 先将bam安装reads名称排序（-n），保证PEreads相邻
samtools sort -n aln.bam aln.qsort
bedtools -i test.bam -fq tmp1.fq -fq2 tmp2.fq
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
