r語言在遺傳統(tǒng)計(jì)學(xué)中的應(yīng)用

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

2、連鎖不平衡的計(jì)算關(guān)聯(lián)研究常用分析方法家系圖的繪制……,數(shù)據(jù)整理,R中的genetics包專門為基因型數(shù)據(jù)提供一個(gè)新的類 — genotypegenotype函數(shù)是genetics包里最基本的函數(shù)，可以將以下四種形式的初始基因型數(shù)據(jù)轉(zhuǎn)換成便于分析的帶有g(shù)enotype類的數(shù)據(jù),以一個(gè)字符分隔的向量 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. 兩個(gè)分開的向量 al

4、lele1 <- c('D','D','D','I','') allele2 <- c('D','I','D','I','') g3 <- genotype(allele1, allele2),4. 數(shù)據(jù)框或矩陣中的兩列 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),獲取位點(diǎn)的基本信息,多態(tài)位點(diǎn)的基本信息包括：位點(diǎn)分型成功率 (call rate)、等位基因頻率、基因型頻率、雜合度和多態(tài)信息含量 (PIC)一個(gè)簡單的例子：#載入popn數(shù)據(jù)data(popn,package="DGCgenetics") #獲取A位點(diǎn)的基本信息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定律是由英國數(shù)學(xué)家哈迪(D.H. Hardy)和德國醫(yī)生溫伯格(W. Weinberg)于1908年分別獨(dú)立發(fā)現(xiàn)的，也稱遺傳平衡定律~(genetic equi

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

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

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

13、quot;)#用圖形顯示結(jié)果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,,關(guān)聯(lián)研究常用分析方法,卡方檢驗(yàn)Logistic回歸線性回歸……,卡方檢驗(yàn),> data(popn,package="DGCgenetics") #首先載入popn數(shù)據(jù)> (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)