生物信息学导论-北大-新一代测序NGS：重测序得回帖和变异鉴定

ref: https://www.coursera.org/learn/sheng-wu-xin-xi-xue/home

本文主要来自本课的讲义。

---

read：测序仪读到的一条短短的DNA片段

fastq格式：序列名称+序列+质量信息

质量：$p$是base callig出错的概率，质量分数Q的计算公式为：$Q=-10\ast lo{g}_{10}(p)$

p	Q
0.1	10
0.01	20
0.001	30
0.0001	40

Q=0-40用下面的符号标识，其中10是+，20是5，30是?，40是I

!”#$%&’()*+,-./0123456789:;≤>?@ABCDEFGHI

Paired-End Reads：从两侧同时读然后拼接

Mapping Quality

假设有一个read x，长度为$l$，参考序列z的长度为L，u是x比对上z的位置，那么事实如此的概率是：

$$p(z|x,u)=\prod _{mismatch}p(Z_i)$$

使用log把连乘转换成连加：（Q是前面的质量分数）

$$SQ(u)=log(p(z|x,u))=\sum _{mismatch}p(z_i)=\sum _{mismatch}Q(z_i)$$

如果假设模型是uniform NULL model（read在任何位置出现的概率相同），那么这个read比对到给定位置u的误差可以这样计算：

$$E(u)=\frac{SQ(u)}{{\sum }_{i}SQ(i)}$$

举例：假设read是ACGT，质量值为30 30 25 20，参考序列是ACGTACGGA
表中*表示空格。

ACGTACGGA	SQ(u)	E(u)
ACGT	0+0+0+0=0	0/415
*ACGT	30+30+25+20=105	105/415
** ACGT	同上	同上
***ACGT	同上	同上
****ACGT	0+0+0+20=20	20/415
*****ACGT	30+30+20=80	80/415

Genetic Variants

• SNV（Single Nucleotide Variant）
  • SNP（Substitution）
  • Indel（insertion/deletion）
• SV（Structual Variation）
  • Large-scale insertion/deletion
  • Inversion
  • Translocation
  • CNV（Copy Number Variation）

SNP calling：找出哪些位置是多态性的，或者与参考序列有至少一个碱基不同

Genotype calling：为每个个体确定基因型，而且只在已经发现SNP或者variant的地方做

做genotype的方法：

1. counting
   • 数high-confident, non-reference的等位基因（比如要求Q≥20）
   • 频率小于20%或大于80%的是纯合的，其他是杂合的
   • 对深度测序区域（deeply sequenced regions, DSR，比如depth>25x）效果好
   • 但是对覆盖度低的区域的、杂合型效果不好
   • 无法度量可信度
2. 基于概率的简单模型
   • 假设一个二倍体基因组，那么等位基因是A和a，可能的基因型是AA，aa和Aa
   • 假设A的数量为k，a的数量为n-k
   • 那么$P(D|AA)={\prod }_{n-k}P(x_i)$，即有n-k个测序错误的概率（a有n-k个，如果全错就是全A），$P(D|aa)={\prod }_{k}P(x_i)$，$P(D|Aa)=1-P(D|AA)-P(D|aa)$
   • 如果能估算出先验概率P(AA)、P(aa)和P(Aa)，那么就可以计算出后验概率P(AA|D)等

做SNP calling和genotyping的模型：

• MAQ
• samtools
• GATK
• SNVMix

Genome Analy ToolKit GATK

• ngs data processing
• variant discovery & genotyping
• integrative analysis