PH525x series - Statistical Models（一）

本章作者介绍了除正态分布以外的几个常用的统计模型，话不多说，开始！

二项分布

再重复的独立试验中，如果每个试验仅可能为两种可能结果之一，而且得到相应结果的概率在每次试验中保持不变，则成这些重复独立试验为Bernoulli 试验。

我们把在成功概率为 $p$ 的 $N$ 次Bernoulli试验中成功次数 $S$ 称为具有参数 $(N,p)$ 的二项分布随机变量。而其概率分布为：
$\mbox{Pr}(S=k) = {N \choose k}p^k (1-p)^{N-k}$

$S/N$ 是成功次数 $S$ 的平均数，因此根据极限中心理论，当 $N$ 很大时， $S$ 大致满足正态分布;

应用

二项分布经常用于变异或基因分型的相关分析中

泊松分布

泊松分布其实大致相当于在大量重复的成功概率很小的独立Bernoulli试验中成功次数的分布，比如说某宇宙粒子击中一个卫星的次数、一种疾病发生的次数等等都属于泊松分布，其概率分布函数为：

$\mbox{Pr}(S=k)=\frac{\lambda^k \exp{-\lambda}}{k!}$

泊松分布的总体均值与方差相等，即 $\lambda$ ,等于 $N*p$ ，在生命科学领域，这种分布可以用在RNAseq分析当中。

举例说明：
假如在RNAseq中样本分为case和control两组，且每组仅有1个样本，现在我们想要研究基因的差异表达倍数。


N=10000##number of genes
lambdas=2^seq(1,16,len=N) ##these are the true abundances of genes
#先假设对于所有基因来说，零假设都未真
y=rpois(N,lambdas)
x=rpois(N,lambdas) 
ind=which(y>0 & x>0)##make sure no 0s due to ratio and log

library(rafalib)
splot(log2(lambdas),log2(y/x),subset=ind)

lambdas1.png

rpois(N,lambdas)函数的作用是生成一组满足泊松分布的数据， $N$ 是试验的次数，在本例中就是基因的个数， $lambdas$ 是成功的次数，在本例中具体是指基因的表达量。

splot()函数是更加智能的绘图plot()函数

从上图中我们可以看到，当零假设全部为真时，基因丰度越低，两样本间的差异表达倍数变化波动就越大，有不少值竟已大于2，这也就是说，在真实转录组分析中，表达倍数大于2的那些低丰度的基因中有很多并不是真正的差异表达基因，而是假阳性。

泊松分布在NGS中的应用

如果我们将泊松分布应用在NGS read count的分析中，就会发现，实际数据的方差比使用模型预测到的方差大，这是由于泊松分布并没有将生物变异相关的因素考虑在内，现在我们就parathyroidSE包中的数据来具体说明：

library(parathyroidSE)
library(rafalib)
library(matrixStats)

data(parathyroidGenesSE)
se <- parathyroidGenesSE
vars=rowVars(assay(se)[,c(2,8,16,21)]) ##we now these four are 4
means=rowMeans(assay(se)[,c(2,8,16,21)]) ##different individulsa

splot(means,vars,log="xy",subset=which(means>0&vars>0)) ##plot a subset of data
abline(0,1,col=2,lwd=2)

parathyroidGenesSE.png

parathyroidGenesSE数据是保存在SummarizedExperiment object当中的，可以用assay()函数进行提取，其中，该数据的行是某种基因组特征，列是样本，每个单元格内的数据是reads数。

我们已知在泊松分布中，均值与方差是相等的，所以满足泊松分布的数据，点 $(mean,var)$ 应该落于对角线上，但实际上却并不是。实际上负二项分布更加适合这种类型的分析，因为这种分布将泊松分布的采样变异与生物变异结合起来了，具体可以看这篇文章。

最大似然估计（MLE）

由于总体均值 $\lambda$ 其实是个未知参数，所以为了了解总体的特征，我们就需要使用最大似然估计法对该参数进行估计：

首先按照如下公式为每个观察值计算似然值：
$\mbox{L}(\lambda; X_1=k_1,\dots,X_n=k_1)=\exp\left\{\sum_{i=1}^n \log \Pr(X_i=k_i;\lambda)\right\}$

然后求出这一组似然值中的最大值，此时的 $\lambda$ 即为MLE。现在我们用一套统计HMCV基因组回文序列坐标的数据进行说明:

datadir="http://www.biostat.jhsph.edu/bstcourse/bio751/data"
x=read.csv(file.path(datadir,"hcmv.csv"))[,2]

breaks=seq(0,4000*round(max(x)/4000),4000)
tmp=cut(x,breaks)
counts=table(tmp)

library(rafalib)
mypar(1,1)
hist(counts)

palindrome.png

很显然，回文序列在基因组上出现的次数是满足泊松分布的，为了求 $\lambda$ 接下来：

l<-function(lambda) sum(dpois(counts,lambda,log=TRUE)) 

lambdas<-seq(3,7,len=100)
ls <- exp(sapply(lambdas,l))

plot(lambdas,ls,type="l")

mle=optimize(l,c(0,10),maximum=TRUE)
abline(v=mle$maximum)

MLE.png

cut(x, breaks)函数的作用是将 $x$ 按照breaks进行划分，返回的是x中的每一个元素在breaks中的区间。

exp(x)函数是求 $e$ 的 $x$ 次方；

optimize(f, interval, ... , maximum = TRUE)函数是求区间interval中使得函数 $f$ 最大的值。

因此，上图中的竖直辅助线所对应的横坐标便是MSE。另外,若你计算过counts的均值的话便会发现MSE其实就是counts的均值。

print(c(mle$maximum, mean(counts)))

## [1] 5.157894 5.157895

通过 QQ图还可以看到，在本例中，泊松分布对于结果的拟合是相当好的：

theoretical<-qpois((seq(0,99)+0.5)/100,mean(counts))
qqplot(theoretical,counts)
abline(0,1)

theoretical.png

PH525x series - Statistical Models（一）

二项分布

泊松分布

泊松分布在NGS中的应用

最大似然估计（MLE）

推荐阅读更多精彩内容