《Modern Statistics for Modern Biology》Chapter 一：離散數據模型的預測（1.1 - 1.3）

Modern Statistics for Modern Biology

書籍目錄

在分子生物學中拙徽，許多情況涉及計數事件：有多少密碼子使用某種拼寫图仓，有多少DNA的reads比對到參考序列上沈撞，DNA序列中CG的含量。這些計數為我們提供了離散的變量，而不是像質量和強度這樣的量是在連續(xù)的尺度上測量的洲脂。
如果我們知道所研究的機制所遵循的規(guī)則，即使結果是隨機的剧包，我們也可以通過計算和標準的概率定律來生成我們感興趣的任何事件的概率恐锦。這是一種自上而下的方法，基于演繹和我們對如何操縱概率的知識疆液。在第2章中一铅，您將看到如何將其與數據驅動(自底向上)統(tǒng)計建模相結合。

1.1 本章的目標

了解如何從給定的模型中獲得所有可能結果的概率堕油，并了解如何將理論頻率與在實際數據中觀察到的頻率進行比較潘飘。
探索一個完整的例子肮之，如何使用泊松分布分析數據的表位檢測。
如何對試驗離散數據生成最有效的模型：泊松卜录、二項和多項分布局骤。
使用R中相關函數計算概率和計數罕見事件。
通過指定的分布生成隨機數暴凑。

1.2 一個真實的案例

涉及markdown語法參考鏈接：

? CSDN-markdown語法之如何使用LaTeX語法編寫數學公式
? Markdown中公式的寫法(Latex)
　　　　　　　　　　　　　　　　　　 $\sqrt{x^3}$ : $\sqrt{x^3}$

我們被告知人類免疫缺陷病毒(HIV)基因組上的突變是隨機發(fā)生的峦甩，每個復制周期的核苷酸數為5×10^?4。這意味著现喳，在一個周期之后凯傲，大約10⁴=10 000個核苷酸的基因組中的突變數量將遵循泊松分布,我們將在稍后給出更多關于這種概率分布的詳細信息。利率為5嗦篱。這說明了什么冰单？這個概率模型預測，在一個復制周期中的突變數量將接近5個灸促，并且這個估計的可變性是 $\sqrt{5}$ (標準誤差)诫欠。我們現在已經有了我們期望在一個典型的HIV毒株中看到的突變數量及其變異性的基線參考值。

? 事實上浴栽，我們可以推斷出更詳細的信息荒叼。如果我們想知道在Poisson(5)模型下發(fā)生3個突變的頻率，我們可以使用一個R函數來生成看到x=3個事件的概率典鸡，取泊松分布的發(fā)生率參數的值被廓，稱為lambda(λ希臘字母，例如λ和μ通常表示表征我們使用的概率分布的重要參數)萝玷。 X表示給定區(qū)間的事件發(fā)生次數嫁乘，例如一個星期內的損壞次數。如果X符合泊松分布球碉，且每個區(qū)間內平均發(fā)生的λ次埠巨，或者說發(fā)生率為λ痴突，則寫作：X ~ Po(λ)

> dpois(x = 3, lambda = 5)
[1] 0.1403739

? 這就是說狼荞，看到三個事件的幾率大約是0.14拾积，也就是大約7個事件中就有1個斯碌。

如果要生成從0到12的所有值的概率傻唾，則不需要編寫循環(huán)冠骄。我們可以簡單地將第一個參數設置為這13個值的向量，使用R的序列操作符扁誓，冒號":"蝗敢。我們可以通過繪制它們來查看概率(圖1.1)前普。

> 0:12
 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12
> dpois(x = 0:12, lambda = 5)
 [1] 0.006737947 0.033689735 0.084224337 0.140373896 0.175467370 0.175467370
 [7] 0.146222808 0.104444863 0.065278039 0.036265577 0.018132789 0.008242177
[13] 0.003434240
> test <- dpois(x = 0:12, lambda = 5)
> barplot(test, names.arg = 0:12, col = "red")    ## 圖1.1

names.arg:  在每個條下出現的名稱的向量

圖1.1

? 數學理論告訴我們，看到值x的泊松概率是由公式

P(X=x)=\frac{\lambda^x}{x!}e^{-\lambda},k=0,1,…

給出的峻厚。在這本書中，我們將不時地討論理論辜王，但更傾向于顯示具體的數值例子和可視化呐馆，如圖 1.1续膳。

1.3 使用離散概率模型

? 點突變可以發(fā)生也可以不發(fā)生坟岔；它是一個二進制事件社付。兩種可能的結果(是, 否)稱為范疇變量的級別瘦穆。

? 并非所有事件都是二進制的。例如熙兔，一個二倍體生物的基因型可以分為三個水平(AA, Aa, aa)住涉。

? 有時舆声，一個分類變量中的水平數非常大；例如蛾找，生物樣本中不同類型的細菌的數量(數百或數千)和由3個核苷酸(64種)組成的密碼子的數量打毛。

? 當我們在一個樣本上測量一個分類變量時，我們經常想要統(tǒng)計一個計數矢量中不同級別的頻率熬甫。R對范疇變量有一種特殊的編碼罗珍，并將其稱為factor。在這里核无，我們抽取了不同的血液基因型的19名受試者繪總成一個表噪沙。

> genotype =  c("AA","AO","BB","AO","OO","AO","AA","BO","BO",
+              "AO","BB","AO","BO","AB","OO","AB","BB","AO","AO")
> table(genotype)
genotype
AA AB AO BB BO OO 
 2  2  7  3  3  2

? 在創(chuàng)建factor時，R會自動檢測局义。可以使用R中“l(fā)evels”函數功能查看level术幔。

> genotypeF = factor(genotype)
> levels(genotypeF)
[1] "AA" "AB" "AO" "BB" "BO" "OO"
> table(genotypeF)
genotypeF
AA AB AO BB BO OO 
 2  2  7  3  3  2

? 問題1.1：如果您想要創(chuàng)建一個尚未在數據中包含某些level的factor，該怎么辦拔妥？
factor(genotypeF, levels = c("AA", "AB", "AO", "BB", "BO", "OO", "new_levels"))

1.3.1 伯努利試驗

? 最為經典的是拋硬幣實驗癌蚁，拋開外在因素的干擾努释，投擲一枚硬幣，落地后的結果只有兩個：正面和反面恩沛。這是一個簡單的伯努利試驗雷客，輸出的結果所謂的伯努利隨機變量搅裙。R語言中binom用來做伯努利試驗計算裹芝。
? 假設我們要模擬15次投擲硬幣的結果部逮，將概率設置為0.5得到15次的結果。

> rbinom(15, prob = 0.5, size = 1)
 [1] 0 1 0 1 0 0 1 0 1 0 0 1 0 0 1

我們將rbinom函數與一組特定的參數66一起用于R函數嫂易，參數也稱為參數：第一個參數是我們要觀察的試驗次數兄朋；在這里我們選擇了15個。我們通過prob參數來指定成功的可能性怜械。按size=1（表示只有兩個結果0或者1）颅和，每個單獨的判斷只包括一次擲硬幣。

? 問題1.2：重復此函數調用多次灼芭。為什么答案不總是一樣的萤衰？
因為這是一個隨機分布函數，隨著重復的次數越來越多，最終的比例越接近剛開始指定的prob.

? 在伯努利試驗中成功和失敗可以為不相等的概率事件，只要概率和為1就行个束。為了模擬12次把球扔進兩個盒子的試驗殿如，如圖1.2所示，在如下右邊的盒子 $2 \over 3$ 和左邊的盒子 $1 \over 3$ 中

圖 1.2

> rbinom(15, prob = 0.5, size = 1)
 [1] 1 0 1 0 0 0 0 1 0 1 1 0 1 0 1

1表示成功讯泣，表示這個球落入到右邊的箱子中拳锚；0表示落入到左邊的箱子中棋凳。

1.3.2 Binomial success counts （翻譯不準勺良，故放原文）

如果我們只關心有多少個球落入到右邊的箱子中逻谦，而不關心它們落入的順序，我們只需就落入到右邊的箱子中的1加起來求和即可础米。因此而不是上面看到的向量: vector，我們只需要單一的一個數字結果惠拭。在R中只需將rbinom函數的size參數設置為12即可职辅。

> rbinom(1, prob = 2/3, size = 12)
[1] 9

這個結果告訴當我們將概率設置為 $2 \over 3$ 時候有多少個球會落到右邊的箱子中域携。當我們輸出結果只有兩種可能時拆内，我們可以使用隨機的two-box模型，比如：頭或者尾瓶颠、成功或者失敗欢搜、CpG或者non-CpG(甲基化或者未甲基化)、M或者F（男性還是女性）谴轮、Y = 嘧啶或者R= 嘌呤炒瘟、患病的或者健康的、正確或者錯誤第步。我們只需要計算成功的可能性p疮装，自然而然就得到了失敗的可能性1 - p。 SSSSSFSSSSFFFSF 總結為 (#Successes=10, #Failures=5), 或者看做x = 10, n = 15
在15次伯努利試驗中粘都，當成功的概率為0.3時結果為成功的次數就稱為一個二項式隨機變量或者一個服從B(15, 0.3)分布的隨機變量廓推。下面在R中使用函數rbinom將試驗次數設置為15：

> set.seed(235569515) ## 上面不是提到不同人運行結果不一樣，這里設置seed后大家運行的結果就一致了翩隧。
> rbinom(1, prob = 0.3, size = 15)
[1] 5

? 問題 1.3 ：將此函數調用重復十次樊展。最常見的結果是什么？
4堆生，因為畢竟這是一個服從B(15, 0.3)的二項分布函數

The complete probability mass distribution is available by typing:

> probabilities = dbinom(0:15, prob = 0.3, size = 15)
> round(probabilities, 2)
 [1] 0.00 0.03 0.09 0.17 0.22 0.21 0.15 0.08 0.03 0.01 0.00 0.00 0.00 0.00 0.00 0.00

繪制上面數據的分布圖

barplot(probabilities, names.arg = 0:15, col = "red")

圖 1.3 可以看到頻率最高的為4

我們通常在R中通過指定size 參數來確定試驗的次數专缠，并且經常用n表示，成功的概率用p表示淑仆。使用數學術語來將X是服從參數(n, p)的二項分布寫作 X ~ B(n, p), X成功的概率為：

? 問題 1.4 當k = 3, p = $2\over3$ , n = 4時涝婉，結果為什么？

> dbinom(3, prob = 2/3, size = 4)
[1] 0.3950617

小插曲: R計算二項Binomial分布的P
pbinom(q,size,prob): 計算累積概率
dbinom(x,size,prob): 計算取某個值的特定概率
rbinom(n, size, prob): 產生n個b(size,prob)的二項分布隨機數
qbinom(p, size, prob): quantile function 分位數函數蔗怠。

#x 
     vector of (non-negative integer) quantiles.
     分位數（非負整數）向量墩弯。
#q 
     vector of quantiles.
     分位數向量。
#p 
     vector of probabilities.
     概率向量寞射。
#n 
     number of random values to return.
     返回隨機值的數量渔工。
#lambda 
     vector of (non-negative) means.

     均值（非負）向量。

#log, log.p 
     logical; if TRUE, probabilities p are given as log(p).
     邏輯型桥温，如果為真引矩，概率p作為log(p)給定。

#lower.tail 
     logical; if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].
     邏輯型策治，如果為真（默認）脓魏，概率是P[X<=x], 否則概率為P[X>x].

假定在一次實驗中，事件A發(fā)生的概率為0.1通惫，那么進行11次這樣的實驗茂翔，觀察到4次的概率是多少？
答案是Binomial Probability: P(X = 4)= dbinom(4,11,0.1)
[1] 0.0157838

觀察到小于等于4次的概率是多少履腋？
Cumulative Probability: P(X < =4)= pbinom(4,11,0.1)
[1] 0.997249

1.3.3 泊松分布

當我們的試驗次數n很大切成功的概率p很小的時候, 二項分布B(n, p)可以完全近似的認為是一個更簡單λ = np的泊松分布珊燎，在圖 1.1 的HIV試驗中涉及惭嚣。
? 問題 1.5
當每個堿基突變的概率為p = 5 * 10^-4 時候，基因組上n=10⁴個堿基觀察到0:12個突變的分布是什么? 當用二項分布B(n, p)和泊松分布（λ = np）為模型時悔政，結果是否相似晚吞？
請注意，與二項分布不同谋国，泊松分布不再依賴于兩個獨立的參數n和p槽地，而只依賴于它們的乘積np。在二項分布的情況下芦瘾，我們也有一個計算泊松概率的數學公式捌蚊。
　　　　　　　　　　　　　　　　 $P(X=k)=\frac{\lambda^k}{k!}e^{-\lambda},k=0,1,…$
例如：將 λ = 5代入，計算P(X = 3)

> 5^3 * exp(-5) / factorial(3)
[1] 0.1403739

> dpois(3, 5)
[1] 0.1403739

可以使用dpois(x, λ)來計算近弟。
同樣的上面堿基突變問題缅糟，模擬10000個位點的突變過程，突變率為5×10^-4祷愉，并計算突變數窗宦。多次重復此操作，并使用barlot函數繪制分布二鳄。

> rbinom(1, prob = 5e-4, size = 10000)
[1] 2
> simulations = rbinom(n = 300000, prob = 5e-4, size = 10000)
> barplot(table(simulations), col = "lavender")

1.3.4 A generative model for epitope detection

在測試某些藥物化合物時赴涵，檢測引起過敏反應的蛋白質是很重要的。負責這些反應的分子位點被稱為表位泥从。表位的技術定義是：
　　抗體與大分子抗原結合的特定部分句占。對于T細胞識別的蛋白質抗原，表位或決定簇是與主要組織相容性復合物(MHC)分子結合以供T細胞受體(TCR)識別的肽段或位點躯嫉。
如果你對免疫學不是很熟悉的話：抗體(如圖1.6所示)是由某些白細胞產生的一種蛋白質，對身體中的一種外來物質作出反應杨拐，這就是所謂的抗原祈餐。

　　

抗體圖，顯示幾種顏色的免疫球蛋白結構域
抗體與其抗原結合(或多或少具有特異性)哄陶。結合的目的是幫助破壞抗原帆阳。抗體可以通過幾種方式發(fā)揮作用屋吨，取決于抗原的性質蜒谤。有些抗體直接破壞抗原。其他人幫助招募白血球來破壞抗原至扰△⒒眨抗原表位，也稱為抗原決定簇敢课，是被免疫系統(tǒng)識別的抗原的一部分阶祭，特別是被抗體绷杜、B細胞或T細胞識別。

已知參數的`ELISA`誤差模型

ELISA酶聯免疫吸附試驗濒募。檢測方法用于檢測蛋白質不同位置的特定表位鞭盟。假設我們使用的ELISA陣列具有以下事實：

每個位置的基線噪聲水平，或者更準確地說瑰剃，假陽性率為1%齿诉。這是表示一個成功的可能性-當沒有的時候我們認為我們有一個表位。我們寫這個P( declare epitope |no epitope)
這種蛋白質是在100個不同的位置測試的晌姚，應該是獨立的鹃两。
我們將對50個病人樣本進行檢查。

一個病人的數據

一個病人的芯片數據如下

##   [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
##  [30] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [59] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
##  [88] 0 0 0 0 0 0 0 0 0 0 0 0 0

其中舀凛，1表示擊中(因此可能發(fā)生過敏反應)俊扳，0表示在該位置沒有反應。
通過仿真驗證了p = 0.01的50個獨立伯努利變量之和與Poission(0.5)隨機變量相同猛遍，具有足夠好的逼近性馋记。

50個芯片數據的結果

我們將研究所有50名患者在100個位置上的數據。如果沒有過敏反應懊烤，假陽性率意味著對于一個病人來說梯醒，每一個位置都有1/100的概率是1。因此腌紧，在對50名患者進行統(tǒng)計后茸习，我們期望在任何給定的位置，50個觀察到的(0壁肋，1)個變量之和有一個參數為0.5的泊松分布号胚。典型的結果可能如圖1.7所示。現在假設我們看到了如圖1.8所示的實際數據浸遗，它作為R對象e100從數據文件e100.RData加載猫胁。

> setwd(dir = "E://compute language/NGS/Modern Statistics for Modern Biology/data")
> load("e100.RData")
> barplot(e100, ylim = c(0, 7), width = 0.7, xlim = c(-0.5, 100.5),
+   names.arg = seq(along = e100), col = "darkolivegreen")

圖1.8中的峰值是驚人的。如果沒有表位跛锌，那么看到7這樣大的值的可能性有多大弃秆？如果我們考慮一個泊松(0.5)隨機變量時，看到一個7(或更大)的數字的概率髓帽，那么答案可以用如下的封閉形式來計算:
這當然與1?P(X≤6)相同菠赚。概率P(X≤6)是6的所謂累積分布函數，R有用來計算它的函數ppois郑藏，我們可以用以下兩種方法中的任何一種來使用它：

> 1 - ppois(6, 0.5)
[1] 1.00238e-06

> ppois(6, 0.5, lower.tail = FALSE)
[1] 1.00238e-06

小插曲： ppois是泊松分布的分布函數（即用來求累計概率）衡查，因為是離散的，所以只會在整數左右有變化, R語言中ppois如何應用
ppois(q, lambda, lower.tail = TRUE, log.p = FALSE)
q：官方幫助文檔說是分位數译秦，指定x軸上的點峡捡。
lambda：就是泊松分布的參數λ
lower.tail：是邏輯變量击碗，當它為真(TRUE，缺省值)時们拙，分布函數的計算公式為

當lower.tail = FALSE時稍途，分布函數的計算公式為

log,log.p是邏輯變量，當它為真(TRUE)時砚婆，函數的返回值不再是泊松分布械拍，而是對數泊松分布.
比如lambda=1的分布函數作圖如下：
在假設沒有表位反應的情況下，看到7的概率是：

Poisson分布的極值分析

停装盯，停止!坷虑。在這種情況下，上述計算不是正確的計算埂奈。

? 問題1.6
如果我們想要計算在沒有表位的情況下觀察到這些數據的概率迄损，你能發(fā)現我們推理中的缺陷嗎？

因此账磺，我們不應該問芹敌，看到泊松(0.5)大到7的概率有多大，而應該問自己垮抗，100泊松(0.5)試驗的最大值達到7的可能性有多大氏捞？我們在這里使用極值分析，這意味著我們對隨機分布的非常大或非常小的值的行為感興趣冒版，例如液茎，最大值或最小值。我們對數據值x₁辞嗡、x₂捆等，x₃,...,x₁₀₀進行排序，并且重新命名為x₍₁₎欲间、x₍₂₎楚里，x₍₃₎,...,x₍₁₀₀₎, 使得x₍₁₎表示100個樹中上最小值，x₍₁₀₀₎表示最大值猎贴。x₍₁₎...x₍₁₀₀₎稱為這個樣本的100個值的秩統(tǒng)計量。
最大值為7是所有100個計數都小于或等于6的互補事件蝴光。兩個互補事件的概率之和為1她渴。因為位置應該是獨立的，所以我們現在可以進行計算了蔑祟。
因為我們假設這100個事件中的每個事件都是獨立的趁耗，所以我們可以使用上面的結果。

真實的計算數字

我們可以讓R計算這個數字的值：(1??)¹⁰⁰疆虚。對于那些對這些計算如何通過近似值快捷的感興趣的人苛败，我們給出了一些細節(jié)满葛。這些可以在第一次閱讀時跳過。
我們從上面回憶起罢屈，? ? 10^?6比1小得多嘀韧。為了近似計算(1??)¹⁰⁰的值，我們可以使用二項式定理缠捌，去掉?的所有“高階”項锄贷，即?²，?³,...的所有項, 因為與其余的(“領先”)術語相比曼月，它們小得可以忽略不計谊却。
另一個等價的路由是使用近似的e^?? ? 1??，與log(1??)???一樣哑芹。因此
因此炎辨，如果沒有表位，在100個位置中看到一些大小或大于7的點擊的正確概率大約是我們以前錯誤計算的概率的100倍聪姿。
10^?6和10^?4的計算概率均小于標準顯著性閾值(0.05碴萧、0.01或0.001)。拒絕任何表位無效的決定也是一樣的咳燕。然而勿决，如果一個人必須在法庭上站起來，像在一些法庭案件中那樣招盲，將p值辯護到8個有效數字(1717年低缩，這發(fā)生在OJ Simpson案的法醫(yī)證據審查中)。那就是另一回事了曹货∨胤保考慮到測試多重性的調整后的p值是應該報告的值，我們將在第6章中回到這一重要問題顶籽。

模擬計算概率

在我們剛才看到的情況下玩般，理論上的概率計算是相當簡單的，我們可以通過一個明確的計算來計算結果礼饱。在實踐中坏为，事情往往更復雜，我們更好地計算我們的概率使用蒙特卡洛方法：一種基于我們的生成模型的計算機模擬镊绪，找到我們感興趣的事件的概率匀伏。下面，我們將生成100000個從100個泊松分布數中選取最大值的實例蝴韭。

> maxes = replicate(100000, {
+   max(rpois(100, 0.5))
+ })
> table(maxes)
maxes
    1     2     3     4     5     6     7     9 
    7 23059 60817 14355  1605   141    15     1

在100000個試驗中够颠，有16個試驗的最大值在7個或7以上。這給出了P(X_max ≥ 7)的以下近似值：

> mean( maxes >= 7 )
[1] 0.00016

這與我們的理論計算基本一致榄鉴。我們已經看到蒙特卡羅模擬的一個潛在限制：模擬結果的“粒度”由模擬次數(100000)的倒數決定履磨，因此將大約為 $10^{-5}$ 蛉抓。任何估計的概率都不能比這個粒度更精確，事實上剃诅，我們估計的精度將是這個粒度的幾倍巷送。到目前為止，我們所做的一切都是可能的综苔，因為我們知道每個位置的假陽性率惩系，我們知道被檢測的病人人數和蛋白質的長度，我們假設我們從模型中得出了相同分布的獨立圖如筛，我們假設我們有來自模型的同分布的獨立畫法堡牡，并且沒有未知的參數。這是一個概率或生成建模的例子：所有的參數都是已知的杨刨，數學理論允許我們以一種自上而下的方式進行推理晤柄。
相反，如果我們在更現實的情況下知道病人的數量和蛋白質的長度妖胀，但不知道數據的分布芥颈，那么我們就必須使用統(tǒng)計模型。這一方法將在第2章中闡述赚抡。我們將看到爬坑，如果我們只有數據開始，我們首先需要適合一個合理的分布來描述它涂臣。然而盾计，在我們討論這個更困難的問題之前，讓我們將離散分布的知識擴展到不僅僅是二進制的赁遗、成功或失敗的結果署辉。

《Modern Statistics for Modern Biology》Chapter 一： 離散數據模型的預測（1.1 - 1.3）