houliping
自我介绍一下,我是医学院的,现在博士在读,研究方向是遗传流行病。接触R完全是拜我研究生师兄所赐。我们科室有几个数学系的学生,老板招他们是想提高我们研究组的统计能力。R也是他们其中的一位介绍给我的。
我接触R的时间算是不短了,已经两年多了。期间断断续续的看了些R网站上的材料。现在已经习惯了用R做数据分析了,并且越来越喜欢用R来做分析了。之前我用过SAS,SPSS也试过Stata,但是这三个软件都没有专门的遗传统计模块(至少国内流行的盗版里没有)。所以和其它专业相比,我想R对我们也许更有用些。
COS论坛里提到R在genetic statistics里的应用的帖子很少。我在这里写一些我平时用到的遗传统计方面的package的说明,一来算是个人小结再者算是抛砖引玉吧,希望COS论坛里的各位多写些相关的东西。
Introduction. CRAN Task View: Statistical Genetics
CRAN 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的计算;
...
zwdbordeaux
加油!
ypchen
好东西
cloud_wei
热烈欢迎!!
houliping
先说一下前面提到的几个包的安装吧,其实很简单。一个一个用install.packages()函数来安装当然可以。相对简单点的方法是用CRAN Task Views里提到的ctv包来批量安装。
install.packages("ctv") #首先安装"ctv" package
library(ctv) #载入ctv package
install.views("Genetics",coreOnly = TRUE) #安装genetics, gap, haplo.stats三个核心包及依赖的包。如果不加"core.only=TRUE"则会安装所有的40个遗传统计相关的package。
DGCgenetics包的下载地址是
http://www-gene.cimr.cam.ac.uk/clayton/software/DGCgenetics_1.0.zip。你需要先下载这个包,然后本地安装。方法大家应该都知道,Rgui的Packages菜单的Install package(s) from local zip files。
pengchy
很好!我主要用R作芯片数据分析,目前还在学习阶段。从正式开始学习R到现在有1年半时间了。
oldbeggar
[quote]引用第5楼pengchy于2008-10-19 11:17发表的“”:
很好!我主要用R作芯片数据分析,目前还在学习阶段。从正式开始学习R到现在有1年半时间了。[/quote]
Bioconductor,做芯片分析简直无敌~~~
比Affymetrix的软件爽多了~
houliping
数据格式(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后分成两个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)
or
data1 <- cbind(allele1 = c("D","D","D","I",""), allele2 = c("D","I","D","I",""))
g6 <- genotype(data1)
asdnwd224
如何在R里面安装这些,这些组件。学习中谢谢2
houliping
数据格式(2)
实际的数据分析过程中建议将每一个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/C
2 39 2 24.5 T/T C/G
3 36 1 23 A/T C/C
4 41 2 26 A/T C/C
5 37 1 29.5 A/A C/C
6 35 2 31 A/T G/G
7 41 2 21.5 A/A C/G
8 43 2 27.5 A/A G/G
9 44 2 24 A/A C/C
10 32 1 19.5 A/T C/C
11 40 2 20 A/A C/G
12 38 2 22.5 T/T G/G
13 42 2 32.5 A/A C/C
14 33 2 25.5 A/T C/C
15 43 1 30.5 A/T G/G
16 35 2 25 A/T C/C
g1 <- 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)
pengchy
参看第4楼
[quote]引用第8楼asdnwd224于2008-10-29 22:52发表的“”:
如何在R里面安装这些,这些组件。学习中谢谢2[/quote]
houliping
获取多态位点的基本信息
我用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 ([u]96.9%[/u]) #call rate
Allele Frequency: (2 alleles) #等位基因频率
Count Proportion
1 1786 [u]0.6[/u]
2 1192 [u]0.4[/u]
NA 94 NA
Genotype Frequency: #基因型频率
Count Proportion
1/2 704 [u]0.47[/u]
2/2 244 [u]0.16[/u]
1/1 541 [u]0.36[/u]
NA 47 NA
Heterozygosity (Hu) = [u]0.4802686[/u] #杂合度
Poly. Inf. Content = [u]0.3648558[/u] #PIC
houliping
Hardy-Weinberg平衡检验
首先简单介绍一下Hardy-Weinberg定律。该定律是由英国数学家哈迪(D.H.Hardy)和德国医生温伯格(W.Weinberg)于1908年分别独立发现的,也称遗传平衡定律(genetic equilibrium law)。该定律可以简单描述为,遗传平衡群体的等位基因频率与基因型频率在世代间维持恒定。该定律的适用条件是:随机婚配,群体足够大,没有突变、选择、迁移和遗传漂变。
在关联研究中Hardy-Weinberg平衡检验常被用来评价基因分型的质量。我们通常对病例和对照组分别进行Hardy-Weinberg平衡检验。如果某一位点在对照组中不符合Hardy-Weinberg平衡,我们通常会怀疑该位点的基因型鉴定的质量。如果该位点在对照组平衡而在病例组出现不平衡,则该位点很可能和疾病有关。
genetics包里面提供两种不同的检验方法:一种是Pearson's chi-square test,可以用HWE.chisq()函数进行该检验,另一种是Fisher exact test,对应于HWE.exact()函数。HWE.chisq()常用于MAF较高、样本量较大的场合。MAF较低的位点建议使用HWE.exact()函数。
<br />
library(genetics)<br />
data(popn, package="DGCgenetics")<br />
control <- popn$affected == "Control"<br />
case <- popn$affected == "Case"<br />
HWE.chisq(popn$A[control])<br />
HWE.exact(popn$A[control])<br />
HWE.chisq(popn$A[case])<br />
HWE.exact(popn$A[case])<br />
macklon
houliping 大侠可以留个联系方式么?我也似医学院的,目前搞这个方向。初入道,硕士是搞动物的,现在博士一下子搞遗传,有点难啊。能帮帮我么啊?先谢谢了!
gelu0
lz辛苦。某小硕,遗传学背景,将来方向之一是统计遗传学。有空多交流。
以前做SNP关联的时候用过SNPassoc这个包,印象中还不错。
最近发现epicalc这个包也挺实用的,虽然这个包不是针对遗传学问题设计的。
houliping
LD 的计算
连锁不平衡是指人群中两个位点处在同一个单体型的频率比期望值高。评价连锁不平衡程度的指标包括D'、r[sup]2[/sup] 等。genetics 包提供计算LD 各种指标的函数,并能以文字和图形两种形式显示位点间的连锁不平衡程度。
<br />
data(popn,package="DGCgenetics") #首先载入popn数据<br />
ldresult <- LD(popn) #用LD函数计算位点间的LD<br />
summary(ldresult, which="D’") #用文字显示D’值<br />
LDtable(ldresult, which = "D’") #用图形显示结果<br />
结果如下:
Pairwise LD
-----------
[table=50%][tr][td]
[/td][td]B
[/td][td]C
[/td][td]D
[/td][/tr][tr][td]A D'
[/td][td]0.978
[/td][td]0.976
[/td][td]0.976
[/td][/tr][tr][td]B D'
[/td][td]
[/td][td]0.998[/td][td]0.991
[/td][/tr][tr][td]C D'
[/td][td]
[/td][td]
[/td][td]0.997
[/td][/tr][/table]
houliping
[quote]引用第14楼gelu0于2009-06-07 17:37发表的 :
lz辛苦。某小硕,遗传学背景,将来方向之一是统计遗传学。有空多交流。
以前做SNP关联的时候用过SNPassoc这个包,印象中还不错。
最近发现epicalc这个包也挺实用的,虽然这个包不是针对遗传学问题设计的。
[/quote]
谢谢,SNPassoc包可以用来做whole genome association study的数据分析,感兴趣的话还可以试试另外一个包--GenABEL。
epicalc是一个用于流行病研究的包,类似的包还有epitools,我个人更喜欢用后面那个。
moxingbo
好强呐!!!
enthumelon
顺便帮楼主完成Haplo的一个函数,用在pool的领域。函数:haplo(n)用于生成n个loci的haplotype configuration matrix;简单的说,就是形如{{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<R.h>
#define DENOT
void haplo(int *q,int *result){
int i,j,tmp;
/*int r=(2<<(*q));cannot use << or >> in R, if q>5 may cause flea*/
int r=1;
for (i=1;i<=*q;i++)r*=2;
for (i=0;i<(*q);i++){
for(j=0;j<r;j++){
result[i*r+j]=0;
}
}
for (j=0;j<r;j++){
tmp=j;
for (i=0;i<(*q);i++){
if (tmp>=1){
result[(*q-i-1)*r+j]=tmp%2;
tmp/=2;
}
}
}
#ifndef DENOT
for (i=0;i<(*q);i++){
printf("\n");
for(j=0;j<r;j++){
printf("%d\t",result[i*r+j]);
}
}
#endif
}
dingpeng
能否介绍一下遗传统计学中的统计问题呢?
似乎蛮有趣的。
尤其是epistasis,呵呵,最近看了一些相关的文献,看不太懂生物学家的语言。