美文网首页生物信息数据科学
57.《Bioinformatics Data Skills》之

57.《Bioinformatics Data Skills》之

作者: DataScience | 来源:发表于2021-08-17 21:31 被阅读0次

    碱基

    FASTA/FASTQ文件中最常见的碱基为大写的A(腺嘌呤),T(胸腺嘧啶)C(胞嘧啶),G(鸟嘌呤),而soft-mask的参考基因组中会以小写字母代表串联重复或者低复杂性序列(由RepeatMasker或Tandem Repeats Finder生成),soft-mask基因组也可以通过将这些小写字母改写为代表任意碱基的N来hard-mask。

    除了N,还有其它字母可以代表两种或者两种以上可能的碱基。国际纯粹与应用化学联合会(IUPAC)给出了所有代表确定或者不确定性碱基的字母定义(见表1):

    IUPAC nucleotide code Base
    A Adenine
    C Cytosine
    G Guanine
    T (or U) Thymine (or Uracil)
    R A or G
    Y C or T
    S G or C
    W A or T
    K G or T
    M A or C
    B C or G or T
    D A or G or T
    H A or C or T
    V A or C or G
    N any base
    . or - gap

    值得注意的是,不同的软件对代表模糊碱基的字母处理方式可能不同,例如BWA会随机取一种可能性,同时生成随机种子以确保index生成的结果一致。

    碱基质量得分

    在FASTQ文件中有一行信息专门代表碱基质量得分(如下第四行):

    @AZ1:233:B390NACCC:2:1203:7689:2153
    GTTGTTCTTGATGAGCCATGAGGAAGGCATGCCAAATTAAAATACTGGTGCGAATTTAAT
    +
    CCFFFFHHHHHJJJJJEIFJIJIJJJIJIJJJJCDGHIIIGIGIJIJIIIIJIJJIJIIH
    

    这个得分使用ASCII码进行编码,本质上ASCII码代表127个数字,不同的测序技术可能会采用不同的映射关系(如表2):

    表2 不同测序技术采用的打分类型

    现在主流采用sanger的计算方式,这里演示一下其转换过程:

    虽然ASCII码有127个,只取能够打印显示的93个(即从33-126,忽略代表空格的32)。每种程序语言都会提供ASCII码与数字之间的转换,python使用的两个函数分别是ordchr。使用ord可以将例子字符串转为数字:

    >>> quality = "CCFFFFHHHHHJJJJJEIFJIJIJJJIJIJJJJCDGHIIIGIGIJIJIIIIJIJJIJIIH"
    >>> [ord(b) for b in quality]
    [67, 67, 70, 70, 70, 70, 72, 72, 72, 72, 72, 74, 74, 74, 74, 74, 69, 73, 70, 74, 73, 74, 73, 74, 74, 74, 73, 74, 73, 74, 74, 74, 74, 67, 68, 71, 72, 73, 73, 73, 71, 73, 71, 73, 74, 73, 74, 73, 73, 73, 73, 74, 73, 74, 74, 73, 74, 73, 73, 72]
    

    根据表2,我们知道质量得分是从33开始的,所以我们需要减去偏移量33得到PHRED质量得分:

    >>> phred = [ord(b) - 33 for b in quality]
    >>> phred
    [34, 34, 37, 37, 37, 37, 39, 39, 39, 39, 39, 41, 41, 41, 41, 41, 36, 40, 37, 41, 40, 41, 40, 41, 41, 41, 40, 41, 40, 41, 41, 41, 41, 34, 35, 38, 39, 40, 40, 40, 38, 40, 38, 40, 41, 40, 41, 40, 40, 40, 40, 41, 40, 41, 41, 40, 41, 40, 40, 39]
    

    PHRED质量得分最早由Phil Green创造用于荧光质量打分,这个得分可以按如下公式转换为碱基出错的可能性P值:
    P = 10^{-Q/10}
    每个碱基的出错可能性为:

    >>> [10**(-q/10) for q in phred]
    [0.00039810717055349735, 0.00039810717055349735, 0.00019952623149688788, 0.00019952623149688788, 0.00019952623149688788, 0.00019952623149688788, 0.00012589254117941674, 0.00012589254117941674, 0.00012589254117941674, 0.00012589254117941674, 0.00012589254117941674, 7.943282347242822e-05, 7.943282347242822e-05, 7.943282347242822e-05, 7.943282347242822e-05, 7.943282347242822e-05, 0.00025118864315095795, 0.0001, 0.00019952623149688788, 7.943282347242822e-05, 0.0001, 7.943282347242822e-05, 0.0001, 7.943282347242822e-05, 7.943282347242822e-05, 7.943282347242822e-05, 0.0001, 7.943282347242822e-05, 0.0001, 7.943282347242822e-05, 7.943282347242822e-05, 7.943282347242822e-05, 7.943282347242822e-05, 0.00039810717055349735, 0.00031622776601683794, 0.00015848931924611142, 0.00012589254117941674, 0.0001, 0.0001, 0.0001, 0.00015848931924611142, 0.0001, 0.00015848931924611142, 0.0001, 7.943282347242822e-05, 0.0001, 7.943282347242822e-05, 0.0001, 0.0001, 0.0001, 0.0001, 7.943282347242822e-05, 0.0001, 7.943282347242822e-05, 7.943282347242822e-05, 0.0001, 7.943282347242822e-05, 0.0001, 0.0001, 0.00012589254117941674]
    

    Illumina标准采用类似的流程,除了偏移量为64。Solexa标准使用的转换公式为Q=(10^{P/10}+1)^{-1},会稍微复杂一些,感兴趣的话可以自行了解一下。

    相关文章

      网友评论

        本文标题:57.《Bioinformatics Data Skills》之

        本文链接:https://www.haomeiwen.com/subject/vbzibltx.html