1、我接触 R 的时间算是不短了,已经两年多了。期间断断续续的看了些 R 网站上的材料。现在已经习惯了用 R 做数据分析了,并且越来越喜欢用 R 来做分析了。之前我用过 SAS,SPSS 也试过 Stata,但是这三个软件都没有专门的遗传统计模块(至少国内流行的盗版里没有)。所以和其它专业相比,我想 R 对我们也许更有用些。COS 论坛里提到 R 在 genetic statistics 里的应用的帖子很少。我在这里写一些我平时用到的遗传统计方面的 package 的说明,一来算是个人小结再者算是抛砖引玉吧,希望 COS 论坛里的各位多写些相关的东西。Introduction. CRAN Task
2、 View: Statistical GeneticsCRAN Task View 当中有一个单独的 Genetics 部分,里面列出了 40 个遗传统计相关的Package 和相关链接。这足可以看出 R 在遗传统计学当中的影响和作用。里面核心的 core package 有以下三个: genetics, gap, 和 haplo.stats。还有一个我经常用到的包是 DGCgenetics,算是对 genetics 包的扩展。以后我会提到以上几个包里面的一些函数。大致包括以下几方面的内容:1. 以上几个 package 对数据格式的要求;2. 多态位点的基本信息(MAF 等);3. Hard
3、y-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 三个核心
4、包及依赖的包。如果不加“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)遗传研究收集的数据有自己的特点。往往
5、是数据集中即包含一般的表型数据(分类和连续变量;如血压水平,BMI 和性别等),又包括基因型数据。分析时往往还需要用到不同的遗传模型,什么显性、隐形、加性模型,或者是按照分类变量来处理(有时候也称为共显性模型)。用SAS 或 SPSS 分析遗传数据时,如果要用不同的遗传模型进行数据分析的话,必须先进行数据转换,过程相对复杂。R 中的 genetics 包专门为基因型数据提供了一个新的 class(类),你可以很方便的用genotype()或 makeGenotypes()函数将不同形式的初始基因型数据转换成基因型数据,并为数据加上 genotype 类属性。genetics 包还提供了相应的
6、summary.genotype()和 plot.genotype()函数。你可以很方便的用 summary()函数获取基因型数据的基因型频率、等位基因频率等基本信息,用 plot()函数做出基因型的柱状图。先说一下 genotype()函数,该函数是 genetics 包里最基本的函数。可以将以下四种形式的初始基因型数据转换成便于分析的带有 genotype class 的数据。1. 以一个字符分隔的向量E.g. g1 #define DENOTvoid haplo(int *q,int *result)int i,j,tmp;/*int r=(2 in R, if q5 may cause
7、 flea*/int r=1;for (i=1;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,
8、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=sampl
9、e.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$
10、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),status=df2$censor)df2 - sample.pedsample.ped$ped=2,names(df2)df2$censor - c(1,1, rep(0, 12)ped2 - pedigree(df2$id, df2$fa
11、ther, 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),status=df2$censor)# create twin relationshipsrelate2 - matrix(c(210,211,1,212,213,3), nrow=2, byrow=TRUE)ped2 - pedigree(df2$id, df2$father, df2$mother,df2$sex, affected=as.matrix(aff2),status=df2$censor,relation=relate2)现在用来做这方面分析的用的比较多的是 GenABEL 这个包