使用bioconductor分析及处理Affymetrix SNP芯片

from Benilton Carvalho, Rafael A. Irizarry, Robert B. Scharpf, Vincent J. Carey. Processing and Analyzing Affymetrix SNP Chips with Bioconductor. Stat Biosci(2009)1:160-180

SNP芯片分析的基础知识
SNPs及Affymetrix SNP Chips
单核苷酸多态性(英语:Single Nucleotide Polymorphism,简称SNP,读作/snip/), 是在同一物种(species)内,不同个体之间或者成对染色体之间发生的单个核甘酸的变化。也就是说对于同一物种,其基因大体是相同的,但是对于个体之间,可能会存在少量碱基的变化。比如一段DNA片段,大多数个体可能是ACTGCG,但是会有少数个体会是ATTGCG。在这个例子中,它们发生了C和T的变化,被称为等位基因(alleles)。在两个等位基因中,大多数个体拥有的那个被称为major allele(主等位基因,因为大多数时候它是显性的,会被称为显性基因),而少数个体拥有的那个被称为minor allele(次等位基因,因为大多数时候它是隐性的,会被称为隐性基因)。

我们假设有一对等位基因,A和B,那么它在人的基因组中可能的组合是AA, AB, BB三种基因型。如果我它与某种表型(比如某种疾病)相关,人们对它所做的统计学分析就会被称为相关性研究(association studies)。但是往往一个SNP并不会有明显的表型,随着高通量技术的日益发展,相关性研究已经不再局限于单个的SNP了,而是多个生物标记(biomarkers)。而这种高通量技术在相关性研究上的应用就被称为Genome Wide Association Studies(GWAS, 基因组范围相关性研究)。从2003年一块SNP芯片上整合一万个生物标记到Affymetrix Genome-Wide Human SNP Array 6.0 (2007年)整合了90万(906600)个生物标记,在通量上已经有了长足的进步。

SNP芯片设计
和其它的基因芯片设计一样,SNP的探针沿用了25bp长,perfect-match(PM)及missmatch(MM)组合,及探针碱基错位爬行的策略来设计。举个实际的例子:

Strand Allele Offset Sequence[PM/MM] [SNP in bold]
Sense A -4 CGGGACCCTTCA[C/G]GTAAGTCACCTT
-3 ACGGGACCCTTC[A/T]GGTAAGTCACCT
-1 GAACGGGACCCT[A/T]CAGGTAAGTCAC
1 AGGAACGGGACC[C/G]TTCAGGTAAGTC
2 AAGGAACGGGAC[C/G]CTTCAGGTAAGT
C -4 CGGGACCCGTCA[C/G]GTAAGTCACCTT
-3 ACGGGACCCGTC[A/T]GGTAAGTCACCT
-1 GAACGGGACCCG[A/T]CAGGTAAGTCAC
1 AGGAACGGGACC[C/G]GTCAGGTAAGTC
2 AAGGAACGGGAC[C/G]CGTCAGGTAAGT
Antisense A -3 ACGGGACCCTTC[A/T]GGTAAGTCAGGT
0 GGAACGGGACCC[A/T]TCAGGTAAGTCA
1 AGGAACGGGACC[C/G]TTCAGGTAAGTC
2 AAGGAACGGGAC[C/G]CTTCAGGTAAGT
4 AGAAGGAACGGG[A/T]CCCTTCAGGTAA
C -3 ACGGGACCCGTC[A/T]GGTAAGTCAGGT
0 GGAACGGGACCC[A/G]TCAGGTAAGTCA
1 AGGAACGGGACC[C/G]GTCAGGTAAGTC
2 AAGGAACGGGAC[C/G]CGTCAGGTAAGT
4 AGAAGGAACGGG[A/T]CCCGTCAGGTAA

需要注意的是,SNP 6.0并不象之前的设计那样区分为PM/MM,而是只包含了PM。而且,对于同一个等位基因,它可能只有一个重复四次的探针,或者三个重复四次的探针。

SNP处理的理论知识
针对SNP芯片算法的发展历史
SNP分析的目的就是能够定量地区分AA, AB,以及BB三种基因型。我们假设\theta_A\theta_B分别是一对等位基因A和B的对数范围内的表达值。假设基因型为AA的话,那么理论上 \theta_A > C_A, \theta_B < C_B, 其中C_A, C_B是理论上相等的A和B同时被PCR扩增后应该被检测的信号。然而实际情况却并不那么简单。

从上图我们可以看出,即使是平行的三组实验,C_A, C_B的值也不能够统一,所以理论上的简单算法是不可行的。

之前介绍的如何处理基因芯片一般都分三步,背景较正(background adjustment), 归一化(标准化)(normalization)和总结(summarization)。实践表明,背景较正对于SNP芯片处理来说并没有太大的意义。因为对于高通量的SNP芯片来说,并不是平行实验中的对应的探针都能够足够强的信号,尤其是对于次等位基因来说。为此,Affymetrix改变了它的传统处理芯片的方式,转向使用动态模型(dynamic model,DM)算法。在这一算法当中,有四个高斯模型,NULL, AA, AB, 以及BB,对每一个SNP都用套用这四个模型,看属于那个分布。需要注意的是DM算法是分别对每一个芯片单独运算的。

但是DM算法并没有带来比传统默认算法更好的结果,还发现它会带来很多误判。为此,很多实验室纷纷开始发展自己的算法。其中robust linear model with Mahalanobis distance(RLMM)算法表现就超越了DM算法。RLMM算法是基于RMA的,它先用RMA做标准化(并没有背景较正步骤),然后使用a supervised learning algorithm similar to linear discriminate analysis来对每个SNP做出总结。线性判别分析(Linear Discriminant Analysis),简称判别分析,是统计学上的一种分析方法,用于在已知的分类之下遇到有新的样本时,选定一个判别标准,以判定如何将新样本放置于哪一个类别之中。它非常类似ANOVA,性线回归以及PCA分析。虽然它的表现远优于DM算法,但是,RLMM算法的一个缺点是它默认芯片之间的差异远大于芯片内部的差异,但实际上芯片内部的差异往往不能忽略。

2006年,Affymetrix公司基于RLMM算法提出了一种改进算法BRLMM(Bayesian Robust Linear with Mahalanobis distance)。这一算法是在RMA的基础上又加上了Bayesian算法。

Irizarry等在RMA的基础上于2007年推出了针对SNP芯片的SNPRMA算法。而这一算法如同他们推出的RMA算法一样,很快成为主流。

SNPRMA算法
SNPRMA算法的步骤如下,

  1. 对于每一个芯片,估算出探针碱基倾向性,SNP所在的位置,及长度带来的影响,以便以后的步骤中移除
  2. 使用四分位标准化(quantile normalization)方法消除芯片与芯片之间的差异
  3. 对于每个SNP,每对alleles,以及每对双链,应用线性模型对PM探针信号做出总结
  4. 对于每条链(strand),对alleles A和B求出信号强度,然后去除第一步中得到的碱基倾向性,SNP所在的位置及长度带来的影响,以及总信号强度的影响(total intensity effect)。
  5. 而后对于每个SNP的平均值及变化使用empirical Bayes或者最大似然比(maximizing the likelihood)进行估计。


更多细节请参见:http://biostatistics.oxfordjournals.org/content/8/2/485.long

分析实例
数据GSE7822拥有了Affymetrix大部分的SNP芯片的数据,包括50K XBA, 50K HIND, 250K STY等等。对于SNP5.0以及SNP6.0操作基本一致。
library("oligo")
library("affxparser")
##调用需要的注释包
annodb <- c('pd.mapping250k.sty', 'pd.mapping250k.nsp', 'pd.mapping50k.xba240', 'pd.mapping50k.hind240')
lapply(annodb, library, character.only=TRUE)
setwd("path/to/cel/dataset/")
##对文件进行分组,因不同类型的SNP文件不能同时进行分析
chiptype <- lapply(fullFilenames, function(.ele){
	header <- readCelHeader(.ele)
	c(header$filename, header$chiptype)
})
chiptype <- do.call(rbind, chiptype)
chiptype <- split(chiptype[, 1], chiptype[, 2])
##准备接受分析好的数据
crlmmOut <- list()
##对每组SNP文件,使用crlmm算法进行分析
for(i in 1:length(chiptype)){
	outputDir <- file.path(getwd(), names(chiptype)[i])
	if (!file.exists(outputDir)) crlmm(chiptype[[i]], outputDir)
	crlmmOut[[i]] <- getCrlmmSummaries(outputDir)
}
names(crlmmOut) <- names(chiptype)
##以下以mapping250k sty芯片为例进行进一步分析
sty <- crlmmOut[[2]]
styCalls <- calls(sty) - 1L ##这里使用calls调取SNP的类型,1为AA,2为AB,3为BB,为了方便,减去1,变成0:AA,1:AB,2:BB
styCallsGermline <- styCalls[,1:12]
styCallsMelanoma <- styCalls[,13:37]
styCountsGermline <- rbind(CC=rowSums(styCallsGermline==0L), CG=rowSums(styCallsGermline==1L), GG=rowSums(styCallsGermline==2L))
styCountsMelanoma <- rbind(CC=rowSums(styCallsMelanoma==0L), CG=rowSums(styCallsMelanoma==1L), GG=rowSums(styCallsMelanoma==2L))
styCountsGermline[, "SNP_A-2264525"]
styCountsMelanoma[, "SNP_A-2264525"]
##至此,我们得到了SNP的类型数据。接下来的工作就是读取注释
con <- pd.mapping250k.sty@getdb()
##显示注释包中所有的表
dbListTables(con)
##使用SQL语句查找需要的注释,以下是两个例子。
dbGetQuery(con, "select * from sequence limit 3")
dbGetQuery(con, "select * from featureSet where gene_assoc LIKE '%miR146a%'")

One thought on “使用bioconductor分析及处理Affymetrix SNP芯片

  1. 请问芯片数据的导入和注释成txt的语言是什么,上面的由于平台太多了所以没有太看懂,现在我就是在做一个pd.mapping250k.nsp芯片的数据,note包能够成功导入了,但是不知道怎么导入cel。尝试过使用readaffy得到的结果好像不太正确(反正用head看不到结果),而且之后也不能够用rma进行注释成txt,希望你能够帮我解答一下,感激不尽!

    > Data b<-rma(Data)
    Error in (function (classes, fdef, mtable) :
    unable to find an inherited method for function ‘rma’ for signature ‘"AffyBatch"’

hzg进行回复 取消回复

电子邮件地址不会被公开。 必填项已用*标注