176-再次理解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

来自:https://www.samformat.info/sam-format-flag

  • 只有一条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 
Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related