R语言的遗传模块.docx

我接触R的时间算是不短了，已经两年多了。期间断断续续的看了些R网站上的材料。现在已经习惯了用R做数据分析了，并且越来越喜欢用R来做分析了。之前我用过SAS，SPSS也试过Stata，但是这三个软件都没有专门的遗传统计模块（至少国内流行的盗版里没有）。所以和其它专业相比，我想R对我们也许更有用些。COS论坛里提到R在genetic statistics里的应用的帖子很少。我在这里写一些我平时用到的遗传统计方面的package的说明，一来算是个人小结再者算是抛砖引玉吧，希望COS论坛里的各位多写些相关的东西。Introduction. CRAN Task View: Statistical GeneticsCRAN Task View当中有一个单独的Genetics部分，里面列出了40个遗传统计相关的Package和相关链接。这足可以看出R在遗传统计学当中的影响和作用。里面核心的core package有以下三个: genetics, gap, 和haplo.stats。还有一个我经常用到的包是DGCgenetics，算是对genetics包的扩展。以后我会提到以上几个包里面的一些函数。大致包括以下几方面的内容：1. 以上几个package对数据格式的要求；2. 多态位点的基本信息（MAF等）；3. Hardy-Weinberg平衡检验；4. LD的计算；5. 关联研究常用检验方法；6. Power的计算；先说一下前面提到的几个包的安装吧，其实很简单。一个一个用install.packages()函数来安装当然可以。相对简单点的方法是用CRAN Task Views里提到的ctv包来批量安装。install.packages(ctv) #首先安装ctv packagelibrary(ctv) #载入ctv packageinstall.views(Genetics,coreOnly = TRUE) #安装genetics, gap, haplo.stats三个核心包及依赖的包。如果不加core.only=TRUE则会安装所有的40个遗传统计相关的package。install.packages (genetics, coreonly = TRUE)DGCgenetics包的下载地址是http:/www-gene.cimr.cam.ac.uk/clayton/software/DGCgenetics_1.0.zip。你需要先下载这个包，然后本地安装。方法大家应该都知道，Rgui的Packages菜单的Install package(s) from local zip files。数据格式（1）遗传研究收集的数据有自己的特点。往往是数据集中即包含一般的表型数据（分类和连续变量；如血压水平，BMI和性别等），又包括基因型数据。分析时往往还需要用到不同的遗传模型，什么显性、隐形、加性模型，或者是按照分类变量来处理（有时候也称为共显性模型）。用SAS或SPSS分析遗传数据时，如果要用不同的遗传模型进行数据分析的话，必须先进行数据转换，过程相对复杂。R中的genetics包专门为基因型数据提供了一个新的class（类），你可以很方便的用genotype()或makeGenotypes()函数将不同形式的初始基因型数据转换成基因型数据，并为数据加上genotype类属性。genetics包还提供了相应的summary.genotype()和plot.genotype()函数。你可以很方便的用summary()函数获取基因型数据的基因型频率、等位基因频率等基本信息，用plot()函数做出基因型的柱状图。先说一下genotype()函数，该函数是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,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） 的格式给出，而不要简单的用数字表示。这样有两个好处，一是可以很方便的用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/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$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数据head(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 Proportion1 1786 0.62 1192 0.4NA 94 NAGenotype Frequency: #基因型频率 Count Proportion1/2 704 0.472/2 244 0.161/1 541 0.36NA 47 NAHeterozygosity (Hu) =0.4802686#杂合度Poly. Inf. Content =0.3648558 #PICHardy-Weinberg平衡检验首先简单介绍一下Hardy-Weinberg定律。该定律是由英国数学家哈迪（D.H.Hardy）和德国医生温伯格（W.Weinberg）于1908年分别独立发现的，也称遗传平衡定律（genetic equilibrium law）。该定律可以简单描述为，遗传平衡群体的等位基因频率与基因型频率在世代间维持恒定。该定律的适用条件是：随机婚配，群体足够大，没有突变、选择、迁移和遗传漂变。在关联研究中Hardy-Weinberg平衡检验常被用来评价基因分型的质量。我们通常对病例和对照组分别进行Hardy-Weinberg平衡检验。如果某一位点在对照组中不符合Hardy-Weinberg平衡，我们通常会怀疑该位点的基因型鉴定的质量。如果该位点在对照组平衡而在病例组出现不平衡，则该位点很可能和疾病有关。genetics包里面提供两种不同的检验方法：一种是Pearsons chi-square test，可以用HWE.chisq()函数进行该检验，另一种是Fisher exact test，对应于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$Acase)HWE.exact(popn$Acase)LD 的计算连锁不平衡是指人群中两个位点处在同一个单体型的频率比期望值高。评价连锁不平衡程度的指标包括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/trtrtdA 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;简单的说，就是形如0，0,0,1,1,0,1,1这样（一阶）所有可能haplotype的矩阵，因为循环次数达到了O(n*2n),所以用C语言写的，用R调用。附件中.so文件是Linux版本，.dll是Windows版本的（今天没有权限再上传附件了把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;for (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 - pedigree(id=sample.ped$id, dadid=sample.ped$father, momid=sample.ped$mother, sex=sample.ped$sex, famid=sample.ped$ped)print(pedAll)ped1basic - pedAll1ped2basic - pedAll2print(ped1basic)print(ped2basic)plot(ped2basic)# plot(ped1basic)kin2 - kinship(ped2basic)kin2pedAll - pedigree(id=sample.ped$id, dadid=sample.ped$father, momid=sample.ped$mother, sex=sample.ped$sex, famid=sample.ped$ped)kinAll - kinship(pedAll)kinAll1:14,1:14kinAll40:43,40:43kinAll42: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, affect