- 57.《Bioinformatics Data Skills》之
- 28.《Bioinformatics-Data-Skills》之
- 18.《Bioinformatics-Data-Skills》之
- 19.《Bioinformatics-Data-Skills》之
- 【shell笔记>生信|专项】生信数据处理技能手札(3)：
- Bioinformatics Data Skills
- 17.《Bioinformatics-Data-Skills》之
- 25.《Bioinformatics-Data-Skills》之
- 25.《Bioinformatics-Data-Skills》之
- 23.《Bioinformatics-Data-Skills》之
碱基
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使用的两个函数分别是ord与chr。使用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值：
每个碱基的出错可能性为：
>>> [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标准使用的转换公式为，会稍微复杂一些，感兴趣的话可以自行了解一下。
网友评论