基因芯片差异性分析

, 01 Nov 2019



起初使用的是oligo与limma结合分析，可惜的是rma校正后limma分析的结果不合实际，所以直接使用GEO2R上的代码，并修改成不在线的读取data.frame。但是仍旧得探究limma直接读取cel文件的情况。其中还遇到了读取CDF文件的情况。

芯片分析步骤

基因芯片的差异表达分析主要有 构建基因表达矩阵、构建实验设计矩阵、构建对比模型（对比矩阵）、线性模型拟合、贝叶斯检验和生成结果报表 六个关键步骤。

文章使用的包多半是 bioconductor 上的包

芯片数据下载

boxplot(exset_rma, value = "pm")

到这基本读取完毕，但是稍后使用这个矩阵得到的结果与GEO2R的结果大相径庭，其错误目前只能怀疑这个oligo包不适用这个芯片，不适用我研究过的五六个芯片。。。有时间需要查看这个包。

在这里提个CDF文件的读取

CDF文件读取

提交者产生CEL文件的时候使用了特殊的CDF文件。所以要用特殊化的方法进行分析。 用affy包和特殊的CDF文件分析,报错The affy package is not designed for this array type, Please use either the oligo or xps package这是因为生成了CDF的custom环境。 custom文件的官网 我使用的特殊文件需使用affy\==1.40.1版本，而这一版本又依赖于R\==3.4版本以下

Version 18 (Data Sources)
Jan 23, 2014, version 18 is released. You might want download our modified affy package, if you got error message "The affy package is not designed for this array type, Please use either the oligo or xps package."
affy_1.40.0.tar.gz or affy_1.40.1.zip

custom使用脚本

#读取CDF文件：
library(makecdfenv)
#载入需要的程辑包：affyio
pkgpath <-'/data3/liaorp/1'
make.cdf.package("hugene20st_Hs_ENTREZG.cdf",cdf.path= '/data3/liaorp/1', compress=FALSE, species = "Homo_sapiens", package.path = pkgpath)
dir(package)
library(hugene20sthsentrezgcdf ) #导入构建好的包
library(affy) #此处affy为1.14.0，寻找CDF文件对应affy包
celFiles <- list.celfiles("/data3/liaorp/MICROARRAY/测试芯片", full.name=TRUE)
celFiles
raw.data=ReadAffy(verbose=TRUE, filenames=celFiles, cdfname="hugene20sthsentrezgcdf")

以上为记录，我对其也不甚了解

使用limma包进行分析

有必要研究limma这个包如何读取原始文件

#使用下载好的基因表达矩阵读取，
gset <- read.table(file = "GSE19332_series_matrix.txt", header =TRUE,
                      comment.char = "!", row.names=1)

if(T){
  ex <- gset[,] #此处可以挑选感兴趣的组
  write.csv(ex, file="yourname.csv", row.names=T)#导出原始数据表格
  # 求log2值
  qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
  LogC <- (qx[5] > 100) || (qx[6]-qx[1] > 50 && qx[2] > 0) || (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
  if (LogC) { ex[which(ex <= 0)] <- NaN
  expmatx <- log2(ex) } #得到log2转换的矩阵
#此处就是GEO2R的代码

构建分组矩阵

suppressMessages(library(limma))
grouplist = c(rep("Control",7), rep("ALS",3))#复制矢量
design <- model.matrix(~0+factor(grouplist)) #创建design
colnames(design)=levels(factor(grouplist)) # 实验设计矩阵的每一行对应一个样品的编码
rownames(design)=colnames(expmatx) # 将行名换成表达矩阵的列名
design

记住，上文的grouplist这个向量应该与表达矩阵的列名是一一对应的，即是ALS对应的是表达矩阵后面的三个列，即样品名

构建比对矩阵

contrast.matrix<-makeContrasts(paste0(unique(grouplist),collapse = "-"),levels = design)
contrast.matrix

线性模型拟合及贝叶斯分析

fit <- lmFit(expmatx,design)
fit2 <- contrasts.fit(fit, contrast.matrix) 
fit2 <- eBayes(fit2,0.01)
tempOutput = topTable(fit2,coef=1,n=Inf,adjust="fdr", sort.by="B", number=1000) #此处的coef= int 对应比对矩阵的列,然后列出1000个基因或者直接全部列出呗
nrDEG = na.omit(tempOutput) # 去除缺失值的行
write.csv(nrDEG,"ALSproject-GSE19332.csv",quote=F)

以上即为GSE19332差异分析的代码

CLL包自带的数据

AffyBatch类是从一个更基础的类eSet类衍生来的。eSet类非常重要，它被写成了一个虚类，衍生出许多非常重要的类，包括ExpressionSet类，SnpSet类以及AffyBatch类等。eSet是Bioconductor为基因表达数据格式所定制的标准，因此非常有必要熟悉eSet及其衍生类。

library(CLL)
library(affy)
data(CLLbatch)
class(CLLbatch)#可以看到CLLbatch为AffyBatch类
library(CLL)
# rma方法进行背景校正【当MM值比PM值还要高时，MM就是杂信号，也就是背景噪声，需要去除】
CLLrma <- rma(CLLbatch)
e_before <- exprs(CLLbatch)
e_after <- exprs(CLLrma)
#对比一下校正前后数据
e_before[1:5,1:5]
e_after[1:5,1:5]
#背景校正，标准化，汇总可以用一个函数来进行处理（三合一函数），但是有文献报道，该方法有缺点！暂时没有着手看这个函数
eset.mas<-expresso(afbatch = CLLbatch,bgcorrect.method = "mas",
                   normalize.method = "constant",pmcorrect.method ="mas",
                   summary.method = "mas")

参数详解：

查看一下各种参数的方法种类：

bgcorrect.methods() normalize.methods(CLLbatch) pmcorrect.methods() express.summary.stat.methods()

help(“AffyBatch”)#查看帮助信息，似乎没有找到rma()作用于AffyBatch这个对象 #但是你会发现可以直接使用boxplot()直接对其绘图

这个数据集有空看看。

数据后续分析

pathway代谢通路:通路分析最常使用的是KEGG【日本主导的一个项目对gene和genome进行了非常详细的注释】

KEGG:Kyoto Encyclopedia of Genes and Genomes: 系统分析基因产物和化合物在细胞中的代谢途径以及这些基因产物的功能的数据库【基于ORA算法】

GO负责分门别类，而pathway负责把每一类对应到具体的代谢网络中。研究pathway的原因是：生物学问题中设定一个“蝴蝶效应”假设：1个Pathway上游基因的改变，会导致下游相关基因改变，从而改变通路中大量基因的表达。

富集分析：把上游的这些差异基因再进行注释、分组，一个类别就相当于一个GO term，然后看这几大类的区别，肯定比看几十甚至上百个基因或蛋白的差异要更加直观。目的就是根据不同功能，把各个分子进行分类，然后使用超几何分布检验进行分析。其中涉及概念，前景基因，背景基因，重点研究的基因集叫做前景基因，需要比对的所有基因集叫背景基因，前景是背景的子集。例如差异性分析数据中的差异基因就是前景基因，所有的表达基因就是背景基因。

基因芯片：以Affymetrix为例，一个芯片可以包括上百万探针（通常由25个碱基组成），被整齐印刷在芯片上。 探针组：一个探针组通常由20对个或11个探针对组成，来自一个基因 探针对：匹配探针（Perfect match， PM）+错配探针（Mismatch， MM）。二者仅仅是中间的那个碱基不同。并非所有芯片都有这两个探针，比如外显子芯片，每个探针组只有4个PM探针，没有MM探针

CEL文件：CEL文件是Affymetrix最常用的格式，它的主要内容就是每个“cell”的灰度信息，而“cell”就是整个芯片图像划分后得到的小网格，每个小网格中的图像被看作来自一个探针。当然，如果要得到芯片上每个探针组对应的表达数据，还需要探针的排布信息（哪个探针对应哪个探针组），这部分信息就存储在CDF文件中。要想知道探针对应的序列信息，就需要用Probe文件 附上一个网址，我还没看过，https://slideplayer.com/slide/4804237/

富集分析算法：大体上富集分析有四类算法：ORA、FCS、PT、NT1 四种算法图 i)ORA(Over Representation Analysis)：过表达分析 它是检验某类功能在一个数据子集中是否表现过度。又称为“2X2方法”，像上图一样，做一个列联表。上图中的ORA中，蓝圈内是感兴趣基因（8个），绿圈内是某个通路的基因（5个）；灰点是既不感兴趣又不在通路内的（6个），蓝点是感兴趣但不在通路内的（5个），绿点是在通路内但不感兴趣的（2个），红点是既感兴趣又在通路内的（3个），于是就能做出来2X2列联表。再利用fisher精确检验或超几何分布得到p值。 简而言之，需要4类数据：总共的基因数（作为背景基因）、总共属于某分类的基因数、样本包含的基因数（也就是用的差异表达基因）、样本中属于某分类的基因数

优点：出现的最早，最常用，有完善的统计学理论基础，结果比较可靠；

缺点：仅仅使用了基因的数目，但是基因的不同表达水平没有考虑，为了得到差异基因，需要人为设置阈值，没有一个设置规定，因此结果因人而异；适用于差异最显著的基因，而差异不显著的基因就会被忽略，检测灵敏度会降低ORA利用统计学假设每个基因相互独立，但是就生物体本身而言，忽略了内部的复杂的相互作用，并且每个基因在不同的生物学过程中发挥的作用大小不一样，同等看待结果可能会不准确

实际上就是把我们感兴趣的基因和背景基因做一个交集。感兴趣的基因也就是差异基因了，包括上调、下调表达的（利用原始表达矩阵中p值和logFC进行筛选），一般人类芯片数据会有几百个;背景基因就是在KEGG等数据库中有注释的基因【人类基因组有2万个左右基因，现在总共有已知功能的是7000左右，随着研究的不断深入，背景基因数量会越来越多，结果也会越来越全面】 KEGG结果

举个例子，KEGG通路hsa05206指的是MicroRNAs in Cancer，包括150个基因，背景基因使用了6517个；GSE17708芯片得到的差异基因数是547个，在KEGG能注释上的有80个，其中就有10个是MicroRNA通路的，概率高达12.5%（enrichKEGG方法都是用能在KEGG注释上的基因，比如这里用80而不是547），那么这个通路是不是在下调基因中被显著改变？需要把全部的80个下调基因，在KEGG的530个通路中注释一遍，再一个一个进行超几何分布检验，得到p值。hsa05206通路在背景基因中查到的概率是150/6517=2.3%，是显著低于12.5%的 >超几何分布属于统计学上一种离散概率分布。它描述了由有限个物件中抽出n个物件，成功抽出指定种类的物件的个数（不归还）。n=1，超几何分布还原为伯努利分布；n接近∞，超几何分布可视为二项分布

差异分析R包

差异分析三R包：Limma、DESeq2、edgeR

limma

芯片数据普遍认为符合正态分布,而正态分布的群体一般就是用t检验（两个样本）或者方差分析（多个样本）,limma采用贝叶斯模型（Empirical Bayesian model），更新的limma-voom适配了转录组数据。

DESeq2

DESeq2：采用负二项分布算法（negative binomial distribution） >RNA-seq中，技术误差是满足泊松分布的，因为期望和方差差不多。但是生物学重复之间的误差不能用泊松分布来描述，因为他的方差可能很大，所以要用负二项分布2

edgeR

使用DEGList读取表达矩阵=》利用count-per-millionCPM严格过滤count值低的数据=〉calcNormFactors函数使用TMM算法对矩阵标准化=》实验设计矩阵Design matrix , model.matrix(~0+group) =〉估算离散值dispersionestimateDisp=》构建比较矩阵makeContrasts 、glmQLFTest=〉提取差异基因decideTestsDGE 、glmTreat

  1. Progress in Gene Functional Enrichment Analysis [return]
  2. https://www.cnblogs.com/leezx/p/6036806.html [return]