高通量芯片数据分析：转录组芯片数据分析

利用R的bioconductor包进行分析。由于安装的是R3.5以上版本所以实际用的是用biomanager指令，其他基本一样。
不同的包有各类坑，具体可以查阅bioconductor官网寻找解决办法。R包安装踩过的各种bug大全见这篇草稿，没有什么是重装不能解决的，有就换环境（暴言）

一、初始配置

1. bioconductor包安装

if (!requireNamespace("BiocManager", quietly = TRUE))
     install.packages("BiocManager")  #3.5之后的新版本安装

2. 所需包的安装和加载

BiocManager::install()  
#BiocManager::install(c("annotate","affycoretools","graph","gcrma","CLL","simpleaffy","affyPLM","RColorBrewer","affy"))
library()

二、数据预处理

以CLL基因表达数据包（详见官方说明书）为例。

1.原始数据导入

library(CLL)              #载入包
data("CLLbatch")     #读入CLLbatch数据集
CLLrma <- rma(CLLbatch)   #这里调用RMA算法来对数据标准化
e <- exprs(CLLrma)  #读取预处理后所有样品基因的表达值矩阵

表达矩阵（局部）如下图所示，行对应Probeset（而不是基因），列代表不同sample。

补充说明：

• rma函数：
  ①CLL包自动调用affy包，affy包含有多种处理函数如rma函数；
  ②rma函数是一体化的标准化函数，芯片其他标准化方法见3
• 查看其他CLL数据集信息
  data.class(CLLbatch) #查看数据类型，发现是AffyBatch类(和ExpressionSet类，SnpSet类都属于eSet类）
  data(disease) 、disease #读入、查看所有样本的状态信息（稳定期、发展期）
  help(AffyBatch)、phenoData(CLLbatch)、featureData(CLLbatch)、protocolData(CLLbatch)、annotation(CLLbatch) #查看"AffyBatch"类的详细介绍
• 表达量矩阵调取：assayData槽与exprs方法
  assayData槽是AffyBatch类必不可少的，其方法返回结果的第一个元素是矩阵类型，用于保存基因表达矩阵。当使用exprs方法时，调取的就是这个基因表达矩阵。→ExpressionSet对象介绍

2.质量评估与控制
三个层次的质控：直观观察（image函数）、平均值方法（simpleaffy包）、数据拟合方法（affyPLM包）

2.1 image函数查看芯片的灰度图像
直接观察法比较粗糙（非量化），但可以对芯片质量产生总体认识。

image(CLLbatch[,1])     #对第一列（一个样本)画图

• 几个灰度图的判断指导：
  ①如果图像特别黑，说明型号强度低；如果图像特别亮，说明信号强度很有可能过饱和；
  ②如果芯片图像有斑块现象就很可能是坏片；
  ③如果四个点看不到印花，说明数据有问题。

2.2 用simpleaffy包对芯片数据进行质量评估
平均方法假设一组实验中的每个芯片数据对于某个平均值指标都相差不大（质量均匀）。

library(simpleaffy)
data(CLLbatch)
Data.qc <- qc(CLLbatch) # 读取芯片数据并获取质量分析报告
plot(Data.qc)  # 图型化显示报告

• 结果图解释
  ①第一列：样品名称。
  ②第二列：检出率、平均背景噪声（往往较高的平均背景噪声都伴随着较低的检出率）
  ③第三列：蓝色正常，红色异常
  尺度因子蓝色实心圆：取值（-3，3）。
  GAPDH 3’/5’ 比值“o”：不大于1.25，否则数据质量不好
  actin 3’/5’比值“△”：不大于3，否则数据质量不好。
  bioB：芯片检测没有达标。

2.3 用affyPLM包对芯片数据进行质量评估
affyPLM包可对芯片原始数据进行拟合回归，最后得到芯片权重（Weights）图、残差（Residuals）图、相对对数表达（RLE，Relative log expression）箱线图、相对标准差（NUSE，Normalized unscaled standard errors）箱线图等等,

library(affyPLM)
data(CLLbatch)     #读取数据
Pset=fitPLM(CLLbatch)  #对数据集做回归分析，结果为一个PLMset类型的对象Pset

image(Pset,type="weights",which=1, main="Weights") #根据计算结果，画权重图
image(Pset,type="resids",which=1, main="Residuals") #根据计算结果，画残差图
image(Pset,type="sign.resids",which=1, main="Residuals.sign") #根据计算结果，画残差符号图

2.4 总体质量评估：RLE和NUSE箱线图
接上得到PLMset结果后，一个样品的所有探针组的分布可用箱型图表示。

RLE：反映基因表达量的一致性趋势，定义为一个探针组在某个样品的表达值除以该探针组在所有样品中表达值的中位数后取对数。
NUSE：定义为一个探针组在某个样品的PM值的标准差除以该探针值在各样品中的PM标准差的中位数。

理论上，随机分布应该画出接近0均匀分布的权重图和接近1均匀分布的残差图。

library(RColorBrewer)  #RColorBrewer包用于提供配色方案。
colors<-brewer.pal(12,"Set3")  #载入一组颜色

#利用上一步的affyPLM包的回归计算结果Pset画箱线图
Mbox(Pset, ylim=c(-1,1),col=colors,main="RLE",las=3)# 绘制RLE(相对对数表达)箱线图
boxplot(Pset,ylim=c(0.95,1.22),col=colors,main="NUSE",las=3) #绘制NUSE（相对标准差）箱线图

CLL1，CLL10 的质量明显异于其他样品，考虑去除

2.5 查看RNA降解曲线

“因为RNA是从5’端开始降解的，理论上5‘端的荧光强度应该低于3’端的荧光强度。如果RNA降解曲线的斜率太小，甚至接近于0，很可能是RNA降解太严重，需要作为坏数据去除。”

library(affy)
# data("CLLbatch")
data.deg <- AffyRNAdeg(CLLbatch) # 获取芯片数据和降解数据
#colors<-brewer.pal(12,"Set3") 
plotAffyRNAdeg(data.deg, col=colors) #利用上面载入过的一组颜色画图
legend("topleft", rownames(pData(CLLbatch)), col=colors, lwd=1, inset=0.05, cex=0.5) #在左上部位加注图

RNA降解图找出了CLL13需要剔除。

综合上述所有数据，把质量差的3个样本的数据从数据集剔除。

CLLbatch <- CLLbatch[, -match(c("CLL10.CEL", "CLL1.CEL", "CLL13.CEL"), sampleNames(CLLbatch))]

2.5 从聚类分析的角度看数据质量
从芯片之间的相互关系来检验芯片的质量，可以画聚类分析树，或做PCA分析。

library(CLL)
library(gcrma)
library(graph)

data(CLLbatch)
data("disease") #数据加载

CLLgcrma<-gcrma(CLLbatch) #使用gcrma算法预处理数据
eset<-exprs(CLLgcrma) #提取基因表达矩阵
pearson_cor<-cor(eset)# 计算样品两两之间的Pearson相关系数
dist.lower<-as.dist(1-pearson_cor)# 得到Pearson距离的下三角矩阵

hc<-hclust(dist.lower,"ave")# 聚类分析并画图
plot(hc)

从聚类分析的结果来看，稳定组（红框）和恶化组分不开。再做主成分分析（PCA）：

library(affycoretools)
samplenames<-sub(pattern="\\.CEL", replacement ="",colnames(eset))  #提取样本名
groups<-factor(disease[,2])  #提取样本状态
plotPCA(eset,addtext=samplenames,groups=groups,groupnames=levels(groups))  #被PCA的是表达值矩阵

结果类似。两者分不开。

注意：
“①只有当聚类图中有明显的类别差异时，才适合考虑去除个别不归类的样品；如果类似CLL数据集来自不同个体差异的原因整体分类被打乱，则不能简单判定所有样品都出了问题。
②使用主成分分析时，还必须考虑前2个主成分是否具有代表性，这要看前2个主成分的累计贡献率，如果低于60%，可以考虑采用多维尺度分析来构建分类图。”

3.背景校正、标准化和汇总*
除了质量控制剔除不合格的样品，还需要关注下面的三步预处理层面：

背景校正background correction：去除背景噪声
标准化normalization：使各次/组测量或各种实验条件下的测量可以相互比较
汇总summarization：使用一定的统计方法将前面得到的荧光强度值从探针水平汇总到探针组水平

• 各类预处理方法

• 方法：
  ①非一体化方法
  有affy包的expresso函数、threestep函数。自定义参数的空间比较大，但是使用比较复杂，详见help()命令。
  ②一体化预处理算法（推荐）
  包括gcrma、rma、和mas5。有时需进行算法比较，进而确定哪种算法最合适

dmas5 <- mas5(data.raw)                        
drma <- rma(data.raw)                                 
dgcrma <- gcrma(data.raw)    

三、基因芯片数据分析

（一）选取差异表达基因

六个关键步骤详解见参考资料。

1.构建基因表达矩阵

library(gcrma)
library(CLL)
data("CLLbatch")
data(disease)      #载入所需包和数据

CLLbatch <- CLLbatch[, -match(c("CLL10.CEL", "CLL1.CEL", "CLL13.CEL"),sampleNames(CLLbatch))]  #移除质控环节决定剔除的数据
CLLgcrma <- gcrma(CLLbatch) #用gcRMA算法进行预处理
#################

sampleNames(CLLgcrma) <- gsub(".CEL$", "",sampleNames(CLLgcrma)) #remove .CEL in sample names
disease <- disease[match(sampleNames(CLLgcrma), disease[, "SampleID"]), ] # remove record in data disease about 异常数据

eset <- exprs(CLLgcrma) #获得余下21个样品的基因表达矩阵

2.构建实验设计矩阵

disease <- factor(disease[, "Disease"])#提取实验条件信息
design <- model.matrix(~-1+disease) #构建实验设计矩阵

3.构建对比模型（对比矩阵）

contrast.matrix <- makeContrasts(contrasts = "diseaseprogres. - diseasestable",levels=design)#构建对比模型比较两个实验条件下表达数据

4.线性模型拟合

fit <- lmFit(eset, design)

5.贝叶斯检验
使用limma包实现贝叶斯检验。limma包对输入数据要求是必须经过对数转换的表达值，所以前面预处理需要经过log2变换的gcRMA算法

 library(limma)
 fit1 <- contrasts.fit(fit, contrast.matrix) #根据对比模型进行差值计算 
 fit2 <- eBayes(fit1)

fit2我遇到一个Warning message:Z ero sample variances detected, have been offset away from zero
据说是因为至少有一个gene它的表达值在所有的样品中没有任何变化，这个gene会在计算中被忽略。不影响我们跑。

6.生成结果报表。

dif <- topTable(fit2, coef="diseaseprogres. - diseasestable", n=nrow(fit2), lfc=log2(1.5))#生成所有基因的检验结果报告
dif <- dif[dif[, "P.Value"]<0.01,] #用P.Value=0.01作为阈值进行筛选，得到全部差异表达基因.阈值可以自行划定。
head(dif) #显示一部分报告结果

→差异表达基因分析补充资料

（二）注释

“注释实质上是一个ID映射过程，这里是将芯片探侦组的ID映射到基因国际标准名称（Gene symbol）和Entrez ID两种ID上。这里映射GI的目的是为了将GI映射到基因本体论（GO），然后做富集分析。”

library(annotate) #加载注释工具包
affydb <- annPkgName(CLLbatch@annotation, type="db") #获得基因芯片注释包名称
#BiocManager::install("affydb") 
library(affydb, character.only = TRUE) #　加载该包
dif$symbols <- getSYMBOL(rownames(dif), affydb) #根据每个探针组的ID获取对应基因Gene Symbol，并作为新的一列
dif$EntrezID <- getEG(rownames(dif), affydb) # 根据探针ID获取对应基因Entrez ID
head(dif)# 显示前几行

（三）统计分析及可视化

差异表达分析可以用富集分析方法做。这里介绍了GO的富集分析和KEGG通路的富集分析，分别由GOstats包和GeneAnswers包实现。
→关于富集分析（enrichment analysis）

1.GOterm检验
（1）结果表生成

library(GOstats)
entrezUniverse <- unique(unlist(mget(rownames(eset), hgu95av2ENTREZID))) # 提取HG_U95Av2芯片中所有探针组对应的EntrezID，unique函数去重
entrezSelected <- unique(dif[!is.na(dif$EntrezID), "EntrezID"])
# 提取差异表达基因及其对应的EntrezID

#设置GO富集分析的所有参数
params <- new("GOHyperGParams", geneIds=entrezSelected, universeGeneIds=entrezUniverse,
              annotation=affydb, ontology="BP", pvalueCutoff=0.001, conditional=FALSE,
              testDirection="over")
             
hgOver <- hyperGTest(params) #对所有的GOterm根据params参数做超几何检验

bp <- summary(hgOver) #生成所有GOterm的检验结果表

htmlReport(hgOver, file="ALL_go.html") # 同时生成所有GOterm的检验结果文件，每个GOterm都有指向官方网站的链接，以获得详细信息
head(bp) #显示局部结果

（2）基因表达热图

#BiocManager::install("pheatmap")
library(pheatmap)
selected <- eset[rownames(dif), ] #从基因表达矩阵中，选取差异表达基因对应的数据
rownames(selected) <- dif$symbols # 将selected矩阵每行的名称由探针组ID转换为对应的基因symbol

pheatmap(selected[1:20,], color = colorRampPalette(c("green", "black", "red"))(100),fontsize_row = 4, scale = "row", border_color = NA) #这里只画前20个基因的热图

（3）DAG关系图

library(Rgraphviz)

ghandle <- goDag(hgOver) #显著富集的GO term的DAG关系图
subGHandle <- subGraph(snodes=as.character(summary(hgOver)[,1]), graph = ghandle)
plot(subGHandle) #选一部分数据构建局部图

2.KEGG通路分析
（1）生成结果

#BiocManager::install("GeneAnswers")
library(GeneAnswers)  #加载所需包

# 选取dif中的三列信息构成新的矩阵，新一列必须是EntrezID
humanGeneInput <- dif[, c("EntrezID", "logFC", "P.Value")]

# 获得humanGeneInput中基因的表达值
humanExpr <- eset[match(rownames(dif), rownames(eset)), ]
humanExpr <- cbind(humanGeneInput[,"EntrezID"], humanExpr)

# 去除NA数据
humanGeneInput <- humanGeneInput[!is.na(humanGeneInput[, 1]),]
humanExpr <- humanExpr[!is.na(humanExpr[,1]),]

# KEGG通路的超几何检验
y <- geneAnswersBuilder(humanGeneInput, "org.Hs.eg.db", categoryType = "KEGG",
testType = "hyperG", pvalueT=0.1, geneExpressionProfile=humanExpr,verbose = FALSE)

getEnrichmentInfo(y)[1:6] #查看部分结果

（2）KEGG关系图

BiocManager::install("KEGG.db") #下载数据
yy<- geneAnswersReadable(y)
geneAnswersConceptNet(yy, colorValueColumn = "logFC", centroidSize = "pvalue", output="interactive")

最后画图一行出现了tcl/tk library not available的问题，没解决。
https://stackoverflow.com/questions/41133847/r-tcltk-not-loading-on-windows

（3）KEGG热图

yyy<- geneAnswersSort(yy,sortBy = "pvalue")
geneAnswersHeatmap(yyy)

引用资料：
https://www.jianshu.com/p/50bcf4ba9d8a
https://www.jianshu.com/p/fb4217512ec0
https://blog.csdn.net/tommyhechina/article/details/80356110
实验课课件