R語言的遺傳模塊

1、我接觸 r 的 算是不短了，已 兩年多了。期 斷斷 的看了些 r 網(wǎng)站上的材料。 在已 了用 r 做數(shù)據(jù)分析了，并且越來越喜 用 r 來做分析了。之前我用 sas ， spss 也試過 stata ，但是 三個 件都沒有 的 模 （至少國內(nèi)流行的盜版里沒有）。所以和其它 相比，我想 r 我 也 更有用些。cos 里提到r 在 方面的packagegenetic statistics里的 用的帖子很少。我在 里寫一些我平 用到的 的 明，一來算是個人小 再者算是拋 引玉吧，希望cos 里的各位多寫些相關(guān)的 西。introduction. cran task view: statistical g

2、eneticscran task view當(dāng)中有一個 獨的package和相關(guān) 接。 足可以看出genetics部分，里面列出了40 個 相關(guān)的r 在 學(xué)當(dāng)中的影響和作用。里面核心的 core package 有以下三個 : genetics, gap, 和 haplo.stats 。 有一個我 常用到的包是 dgcgenetics ，算是 genetics 包的 展。以后我會提到以上幾個包里面的一些函數(shù)。大致包括以下幾方面的內(nèi)容：1. 以上幾個 package 數(shù)據(jù)格式的要求；2. 多 位點的基本信息（ maf 等）；3. hardy-weinberg 平衡 ；4. ld 的 算；5. 關(guān)

3、研究常用 方法；6. power 的 算；先 一下前面提到的幾個包的安裝吧，其 很 。一個一個用當(dāng)然可以。相 點的方法是用cran task views里提到的install.packages()ctv 包來批量安裝。函數(shù)來安裝install.packages(ctv) #首先安裝 ctv packagelibrary(ctv) # 入 ctv packageinstall.views(genetics,coreonly = true) #安裝的包。如果不加core.only=true 會安裝所有的genetics, gap, haplo.stats40 個 相關(guān)的三個核心包及依 packag

4、e 。install.packages (genetics, coreonly = true)dgcgenetics包的下 地址是oftware/dgcgenetics_1.0.zip。你需要先下 個包，然后本地安裝。方法大家 都知道，數(shù)據(jù)格式（ 1）rgui的packages菜 的install package(s) from local zip files。遺傳研究收集的數(shù)據(jù)有自己的特點。往往是數(shù)據(jù)集中即包含一般的表型數(shù)據(jù)（分類和連續(xù)變量；如血壓水平， bmi 和性別等），又包括基因型數(shù)據(jù)。分析時往往還需要用到不同的遺傳模型，什么顯性、隱形、加性模型，或者是按照分類變量來處理（有時候也稱為

5、共顯性模型）。用 sas 或 spss 分析遺傳數(shù)據(jù)時，如果要用不同的遺傳模型進行數(shù)據(jù)分析的話，必須先進行數(shù)據(jù)轉(zhuǎn)換，過程相對復(fù)雜。r 中的genetics包專門為基因型數(shù)據(jù)提供了一個新的class（類），你可以很方便的用genotype() 或 makegenotypes()函數(shù)將不同形式的初始基因型數(shù)據(jù)轉(zhuǎn)換成基因型數(shù)據(jù)，并為數(shù)genotypegeneticssummary.genotype()plot.genotype()函數(shù)。你可以很方便的用summary() 函數(shù)獲取基因型數(shù)據(jù)的基因型頻率、等位基因頻率等基本信息，用plot() 函數(shù)做出基因型的柱狀圖。先說一下 genotype() 函

6、數(shù)，該函數(shù)是基因型數(shù)據(jù)轉(zhuǎn)換成便于分析的帶有g(shù)enetics包里最基本的函數(shù)?？梢詫⒁韵滤姆N形式的初始genotype class的數(shù)據(jù)。1. 以一個字符分隔的向量e.g.g1 - genotype(c(d/d,d/i,d/d,i/i,d/d,na)g2 - genotype(c(c-c,c-t,c-c,t-t,c-c,),sep=-)2. 可以按某一位置分隔的向量e.g.g3 - genotype(c(dd,di,dd,ii,),sep=1)#sep=1表示在位置1 后分成兩個allele3. 兩個分開的向量e.g.allele1 - c(d,d,d,i,)allele2 - c(d,i,d,

7、i,)g4 - genotype(allele1, allele2)4. 數(shù)據(jù)框或矩陣中的兩列data - data.frame(allele1 = c(d,d,d,i,), allele2 = c(d,i,d,i,)g5 - genotype(data$allele1,data$allele2)ordata1 - cbind(allele1 = c(d,d,d,i,), allele2 = c(d,i,d,i,)g6 - genotype(data1)實際的數(shù)據(jù)分析過程中建議將每一個snp 位點的基因型數(shù)據(jù)按照等位基因 /等位基因（e.g.a/a, a/t, t/t） 的格式給出，而不要簡單

8、的用數(shù)字表示。這樣有兩個好處，一是可以很方便的用 makegenotypes() 函數(shù)一次性地將多個位點的原始基因型數(shù)據(jù)轉(zhuǎn)換成帶有 genotype 類屬性的基因型數(shù)據(jù)，二是便于數(shù)據(jù)分析過程中了解具體是哪一個等位基因和疾病或性狀有關(guān)。如果用數(shù)字的話，位點數(shù)目一多，可能就完全糊涂了。舉個實例演示一下：library(genetics)#用 scan() 函數(shù)讀入16 個人的數(shù)據(jù)g1 - scan(nline = 16, what=list(0,0,0,0,)1 45 1 31.5 a/a c/c2 39 2 24.5 t/t c/g3 36 1 23 a/t c/c4 41 2 26 a/t c

9、/c5 37 1 29.5 a/a c/c6 35 2 31 a/t g/g7 41 2 21.5 a/a c/g8 43 2 27.5 a/a g/g9 44 2 24 a/a c/c10 32 1 19.5 a/t c/c11 40 2 20 a/a c/g12 38 2 22.5 t/t g/g13 42 2 32.5 a/a c/c14 33 2 25.5 a/t c/c15 43 1 30.5 a/t g/g16 35 2 25 a/t c/cg1 - as.data.frame(g1)names(g1) - c(id, age, gender, bmi, snp1, snp2)g1

10、$gender - factor(g1$gender, labels=c(male,female)#用 makegenotypes()函數(shù)將g1 中的兩列基因型數(shù)據(jù)附上genotype類屬性g2 - makegenotypes(g1)#大功告成，可以用str() 和 summary() 看看 g1 和 g2 的區(qū)別str(g1); str(g2)summary(g1); summary(g2)獲取多態(tài)位點的基本信息我用 dgcgenetics包里面的popn 數(shù)據(jù)為例子，介紹一下獲取多態(tài)位點基本信息的函數(shù)。data(popn,package=dgcgenetics) #首先載入popn 數(shù)據(jù)h

11、ead(popn) # 該數(shù)據(jù)包含四個多態(tài)位點（ a, b, c, and d ）、性別（ sex ）、疾病狀態(tài)（ affect ）及 id 號（ subject ）。summary(popn$a) #獲取 a 位點的基本信息，包括該位點分型成功率（call rate ）、等位基因頻率、基因型頻率、雜合度和多態(tài)信息含量（pic ）number of samples typed: 1489 (96.9%) #call rateallele frequency: (2 alleles) #等位基因頻率count proportion117860.6211920.4na94nagenotype fr

12、equency: #基因型頻率count proportion1/27040.472/22440.161/15410.36na47naheterozygosity (hu) = 0.4802686 #雜合度poly. inf. content= 0.3648558#pichardy-weinberg平衡檢驗首先簡單介紹一下hardy-weinberg定律。該定律是由英國數(shù)學(xué)家哈迪（d.h.hardy ）和德國醫(yī)生溫伯格（ w.weinberg）于 1908 年分別獨立發(fā)現(xiàn)的，也稱遺傳平衡定律（geneticequilibrium law ）。該定律可以簡單描述為，遺傳平衡群體的等位基因頻率與基

13、因型頻率在世代間維持恒定。該定律的適用條件是：隨機婚配，群體足夠大，沒有突變、選擇、遷移和遺傳漂變。在關(guān)聯(lián)研究中hardy-weinberg平衡檢驗常被用來評價基因分型的質(zhì)量。我們通常對病例和對照組分別進行hardy-weinberg平衡檢驗。如果某一位點在對照組中不符合hardy-weinberg平衡，我們通常會懷疑該位點的基因型鑒定的質(zhì)量。如果該位點在對照組平衡而在病例組出現(xiàn)不平衡，則該位點很可能和疾病有關(guān)。genetics包里面提供兩種不同的檢驗方法：一種是pearsons chi-square test，可以用hwe.chisq()函數(shù)進行該檢驗，另一種是fisher exact te

14、st，對應(yīng)于hwe.exact()函數(shù)。hwe.chisq() 常用于 maf 較高、樣本量較大的場合。 maf 較低的位點建議使用 hwe.exact() 函數(shù)。library(genetics)data(popn, package=dgcgenetics)control - popn$affected=controlcase - popn$affected=casehwe.chisq(popn$acontrol)hwe.exact(popn$acontrol)hwe.chisq(popn$a case )hwe.exact(popn$a case )ld的計算連鎖不平衡是指人群中兩個位點處

15、在同一個單體型的頻率比期望值高。評價連鎖不平衡程度的指標包括 d 、 r2 等。 genetics包提供計算ld 各種指標的函數(shù)，并能以文字和圖形兩種形式顯示位點間的連鎖不平衡程度。data(popn,package=dgcgenetics)# 首先載入popn 數(shù)據(jù)ldresult - ld(popn)#用 ld 函數(shù)計算位點間的ldsummary(ldresult, which=d ) # 用文字顯示d值ldtable(ldresult, which =d ) # 用圖形顯示結(jié)果結(jié)果如下：pairwise ldtable=50%trtd/tdtdb/tdtdc/tdtdd/td/trtrt

16、da d/tdtd0.978/tdtd0.976/tdtd0.976/td/trtrtdb d/tdtd/tdtd0.998/tdtd0.991/td/trtrtdc d/tdtd/tdtd/tdtd0.997/td/tr/table 便幫樓主完成haplo的一個函數(shù)，用在pool的 域。函數(shù)：haplo(n)用于生成n 個loci的haplotype configuration matrix;（一 ）所有可能haplotype00,0,1,1,0,1,1的矩 ，因 循 次數(shù)達到了o(n*2n), 所以用 c 言寫的，用r 用。附件中.so 文件是 linux 版本， .dll 是 windo

17、ws版本的（今天沒有 限再上 附件了把 c 代 附加在最后）。r 用的代 如下：dyn.load(/code/haplo.so)#用者自酌haplo-function(n)densa- .c(haplo,eger(n),result=eger(vector(integer,n*2n)densaresultc 的代 如下：#include#define denotvoid haplo(int *q,int *result)int i,j,tmp;/*int r=(2(*q);cannot use in r, if q5 may cause flea*/int r=1;fo

18、r (i=1;i=*q;i+)r*=2;for (i=0;i(*q);i+)for(j=0;jr;j+)resulti*r+j=0;for (j=0;jr;j+)tmp=j;for (i=0;i=1)result(*q-i-1)*r+j=tmp%2;tmp/=2;#ifndef denotfor (i=0;i(*q);i+)printf(n);for(j=0;jr;j+)printf(%dt,resulti*r+j);#endifkinship2基因遺傳學(xué) s:11library(kinship2)data(sample.ped)sample.ped1: 10 ,pedall - pedigr

19、ee(id=sample.peddadid=sample.pedsex=sample.pedprint(pedall)ped1basic - pedall1ped2basic - pedall2print(ped1basic)print(ped2basic)plot(ped2basic)# plot(ped1basic)kin2 - kinship(ped2basic)kin2pedall - pedigree(id=sample.ped$id ,$father, momid=sample.ped$sex , famid=sample.ped$ped )$id ,$mother ,dadid=

20、sample.pedsex=sample.ped$father, momid=sample.ped$sex , famid=sample.ped$ped )$mother ,kinall - kinship(pedall)kinall1: 14 , 1: 14kinall40: 43, 40 : 43 kinall42: 46,42: 46df2 - sample.pedsample.ped$ped = 2,names(df2)df2 $censor- c(1, 1, rep(0,12)ped2 - pedigree(df2$id , df2$father, df2$mother ,df2$sex , status=df2$censor)ped2 - pedigree(df2$id , df2$father, df2$mother ,df2$sex , affected=df2$affected,status=df2$censor)aff2 - data.frame(blue=df2$affected,bald=c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1)ped2 - pedigree(df2$id , df2$father, df2$mother ,df2$sex , affected=as.matrix(aff2),