R语言的遗传模块.docx

资源描述

R语言的遗传模块.docx

《R语言的遗传模块.docx》由会员分享，可在线阅读，更多相关《R语言的遗传模块.docx（10页珍藏版）》请在冰豆网上搜索。

R语言的遗传模块.docx

R语言的遗传模块

我接触R的时间算是不短了，已经两年多了。

期间断断续续的看了些R网站上的材料。

现在已经习惯了用R做数据分析了，并且越来越喜欢用R来做分析了。

之前我用过SAS，SPSS也试过Stata，但是这三个软件都没有专门的遗传统计模块（至少国内流行的盗版里没有）。

所以和其它专业相比，我想R对我们也许更有用些。

COS论坛里提到R在geneticstatistics里的应用的帖子很少。

我在这里写一些我平时用到的遗传统计方面的package的说明，一来算是个人小结再者算是抛砖引玉吧，希望COS论坛里的各位多写些相关的东西。

Introduction.CRANTaskView:

StatisticalGenetics

CRANTaskView当中有一个单独的Genetics部分，里面列出了40个遗传统计相关的Package和相关链接。

这足可以看出R在遗传统计学当中的影响和作用。

里面核心的corepackage有以下三个:

genetics,gap,和haplo.stats。

还有一个我经常用到的包是DGCgenetics，算是对genetics包的扩展。

以后我会提到以上几个包里面的一些函数。

大致包括以下几方面的内容：

1.以上几个package对数据格式的要求；

2.多态位点的基本信息（MAF等）；

3.Hardy-Weinberg平衡检验；

4.LD的计算；

5.关联研究常用检验方法；

6.Power的计算；

…

先说一下前面提到的几个包的安装吧，其实很简单。

一个一个用install.packages（）函数来安装当然可以。

相对简单点的方法是用CRANTaskViews里提到的ctv包来批量安装。

install.packages（"ctv"）#首先安装"ctv"package

library（ctv）#载入ctvpackage

install.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菜单的Installpackage（s）fromlocalzipfiles。

数据格式

（1）

遗传研究收集的数据有自己的特点。

往往是数据集中即包含一般的表型数据（分类和连续变量；如血压水平，BMI和性别等），又包括基因型数据。

分析时往往还需要用到不同的遗传模型，什么显性、隐形、加性模型，或者是按照分类变量来处理（有时候也称为共显性模型）。

用SAS或SPSS分析遗传数据时，如果要用不同的遗传模型进行数据分析的话，必须先进行数据转换，过程相对复杂。

R中的genetics包专门为基因型数据提供了一个新的class（类），你可以很方便的用genotype（）或makeGenotypes（）函数将不同形式的初始基因型数据转换成基因型数据，并为数据加上genotype类属性。

genetics包还提供了相应的summary.genotype（）和plot.genotype（）函数。

你可以很方便的用summary（）函数获取基因型数据的基因型频率、等位基因频率等基本信息，用plot（）函数做出基因型的柱状图。

先说一下genotype（）函数，该函数是genetics包里最基本的函数。

可以将以下四种形式的初始基因型数据转换成便于分析的带有genotypeclass的数据。

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后分成两个allele

3.两个分开的向量

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）

获取多态位点的基本信息

我用DGCgenetics包里面的popn数据为例子，介绍一下获取多态位点基本信息的函数。

data（popn,package="DGCgenetics"）#首先载入popn数据

head（popn）#该数据包含四个多态位点（A,B,C,andD）、性别（sex）、疾病状态（affect）及ID号（subject）。

summary（popn$A）#获取A位点的基本信息，包括该位点分型成功率（callrate）、等位基因频率、基因型频率、杂合度和多态信息含量（PIC）

Numberofsamplestyped:

1489（96.9%）#callrate

AlleleFrequency:

（2alleles）#等位基因频率

CountProportion

1 1786 0.6

2 1192 0.4

NA 94 NA

GenotypeFrequency:

#基因型频率

CountProportion

1/2 704 0.47

2/2 244 0.16

1/1 541 0.36

NA 47 NA

Heterozygosity（Hu）= 0.4802686 #杂合度

Poly.Inf.Content = 0.3648558 #PIC

Hardy-Weinberg平衡检验

首先简单介绍一下Hardy-Weinberg定律。

该定律是由英国数学家哈迪（D.H.Hardy）和德国医生温伯格（W.Weinberg）于1908年分别独立发现的，也称遗传平衡定律（geneticequilibriumlaw）。

该定律可以简单描述为，遗传平衡群体的等位基因频率与基因型频率在世代间维持恒定。

该定律的适用条件是：

随机婚配，群体足够大，没有突变、选择、迁移和遗传漂变。

在关联研究中Hardy-Weinberg平衡检验常被用来评价基因分型的质量。

我们通常对病例和对照组分别进行Hardy-Weinberg平衡检验。

如果某一位点在对照组中不符合Hardy-Weinberg平衡，我们通常会怀疑该位点的基因型鉴定的质量。

如果该位点在对照组平衡而在病例组出现不平衡，则该位点很可能和疾病有关。

genetics包里面提供两种不同的检验方法：

一种是Pearson'schi-squaretest，可以用HWE.chisq（）函数进行该检验，另一种是Fisherexacttest，对应于HWE.exact（）函数。

HWE.chisq（）常用于MAF较高、样本量较大的场合。

MAF较低的位点建议使用HWE.exact（）函数。

library（genetics）

data（popn,package="DGCgenetics"）

control<-popn$affected=="Control"

case<-popn$affected=="Case"

HWE.chisq（popn$A[control]）

HWE.exact（popn$A[control]）

HWE.chisq（popn$A[case]）

HWE.exact（popn$A[case]）

LD的计算

连锁不平衡是指人群中两个位点处在同一个单体型的频率比期望值高。

评价连锁不平衡程度的指标包括D'、r2 等。

genetics包提供计算LD各种指标的函数，并能以文字和图形两种形式显示位点间的连锁不平衡程度。

data（popn,package="DGCgenetics"）#首先载入popn数据

ldresult<-LD（popn）#用LD函数计算位点间的LD

summary（ldresult,which="D’"）#用文字显示D’值

LDtable（ldresult,which="D’"）#用图形显示结果

结果如下：

PairwiseLD

———–

[table=50%][tr][td]

[/td][td]B

[/td][td]C

[/td][td]D

[/td][/tr][tr][td]AD'

[/td][td]0.978

[/td][td]0.976

[/td][/tr][tr][td]BD'

[/td][td]

[/td][td]0.998[/td][td]0.991

[/td][/tr][tr][td]CD'

[/td][td]

[/td][td]0.997

[/td][/tr][/table]

顺便帮楼主完成Haplo的一个函数，用在pool的领域。

函数：

haplo（n）用于生成n个loci的haplotypeconfigurationmatrix;简单的说，就是形如{｛0，0｝,｛0,1｝,｛1,0｝,｛1,1｝}这样（一阶）所有可能haplotype的矩阵，因为循环次数达到了O（n*2^n）,所以用C语言写的，用R调用。

附件中.so文件是Linux版本，.dll是Windows版本的（今天没有权限再上传附件了…把C代码附加在最后）。

R调用的代码如下：

dyn.load（"~/code/haplo.so"）#用者自酌

haplo<-function（n）{

densa<-.C（"haplo",as.integer（n）,result=as.integer（vector（"integer",n*2^n）））

densa[["result"]]

}

C的代码如下：

#include

#defineDENOT

voidhaplo（int*q,int*result）{

inti,j,tmp;

/*intr=（2<<（*q））;cannotuse<>inR,ifq>5maycauseflea*/

intr=1;

for（i=1;i<=*q;i++）r*=2;

for（i=0;i<（*q）;i++）{

for（j=0;j

result[i*r+j]=0;

}

for（j=0;j

tmp=j;

for（i=0;i<（*q）;i++）{

if（tmp>=1）{

result[（*q-i-1）*r+j]=tmp%2;

tmp/=2;

}

#ifndefDENOT

for（i=0;i<（*q）;i++）{

printf（"\n"）;

for（j=0;j

printf（"%d\t",result[i*r+j]）;

}

#endif

}

kinship2基因遗传学[s:

11]

library（kinship2）

data（sample.ped）

sample.ped[1:

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<-pedAll['1']

ped2basic<-pedAll['2']

print（ped1basic）

print（ped2basic）

plot（ped2basic）

#plot（ped1basic）

kin2<-kinship（ped2basic）

kin2

pedAll<-pedigree（id=sample.ped$id,

dadid=sample.ped$father,momid=sample.ped$mother,

sex=sample.ped$sex,famid=sample.ped$ped）

kinAll<-kinship（pedAll）

kinAll[1:

14,1:

14]

kinAll[40:

43,40:

43]

kinAll[42:

46,42:

46]

df2<-sample.ped[sample.ped$ped==2,]

names（df2）

df2$censor<-c（1,1,rep（0,12））