比较不同样本同一区域测序深度的R实现
一、测序深度概率分布模型
只有理解测序深度分布模型，才能更好地理解测序数据，理解测序深度相关的分析，比如拷贝数变异检测、差异表达分析等。下文基于参考资料，加上自己的理解写成，理解的也不是很正确，欢迎留言指教。
1-多项分布
ki为某个region上比对的reads总数，一次测序好比从一个瓮中抽一个彩球，一个瓮就是一个pool，抽某个颜色的球的概率pi取决于彩球在瓮中的比例，如果抽样次数即测序总reads数N远远小于DNA模板数（即瓮中彩球总数），那么每种球抽到的个数（每个区域上的片段总数ki的向量)符合多项分布，每种球被抽到的概率pi是不会随着抽样次数N而改变的。但是随着我们测序通量的加大，测序的reads数有时候甚至能够超过模板数，因此这个模型可能不是太适用了。
2-二项分布和泊松分布
如果单独地考虑一个片段被测序的次数（一个球是否被抽到，不放回抽样），是否被测序（抽到没抽到）只有是否2个选项，测序彼此之间相互独立（每次抽相互独立，球与球之间不会互相干扰），区域i上测序reads数ki即符合二项分布，试验次数N即测序总reads数，成功概率取决于region i的片段占中片段的比例和测序的容易程度，pi会随着N而变化。N变大，pi变小，二项分布会收敛于泊松分布。实际的测序比较符合后者，因为现在的NGS里测序reads数都很大。
3-过度离散的泊松分布（复合泊松分布）
上述的理想实验描述了从一个非常大的DNA片段总体中重复抽样的过程，但是，重复抽样虽然可以作为技术重复的假设，却没办法作为生物学重复的假设。一个实验里的生物学重复就意味着生成了一个新的DNA片段总体，在这个新的样本总体里，基因组上不同区域上被抽样的概率p不是完全相同。因此虽然泊松近似的二项分布依然可以来描述一个单独的抽样，因为N依然很大，而pi很小，但是样本j在特征（比如区域）i上读断数的期望值E（Kij）会随着样本j而变化。测序深度符合泊松分布，测序深度的均值也符合泊松分布，所以就是复合泊松分布。而且因为方差大于均值，所以这个泊松分布是过度离散的。就是说不同样本同一个区域之间的reads数方差大于他们的均值。
为了构建这样一个模型，模型需要有以下的特征：
1）模型需要能够支持0~∞即非负实数，因为它代表的是测序的reads数，而这个数是非负的；
2）需要至少有2个参数来指定参数的均值和方差，当均值的方差为0时，就是说均值不变的时候，就是泊松分布。
4-泊松-伽马混合分布/负二项分布
伽马分布能够满足上述模型的两个特征。泊松-伽马混合分布，即泊松分布含有一个符合伽马分布的均值参数，通常被称为负二项分布。“负二项分布”说明了这个分布和数伯努利试验成功次数一样，即成功为止需要的实验次数的分布。随机变量的密度K~NB（μ，α），其中均值参数μ>0，分散参数α>0，定义：
5-对数正态分布
另一种泊松的平均参数可能的分布是对数正态分布，即泊松对数正态分布。但是泊松对数正态不像负二项分布一样存在分布边界，随机样本可以通过正态随机变量来快速地建立，Z~N(m, σ2)，并且数量K~Pois（e^Z）。注意，E(e^Z)=e^(m+σ^2/2)，即均值e^Z不仅仅指正态随机变量Z的期望均值。但是负二项分布和泊松对数正态分布在α和σ2很小的时候非常近似，当这两个数变大的时候模型就不一样了。当分散系数提高，负二项分布在概率上的权重提高，P（k=0）,可是泊松对数分布给右边尾巴的权重更重。
二 R 实现
学习了上述材料之后，尝试自己构建一个模型，用R语言实现。
虽然上述公式很复杂，但是我看了好几个例子，里面都是用非常直观的方法算出来的，甚至不需要搞清楚公式的含义——首先确定深度的期望公式，然后把每部分的对数似然概率相加，然后用最大似然法求出参数，再把参数带入模型即可。
下面是我的测试数据：
数据集如下，37个样本某个区域上的平均测序深度，深度已经用样本总碱基数进行了均一化，这37个测序数据来自三个批次（理论上应该是3个pool，但是我没有去区分他们），用的是相同的捕获探针。其中，已知S21样本纯合缺失，S4和S8是S21的父母，预期是杂合缺失，其他样本为同一批上机的其他数据，是否缺失是未知的。但是该区域的缺失已知只在患有相关疾病的人群里出现，是致病的缺失，预期应该是少数。
hotspot=data.frame(
id=c("S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9", "S10", "S11", "S12", "S13", "S14", "S15", "S16", "S17", "S18", "S19", "S20", "S21", "S22", "S23", "S24", "S25", "S26", "S27", "S28", "S29", "S30", "S31", "S32", "S33", "S34", "S35", "S36", "S37"),
depth=c(89,75,89,44,92,75,89,43,74,80,66,118,73,107,99,90,39,82,75,62,2,79,94,93,51,75,57,94,73,57,116,104,85,103,84,111,85))
先计算一下depth的均值和方差，可以看到均值小于方差。
> mean(hotspot$depth)
[1] 79.02703
> var(hotspot$depth)
[1] 557.3048
depth的分布如下
hist(hotspot$depth)
首先，我们要先写出深度的公式：深度可能的因素有测序总数据量，区域长度，区域GC含量等，如果比较的是不同的区域，需要考虑一下这些因素。然后写出深度期望的公式，比如：
μ ~ θ1 * len + θ2*GC_content + θ3*total_bases + θ0。然后我们考虑一下每个参数都符合什么概率分布，depth满足depth~P(μ）或者depth~NB (μ, α），theta0满足θ0~P(λ)，len、GC、total base都是常量。根据概率分布的可加性，然后把每个参数的对数似然概率加起来。然后写出对数似然函数，再用optim()或者nlm()函数求出最优参数。
我们要使用的数据集是不同样本同一个区域的深度，depth的公式为 μ = depth上划线， μ是depth的期望。负二项分布有两个参数，mu和分散度size，μ = depth'*size。
1- 泊松分布的似然函数
negpois.LL<-function(p){
mu <- p[1]
LL <- sum(dpois(hotspot$depth, lambda = mu, log=T))
-LL
}
用nlm()函数进行最大似然法参数估计，2是初始值，可能还需要规定上下限，否则会警告或者收敛失败。
out.pois1 <- nlm(negpois.LL, 2)
然后把估计的参数带入概率分布模型，求得概率：
indivi.p0 <- dpois(hotspot$depth, lambda=out.pois1$estimate[1])
然后我们按照概率大小进行排序：
hotspot[order(indivi.p0),]
> hotspot[order(indivi.p0),]
id depth
21 S21 2
17 S17 39
8 S8 43
4 S4 44
12 S12 118
31 S31 116
36 S36 111
25 S25 51
14 S14 107
32 S32 104
34 S34 103
27 S27 57
给概率设置一个阈值，能够划定缺失界限。当然也可以用更加复杂的算法，比如HMM来找可能的缺失或者扩增。但是我还不会自己码这种代码OTL。
2- 负二项分布的似然函数
negNB.LL <- function(p){
mu <- p[1]
theta <- p[2]
LL <- sum(dnbinom(hotspot$depth, mu=mu, size=theta, log=T))
-LL
}
进行同样的操作
out.NB1 <- nlm(negNB.LL, c(2,1))
indivi.p0 <- dnbinom(hotspot$depth, mu=out.NB1$estimate[1], size=out.NB1$estimate[2])
hostpot[order(indivi.p0),]
> hotspot[order(indivi.p0),]
id depth
21 S21 2
12 S12 118
31 S31 116
36 S36 111
14 S14 107
17 S17 39
32 S32 104
34 S34 103
8 S8 43
15 S15 99
4 S4 44
23 S23 94
28 S28 94
24 S24 93
发现泊松比NB要好，可能是这个情景下泊松对于均值左侧的更加敏感？
不过后续如果用诸如HMM或者贝叶斯之类的算法来找缺失和重复的话，应该用这两种模型都是没有问题的。
但是最近在使用的时候，发现这个算法在深度很低的时候，主要是检测杂合的时候，灵敏度降低。所以接下来又学习了一些东西，如果解决了这个问题，将写一篇新的博文。
参考文献：
[1] Statistical analysis of high-throughput sequencing count data
[2] 最大似然法
[3] 测序读断的负二项分布模型
[4] 负二项分布R实现
只有理解测序深度分布模型，才能更好地理解测序数据，理解测序深度相关的分析，比如拷贝数变异检测、差异表达分析等。下文基于参考资料，加上自己的理解写成，理解的也不是很正确，欢迎留言指教。
1-多项分布
ki为某个region上比对的reads总数，一次测序好比从一个瓮中抽一个彩球，一个瓮就是一个pool，抽某个颜色的球的概率pi取决于彩球在瓮中的比例，如果抽样次数即测序总reads数N远远小于DNA模板数（即瓮中彩球总数），那么每种球抽到的个数（每个区域上的片段总数ki的向量)符合多项分布，每种球被抽到的概率pi是不会随着抽样次数N而改变的。但是随着我们测序通量的加大，测序的reads数有时候甚至能够超过模板数，因此这个模型可能不是太适用了。
2-二项分布和泊松分布
如果单独地考虑一个片段被测序的次数（一个球是否被抽到，不放回抽样），是否被测序（抽到没抽到）只有是否2个选项，测序彼此之间相互独立（每次抽相互独立，球与球之间不会互相干扰），区域i上测序reads数ki即符合二项分布，试验次数N即测序总reads数，成功概率取决于region i的片段占中片段的比例和测序的容易程度，pi会随着N而变化。N变大，pi变小，二项分布会收敛于泊松分布。实际的测序比较符合后者，因为现在的NGS里测序reads数都很大。
3-过度离散的泊松分布（复合泊松分布）
上述的理想实验描述了从一个非常大的DNA片段总体中重复抽样的过程，但是，重复抽样虽然可以作为技术重复的假设，却没办法作为生物学重复的假设。一个实验里的生物学重复就意味着生成了一个新的DNA片段总体，在这个新的样本总体里，基因组上不同区域上被抽样的概率p不是完全相同。因此虽然泊松近似的二项分布依然可以来描述一个单独的抽样，因为N依然很大，而pi很小，但是样本j在特征（比如区域）i上读断数的期望值E（Kij）会随着样本j而变化。测序深度符合泊松分布，测序深度的均值也符合泊松分布，所以就是复合泊松分布。而且因为方差大于均值，所以这个泊松分布是过度离散的。就是说不同样本同一个区域之间的reads数方差大于他们的均值。
为了构建这样一个模型，模型需要有以下的特征：
1）模型需要能够支持0~∞即非负实数，因为它代表的是测序的reads数，而这个数是非负的；
2）需要至少有2个参数来指定参数的均值和方差，当均值的方差为0时，就是说均值不变的时候，就是泊松分布。
4-泊松-伽马混合分布/负二项分布
伽马分布能够满足上述模型的两个特征。泊松-伽马混合分布，即泊松分布含有一个符合伽马分布的均值参数，通常被称为负二项分布。“负二项分布”说明了这个分布和数伯努利试验成功次数一样，即成功为止需要的实验次数的分布。随机变量的密度K~NB（μ，α），其中均值参数μ>0，分散参数α>0，定义：
在这个参数公式里，分散参数α等于失败次数的倒数，或者说等于伯努利试验的次数。均值和方差为：
负二项分布首先在生态学上被描述，即一个特定物种在不同区域或者随着时间的数量变化，这个分散参数最早就是在生态学场景下被定义的。Robinson等人认为负二项分布可以用于测序数据，包括RNA-seq。当α→0时，负二项分布等于泊松分布，满足特征（2）。
5-对数正态分布
另一种泊松的平均参数可能的分布是对数正态分布，即泊松对数正态分布。但是泊松对数正态不像负二项分布一样存在分布边界，随机样本可以通过正态随机变量来快速地建立，Z~N(m, σ2)，并且数量K~Pois（e^Z）。注意，E(e^Z)=e^(m+σ^2/2)，即均值e^Z不仅仅指正态随机变量Z的期望均值。但是负二项分布和泊松对数正态分布在α和σ2很小的时候非常近似，当这两个数变大的时候模型就不一样了。当分散系数提高，负二项分布在概率上的权重提高，P（k=0）,可是泊松对数分布给右边尾巴的权重更重。
学习了上述材料之后，尝试自己构建一个模型，用R语言实现。
虽然上述公式很复杂，但是我看了好几个例子，里面都是用非常直观的方法算出来的，甚至不需要搞清楚公式的含义——首先确定深度的期望公式，然后把每部分的对数似然概率相加，然后用最大似然法求出参数，再把参数带入模型即可。
下面是我的测试数据：
数据集如下，37个样本某个区域上的平均测序深度，深度已经用样本总碱基数进行了均一化，这37个测序数据来自三个批次（理论上应该是3个pool，但是我没有去区分他们），用的是相同的捕获探针。其中，已知S21样本纯合缺失，S4和S8是S21的父母，预期是杂合缺失，其他样本为同一批上机的其他数据，是否缺失是未知的。但是该区域的缺失已知只在患有相关疾病的人群里出现，是致病的缺失，预期应该是少数。
hotspot=data.frame(
id=c("S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8", "S9", "S10", "S11", "S12", "S13", "S14", "S15", "S16", "S17", "S18", "S19", "S20", "S21", "S22", "S23", "S24", "S25", "S26", "S27", "S28", "S29", "S30", "S31", "S32", "S33", "S34", "S35", "S36", "S37"),
depth=c(89,75,89,44,92,75,89,43,74,80,66,118,73,107,99,90,39,82,75,62,2,79,94,93,51,75,57,94,73,57,116,104,85,103,84,111,85))
先计算一下depth的均值和方差，可以看到均值小于方差。
> mean(hotspot$depth)
[1] 79.02703
> var(hotspot$depth)
[1] 557.3048
depth的分布如下
hist(hotspot$depth)
首先，我们要先写出深度的公式：深度可能的因素有测序总数据量，区域长度，区域GC含量等，如果比较的是不同的区域，需要考虑一下这些因素。然后写出深度期望的公式，比如：
μ ~ θ1 * len + θ2*GC_content + θ3*total_bases + θ0。然后我们考虑一下每个参数都符合什么概率分布，depth满足depth~P(μ）或者depth~NB (μ, α），theta0满足θ0~P(λ)，len、GC、total base都是常量。根据概率分布的可加性，然后把每个参数的对数似然概率加起来。然后写出对数似然函数，再用optim()或者nlm()函数求出最优参数。
我们要使用的数据集是不同样本同一个区域的深度，depth的公式为 μ = depth上划线， μ是depth的期望。负二项分布有两个参数，mu和分散度size，μ = depth'*size。
1- 泊松分布的似然函数
negpois.LL<-function(p){
mu <- p[1]
LL <- sum(dpois(hotspot$depth, lambda = mu, log=T))
-LL
}
用nlm()函数进行最大似然法参数估计，2是初始值，可能还需要规定上下限，否则会警告或者收敛失败。
out.pois1 <- nlm(negpois.LL, 2)
然后把估计的参数带入概率分布模型，求得概率：
indivi.p0 <- dpois(hotspot$depth, lambda=out.pois1$estimate[1])
然后我们按照概率大小进行排序：
hotspot[order(indivi.p0),]
> hotspot[order(indivi.p0),]
id depth
21 S21 2
17 S17 39
8 S8 43
4 S4 44
12 S12 118
31 S31 116
36 S36 111
25 S25 51
14 S14 107
32 S32 104
34 S34 103
27 S27 57
给概率设置一个阈值，能够划定缺失界限。当然也可以用更加复杂的算法，比如HMM来找可能的缺失或者扩增。但是我还不会自己码这种代码OTL。
2- 负二项分布的似然函数
negNB.LL <- function(p){
mu <- p[1]
theta <- p[2]
LL <- sum(dnbinom(hotspot$depth, mu=mu, size=theta, log=T))
-LL
}
进行同样的操作
out.NB1 <- nlm(negNB.LL, c(2,1))
indivi.p0 <- dnbinom(hotspot$depth, mu=out.NB1$estimate[1], size=out.NB1$estimate[2])
hostpot[order(indivi.p0),]
> hotspot[order(indivi.p0),]
id depth
21 S21 2
12 S12 118
31 S31 116
36 S36 111
14 S14 107
17 S17 39
32 S32 104
34 S34 103
8 S8 43
15 S15 99
4 S4 44
23 S23 94
28 S28 94
24 S24 93
发现泊松比NB要好，可能是这个情景下泊松对于均值左侧的更加敏感？
不过后续如果用诸如HMM或者贝叶斯之类的算法来找缺失和重复的话，应该用这两种模型都是没有问题的。
但是最近在使用的时候，发现这个算法在深度很低的时候，主要是检测杂合的时候，灵敏度降低。所以接下来又学习了一些东西，如果解决了这个问题，将写一篇新的博文。
参考文献：
[1] Statistical analysis of high-throughput sequencing count data
[2] 最大似然法
[3] 测序读断的负二项分布模型
[4] 负二项分布R实现
评论
发表评论