r語言在遺傳統計學中的應用

1、遺傳與疾病,人類的某些性狀及部分疾病與人體的遺傳因素密切相關闡明遺傳因素與人體疾病或健康狀態(tài)的關系有非常重要的意義遺傳統計學在這其中起著至關重要的作用,遺傳與疾病,疾病的易感基因研究,研究的特點,收集的數據即包含一般的表型數據也包含基因型數據數據分析時需要用到不同的遺傳模型需要一些遺傳統計特有的分析方法：LD的計算，家系圖的繪制等,R在遺傳統計中的應用,數據整理獲取位點的基本信息Hardy-Weinberg平衡檢驗

2、連鎖不平衡的計算關聯研究常用分析方法家系圖的繪制……,數據整理,R中的genetics包專門為基因型數據提供一個新的類 — genotypegenotype函數是genetics包里最基本的函數，可以將以下四種形式的初始基因型數據轉換成便于分析的帶有genotype類的數據,以一個字符分隔的向量 g1 <- genotype(c('C-C','C-T','C-C',&

3、#39;T-T', 'C-C',''),sep='-')2. 可以按某一位置分隔的向量 g2 <- genotype(c('DD','DI','DD','II',''),sep=1)3. 兩個分開的向量 al

4、lele1 <- c('D','D','D','I','') allele2 <- c('D','I','D','I','') g3 <- genotype(allele1, allele2),4. 數據框或矩陣中的兩列 data <

5、;- data.frame(allele1 = c('D','D','D','I',''), allele2 = c('D','I','D','I','')) g4 <- genotype(da

6、ta$allele1,data$allele2) 或 data1 <- cbind(allele1 = c('D','D','D','I',''), allele2 = c('D','I','D','I','&

7、#39;)) g5 <- genotype(data1),獲取位點的基本信息,多態(tài)位點的基本信息包括：位點分型成功率 (call rate)、等位基因頻率、基因型頻率、雜合度和多態(tài)信息含量 (PIC)一個簡單的例子：#載入popn數據data(popn,package="DGCgenetics") #獲取A位點的基本信息summary(popn$A),Number of samp

8、les typed: 1489 (96.9%)Allele Frequency: (2 alleles) Count Proportion1 1786 0.62 1192 0.4NA 94 NAGenotype Frequency: Count Proportion1/2 704 0.472/2 244 0.16

9、1/1 541 0.36NA 47 NAHeterozygosity (Hu) = 0.4802686Poly. Inf. Content = 0.3648558,Hardy-Weinberg定律,Hardy-Weinberg定律是由英國數學家哈迪(D.H. Hardy)和德國醫(yī)生溫伯格(W. Weinberg)于1908年分別獨立發(fā)現的，也稱遺傳平衡定律~(genetic equi

10、librium law)該定律可以簡單描述為，遺傳平衡群體的等位基因頻率與基因型頻率在世代間維持恒定該定律的適用條件是：隨機婚配，群體足夠大，沒有突變、選擇、遷移和遺傳漂變,Hardy-Weinberg平衡檢驗,關聯研究中Hardy-Weinberg平衡檢驗常被用來評價基因分型的質量。我們通常對病例和對照組分別進行Hardy-Weinberg平衡檢驗如果某一位點在對照組中不符合Hardy-Weinberg平衡，我們通常會懷疑

11、該位點的基因型鑒定的質量如果該位點在對照組平衡而在病例組出現不平衡，則該位點可能和疾病有關,Hardy-Weinberg平衡檢驗,genetics包里面提供兩種不同的檢驗方法一種是Pearson‘s chi-square test，可以用HWE.chisq函數進行該檢驗，另一種是Fisher exact test，對應于HWE.exact函數HWE.chisq常用于MAF較高、樣本量較大的場合；MAF較低的位點建議使用HWE

12、.exact函數,LD的計算,連鎖不平衡則是指人群中兩個位點處在同一個單體型的頻率比期望值高評價連鎖不平衡程度的指標包括D'、r2等genetics包提供計算LD各種指標的函數，并能以文字和圖形兩種形式顯示位點間的連鎖不平衡程度,LD的計算,#用LD函數計算位點間的LDldresult <- LD(popn) #用文字顯示D'值summary(ldresult, which="D'&

13、quot;)#用圖形顯示結果LDtable(ldresult, which = "D'")Pairwise LD----------- B C DA D' 0.979 0.976 0.976B D' 0.998 0.991C D'

14、 0.997,,關聯研究常用分析方法,卡方檢驗Logistic回歸線性回歸……,卡方檢驗,> data(popn,package="DGCgenetics") #首先載入popn數據> (geno chisq.test(geno)Pearson's Chi-squared testdata: genoX-squared = 23.7385, df = 2, p-va

15、lue = 7.003e-06,> (alle chisq.test(alle) Pearson's Chi-squared test with Yates' continuity correctiondata: alleX-squared = 23.6881, df = 1, p-value = 1.133e-06,Logistic回歸,1. 共顯性模型> summary(glm

16、(affected ~ A + sex, family=binomial, data=popn))Call:glm(formula = affected ~ A + sex, family = binomial, data = popn)Deviance Residuals: Min 1Q Median 3Q Max -1.4081 -1.2428 -0.6515 1

17、.1134 1.8190 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.6589 0.1354 -4.868 1.13e-06 *** A1/2 -0.3752 0.1234 -3.041 0.00236 **

18、A2/2 -0.7832 0.1695 -4.620 3.84e-06 ***sexFemale 1.1866 0.1335 8.890 < 2e-16 ***,Logistic回歸,2. 加性模型>summary(glm(affected ~ allele.count(A,'2') + sex, family=binomial, data=popn))

19、Call:glm(formula = affected ~ allele.count(A, "2") + sex, family = binomial, data = popn)Deviance Residuals: Min 1Q Median 3Q Max -1.410 -1.239 -0.655 1.117 1.814 Coeffici

20、ents: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.65396 0.13054 -5.010 5.45e-07 ***allele.count(A, "2") -0.38817 0.08107 -4.788 1.69e-06 ***s

21、exFemale 1.18676 0.13349 8.890 < 2e-16 ***,Logistic回歸,3. 顯性或隱性模型> summary(glm(affected ~ carrier(A,'2') + sex, family=binomial, data=popn))Call:glm(formula = affected ~ carrier(A, &q

22、uot;2") + sex, family = binomial, data = popn)Deviance Residuals: Min 1Q Median 3Q Max -1.4078 -1.1979 -0.7465 1.1571 1.6817 Coefficients:

23、Estimate Std. Error z value Pr(>|z|) (Intercept) -0.6566 0.1352 -4.857 1.19e-06 ***carrier(A, "2")TRUE -0.4788 0.1164 -4.115 3.87e-05 ***sexFemale 1.1835

24、0.1332 8.884 < 2e-16 ***,家系圖的繪制,library(kinship) #載入kinship包p1 <- scan(nlines=6,what=list(0,0,0,0,0,0)) 1 100 101 102 1 11 101 0 0 1 21 102 0 0 2 11 103 101 102 1 21 104 101 102 1 11 105 101 102 1 1ped &

25、lt;- as.data.frame(p1)names(ped) <- c("famid","id","dadid","monid","sex","aff")par(xpd=TRUE)ped1=pedigree(ped$id, ped$dadid, ped$monid, ped$sex, ped$aff)