R语言的遗传模块

1、我接触 r 的 算是不短了，已 两年多了。期 断断 的看了些 r 网站上的材料。 在已 了用 r 做数据分析了，并且越来越喜 用 r 来做分析了。之前我用 sas ， spss 也试过 stata ，但是 三个 件都没有 的 模 （至少国内流行的盗版里没有）。所以和其它 相比，我想 r 我 也 更有用些。cos 里提到r 在 方面的packagegenetic statistics里的 用的帖子很少。我在 里写一些我平 用到的 的 明，一来算是个人小 再者算是抛 引玉吧，希望cos 里的各位多写些相关的 西。introduction. cran task view: statistical g

2、eneticscran task view当中有一个 独的package和相关 接。 足可以看出genetics部分，里面列出了40 个 相关的r 在 学当中的影响和作用。里面核心的 core package 有以下三个 : genetics, gap, 和 haplo.stats 。 有一个我 常用到的包是 dgcgenetics ，算是 genetics 包的 展。以后我会提到以上几个包里面的一些函数。大致包括以下几方面的内容：1. 以上几个 package 数据格式的要求；2. 多 位点的基本信息（ maf 等）；3. hardy-weinberg 平衡 ；4. ld 的 算；5. 关

3、研究常用 方法；6. power 的 算；先 一下前面提到的几个包的安装吧，其 很 。一个一个用当然可以。相 点的方法是用cran task views里提到的install.packages()ctv 包来批量安装。函数来安装install.packages(ctv) #首先安装 ctv packagelibrary(ctv) # 入 ctv packageinstall.views(genetics,coreonly = true) #安装的包。如果不加core.only=true 会安装所有的genetics, gap, haplo.stats40 个 相关的三个核心包及依 packag

4、e 。install.packages (genetics, coreonly = true)dgcgenetics包的下 地址是oftware/dgcgenetics_1.0.zip。你需要先下 个包，然后本地安装。方法大家 都知道，数据格式（ 1）rgui的packages菜 的install package(s) from local zip files。遗传研究收集的数据有自己的特点。往往是数据集中即包含一般的表型数据（分类和连续变量；如血压水平， bmi 和性别等），又包括基因型数据。分析时往往还需要用到不同的遗传模型，什么显性、隐形、加性模型，或者是按照分类变量来处理（有时候也称为

5、共显性模型）。用 sas 或 spss 分析遗传数据时，如果要用不同的遗传模型进行数据分析的话，必须先进行数据转换，过程相对复杂。r 中的genetics包专门为基因型数据提供了一个新的class（类），你可以很方便的用genotype() 或 makegenotypes()函数将不同形式的初始基因型数据转换成基因型数据，并为数genotypegeneticssummary.genotype()plot.genotype()函数。你可以很方便的用summary() 函数获取基因型数据的基因型频率、等位基因频率等基本信息，用plot() 函数做出基因型的柱状图。先说一下 genotype() 函

6、数，该函数是基因型数据转换成便于分析的带有genetics包里最基本的函数。可以将以下四种形式的初始genotype class的数据。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. 数据框或矩阵中的两列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)实际的数据分析过程中建议将每一个snp 位点的基因型数据按照等位基因 /等位基因（e.g.a/a, a/t, t/t） 的格式给出，而不要简单

8、的用数字表示。这样有两个好处，一是可以很方便的用 makegenotypes() 函数一次性地将多个位点的原始基因型数据转换成带有 genotype 类属性的基因型数据，二是便于数据分析过程中了解具体是哪一个等位基因和疾病或性状有关。如果用数字的话，位点数目一多，可能就完全糊涂了。举个实例演示一下：library(genetics)#用 scan() 函数读入16 个人的数据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()函数将g1 中的两列基因型数据附上genotype类属性g2 - makegenotypes(g1)#大功告成，可以用str() 和 summary() 看看 g1 和 g2 的区别str(g1); str(g2)summary(g1); summary(g2)获取多态位点的基本信息我用 dgcgenetics包里面的popn 数据为例子，介绍一下获取多态位点基本信息的函数。data(popn,package=dgcgenetics) #首先载入popn 数据h

11、ead(popn) # 该数据包含四个多态位点（ a, b, c, and d ）、性别（ sex ）、疾病状态（ affect ）及 id 号（ subject ）。summary(popn$a) #获取 a 位点的基本信息，包括该位点分型成功率（call rate ）、等位基因频率、基因型频率、杂合度和多态信息含量（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定律。该定律是由英国数学家哈迪（d.h.hardy ）和德国医生温伯格（ w.weinberg）于 1908 年分别独立发现的，也称遗传平衡定律（geneticequilibrium law ）。该定律可以简单描述为，遗传平衡群体的等位基因频率与基

13、因型频率在世代间维持恒定。该定律的适用条件是：随机婚配，群体足够大，没有突变、选择、迁移和遗传漂变。在关联研究中hardy-weinberg平衡检验常被用来评价基因分型的质量。我们通常对病例和对照组分别进行hardy-weinberg平衡检验。如果某一位点在对照组中不符合hardy-weinberg平衡，我们通常会怀疑该位点的基因型鉴定的质量。如果该位点在对照组平衡而在病例组出现不平衡，则该位点很可能和疾病有关。genetics包里面提供两种不同的检验方法：一种是pearsons chi-square test，可以用hwe.chisq()函数进行该检验，另一种是fisher exact te

14、st，对应于hwe.exact()函数。hwe.chisq() 常用于 maf 较高、样本量较大的场合。 maf 较低的位点建议使用 hwe.exact() 函数。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 各种指标的函数，并能以文字和图形两种形式显示位点间的连锁不平衡程度。data(popn,package=dgcgenetics)# 首先载入popn 数据ldresult - ld(popn)#用 ld 函数计算位点间的ldsummary(ldresult, which=d ) # 用文字显示d值ldtable(ldresult, which =d ) # 用图形显示结果结果如下：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的一个函数，用在pool的 域。函数：haplo(n)用于生成n 个loci的haplotype configuration matrix;（一 ）所有可能haplotype00,0,1,1,0,1,1的矩 ，因 循 次数达到了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基因遗传学 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),