zoukankan      html  css  js  c++  java
  • samtools flags 的含义

    对于双端比对的数据,生成的BAM文件中,R1端序列和R2端序列的标识符是一样的,之前一直不知道如何根据bam文件区分哪条序列是R1端,哪条序列是R2端,昨天仔细研究了一下,原来代表R1端和R2端的信息都存储在flag中,即bam文件的第二列;

    在bam文件格式中定义了各种flag代表的意思

    /*! @abstract the read is paired in sequencing, no matter whether it is mapped in a pair */
    #define BAM_FPAIRED        1
    /*! @abstract the read is mapped in a proper pair */
    #define BAM_FPROPER_PAIR   2
    /*! @abstract the read itself is unmapped; conflictive with BAM_FPROPER_PAIR */
    #define BAM_FUNMAP         4
    /*! @abstract the mate is unmapped */
    #define BAM_FMUNMAP        8
    /*! @abstract the read is mapped to the reverse strand */
    #define BAM_FREVERSE      16
    /*! @abstract the mate is mapped to the reverse strand */
    #define BAM_FMREVERSE     32
    /*! @abstract this is read1 */
    #define BAM_FREAD1        64
    /*! @abstract this is read2 */
    #define BAM_FREAD2       128
    /*! @abstract not primary alignment */
    #define BAM_FSECONDARY   256
    /*! @abstract QC failure */
    #define BAM_FQCFAIL      512
    /*! @abstract optical or PCR duplicate */
    #define BAM_FDUP        1024
    /*! @abstract supplementary alignment */
    #define BAM_FSUPPLEMENTARY 2048

    1 : 代表这个序列采用的是PE双端测序

    2: 代表这个序列和参考序列完全匹配,没有插入缺失

    4: 代表这个序列没有mapping到参考序列上

    8: 代表这个序列的另一端序列没有比对到参考序列上,比如这条序列是R1,它对应的R2端序列没有比对到参考序列上

    16:代表这个序列比对到参考序列的负链上

    32 :代表这个序列对应的另一端序列比对到参考序列的负链上

    64 : 代表这个序列是R1端序列, read1;

    128 : 代表这个序列是R2端序列,read2;

    256: 代表这个序列不是主要的比对,一条序列可能比对到参考序列的多个位置,只有一个是首要的比对位置,其他都是次要的

    512: 代表这个序列在QC时失败了,被过滤不掉了(# 这个标签不常用)

    1024: 代表这个序列是PCR重复序列(#这个标签不常用)

    2048: 代表这个序列是补充的比对(#这个标签具体什么意思,没搞清楚,但是不常用)

    上面的这几个标签都是2的n次方,这样的数列有一个特点,就是随机挑选其中的几个,它们的和是唯一的,比如

    65 只能是1 和 64 组成,代表这个序列是双端测序,而且是read1

    所以在bam文件中的第二列,即flag列的值代表这条序列符合上述所有条件的值的和,所以根据这个flag我们可以确定这条序列究竟是read1 还是read2

    在新版的samtools中提供了一个flags 功能,可以查看flag 代表的含义,比如下面的sam文件

    NB500986:26:HGJ2VBGXX:3:13403:17250:17986       419     chr1    12013   1       131M    =       12016   131     TCTGACTTCCAGCAACTGCTGGCCTGTGCCAGGGTGCAAGCTGAG
    CACTGGAGTGGAGTTTTCCTGTGGAGAGGAGCCATGCCTAGAGTGGGATGGGCCATTGTTCATCTTCTGGCCCCTGTTGTCTGCAT  AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEAEEEEAEEEEEEEEEEEEEEEEEEEE
    EEEEEEEEEEEEEEEEEEEEEEEEE/EE/EE<EEEEEEEEEEEEEEEEEEEEEEEEEEEEEE  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:131        YT:Z:UU XS:A:+  NH:i:4  CC:Z:
    chr15   CP:i:102519027  HI:i:0
    NB500986:26:HGJ2VBGXX:3:13403:17250:17986       339     chr1    12016   1       128M    =       12013   -131    GACTTCCAGCAACTGCTGGCCTGTGCCAGGGTGCAAGCTGAGCAC
    TGGAGTGGAGTTTTCCTGTGGAGAGGAGCCATGCCTAGAGTGGGATGGGCCATTGTTCATCTTCTGGCCCCTGTTGTCTGCAT     EEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEE<AEEEEEEEEEEEAEEEAEEEEEEEEEEEAEEEEE
    EEEEEEEEEEEEEEEEEEEEEAAEEAEEEEEAEEEEEEEEEE/EEEEEEEEEEEAAAAA     AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:128        YT:Z:UU XS:A:+  NH:i:4  CC:Z:
    chr15   CP:i:102519027  HI:i:0

    一共有两个序列,他们的标识符是一样的,那怎么区分那个是R1端,那个是R2端呢?

    直接根据flag 查看

    samtools flags 419
    0x1a3   419     PAIRED,PROPER_PAIR,MREVERSE,READ2,SECONDARY
    samtools flags 339
    0x153   339     PAIRED,PROPER_PAIR,REVERSE,READ1,SECONDARY

    根据上面给结果可以看出,flag为419代表PAIREED(1), PROPER_PAIR(2),MREVRSE(32),READ2(128),SECONDARY(256),其中

    PAIRED 代表这条序列采用双端测序, 其值为1;

    PROPER_PAIR 代表这条序列完全匹配, 其值为2;

    MREVRSE 代表这条序列对应的另一端序列比对到参考序列的负链上,其值为32;

    READ2 代表这条序列是R2端序列,其值为128

    SECONDARY 代表这条序列不是primary alignment, 其值为256

    1 + 2 +32 +128 + 256 = 419

    所以419对应的序列是R2端序列, 同理,339代表的序列是R1段序列

    搞清楚了序列究竟是R1端还是R2端之后,还有一点值得注意的地方,samtools 还提供了一个bam2fq 的功能,根据bam文件抽取比对时采用的fastq序列

    在bam文件中,第10列代表的是序列,第11列代表的是序列的质量,根据第二列的flag还可以确定序列是R1还是R2,

    但是当序列本身比对到参考序列的负链时,即flag 包含16时,bam文件中记录的序列是原始序列的反向互补序列,而且质量值也是反向的,所以根据这样的序列还原时要小心一点,以上面flag为339的序列为例

    因为339包含了REVERSE,所以对应的序列应该是第10列序列的反向互补序列,碱基质量值为第11列的反向序列,

    搞清楚了这些,对于bam文件的理解又更清晰了一些。

  • 相关阅读:
    4.28综合练习
    团队项目第一阶段冲刺第六天
    4.27防盗链和代理
    梦断代码阅读笔记3
    团队项目第一阶段冲刺第五天
    4.26抓取猪⼋戒数据
    团队项目第一阶段冲刺第四天
    4.25xpath解析
    4.24aiohttp模块学习
    如何将类数组转化为数组?
  • 原文地址:https://www.cnblogs.com/xudongliang/p/5437850.html
Copyright © 2011-2022 走看看