Bioconductor学习笔记
http://leweibo.blog.163.com/blog/static/179800089201271110403343/
rm(list = ls())
setwd("G:/EverBox/Cel/IgAN")
library(ALL)
data("ALL")
##----------安装相关包-----
source("http://www.bioconductor.org/biocLite.R")
biocLite("hgu95av2.db")
##-----------ExpressionSet及AffyBatch数据的基本操作----------------------------
#查看样本编码
sampleNames(ALL)
#查看ALL的表型变量
varMetadata(ALL)
#查看某表型变量
ALL$sex[1:10]
#查看芯片基因编码
featureNames(ALL)[1:100]
#提取表达谱数据为matrix结构
mat = exprs(ALL)
#提取出表型数据，mat数据对应基因的表达值列表
adf = phenoData(ALL)
#将数据分成亚组，并提取出来
vv = ALL[, 1:3]
#提取仅包含男性患者的数据
males = ALL[, ALL$sex == "M"]
###----------------A?ymetrix表达芯片数据处理----------------------------------
rm(list = ls())
setwd("G:/EverBox/Cel/IgAN")
setwd("G:/GSE11787")
##-----读入数据------
library("affy")
#使用read.AnnotatedDataFrame读入样本的表型变量数据
pd <- read.AnnotatedDataFrame("pd.txt",header=TRUE)
pData(pd)
#读入CEL结构的芯片数据，返回AffyBatch类数据
#ReadAffy()如果不带任何参数，将默认读入所有工作目录下的CEL文件
#可以指定读入的CEL文件：myAB = ReadAffy(filenames=c("a1.cel", "a2.cel", "a3.cel"))
Data <- ReadAffy(filenames=rownames(pData(pd)),phenoData=pd,verbose = TRUE)
#查看样本名称
sampleNames(Data)
PM = pm(Data)
MM = mm(Data)
gnames <- geneNames(Data)
##-----质量评价与控制------
#尺度因子（scaling factor）=100/芯片上所有基因定量的平均值
#Affymetrix公司标准：用于比较的芯片之间的尺度因子的比例必须小于3
#P/M/A分类：P(present)就是PM和MM的值有显著差别，指代该基因被检测到；
#?????????? A(absent)就是PM和MM的值没有显著差别，指代该基因未被检测到；
#?????????? M(marginal present)就是介于这两者之间的临界状态。
#如果大部分的基因都未被检测到，说明实验出现了问题。
#?? 而在多组平行实验中，如果其中一组的被检测到的基因和其它组有显著的差别，那说明该实验可能出现了问题。
#对于每一块芯片，所有的MM值做出统计可以得到背景噪音的平均值，最小值和最大值。
#?? 如果背景噪音的相对于其它平行组来说平均值过大，那也说明该实验可能出现了问题。
#?? 而往往高的噪音都伴随着低的被检测到的基因比例(low percentage present)
#?? 所以这两个可以做为判断实验是否合格的一个指标。
#
#Affymetrix内参:(1)β-actin和GAPDH:RNA降解程度的内参;
#?????????????? (2) 嵌入探针(spike-in probe):BioB(浓度下限)，BioC，BioD 和CreX(浓度上限)。它们的RNA是在标记的过程中加入样品体系的。如果BioB不能被MAS5算法标记为P，说明该芯片的敏感度没有达标。这很有可能是芯片本身的问题。
#如果有一个芯片各项指标都不太正常，尤其是BioB无法检测到，可以判定为芯片故障。
library(affy)
library(simpleaffy)
#library(AnnotationDbi)
Data.qc = qc(Data)
#显示质量分析结果QCStats
Data.qc
#以下为逐一显示
#质量探针actin、gapdh 5’3’或者5’M比值。5’3’比值过大，也说明实验存在着质量问题。
ratios(Data.qc)
#质量探针
qcProbes(Data.qc)
#嵌入探针
spikeInProbes(Data.qc)
#尺度因子
sfs(Data.qc)
#--------质量控制的结果绘图展示：简单地讲，如果标记为蓝色，说明正常，如果标记为红色，说明可能存在质量问题
#浅蓝色的竖条代表着尺度因子正常的取值范围，它会依照实验具体数据来计算出这个范围，通常它应该是在三倍以内
#横轴所标记的数字就是尺度因子的座标
#如果所有的一组需要相互比较的芯片间的尺度因子都落在了蓝色范围内，它会以蓝色线条及蓝色端点显示，表明这些芯片可以相互比较，如果标记为红色（比如说这个示例），那就意味着它们不能相互比较。
#??? 本组实验C3超出了此范围，因此尺度因子被标记为红色
#最左侧是样品的名字，而后是两个数字，上面的以百分比形式出现的是P所占比重，下面的数字表明平均背景噪音。如果它们标记为红色，说明存在质量问题。
#而P所占的比重应该在平行实验间较为一致。而过低的P比重（<20%）说明制样过程可能存在问题
#如果图中出现红色的BioB字样，说明该样品嵌入探针未能检测到BioB
#actin和GAPDH 3’/5’比值 也分别以△和○表示出来。对于actin的3’/5’应该落在3以内，而对于GAPDH应该落在1左右。如果超过了设定的标准，就会以红色显示
#绘图
plot(Data.qc)
#--------使用FitPLM生成权重，残差及NUSE图像
#导入affyPLM库。这个库会基于RMA算法来对数据进行预处理，然后给出与质量控制相关的直观图像
#权重图中，绿色代表较低的权重（接近0），白色、灰色代表较高的权重（接近1）
#残差图中，红色代表正的高残差，白色代表低残差，蓝色代表负残差
#sign.resids图使用红色表达正的残差，蓝色表达负的残差
#权重和残差都是随机分布的，所以应该看到的是绿色均匀分布的权重图，或者红蓝均匀分布的残差图
library(affyPLM)
Pset <- fitPLM(Data)
x11()
par(mfrow=c(2,2))
image(Data[,2]);image(Pset,type="weights",which=2,main="Weights");
image(Pset, type="resids", which=2, main="Residuals");
image(Pset, type="sign.resids", which=2, main="Residuals sign");
#相对对数表达（Relative Log Expression(RLE)）箱线图
#如果有一个芯片的表现和其它的平行组都很不同，那说明它可能出现了质量问题
library(RColorBrewer)
colors <- brewer.pal(12,"Set3")
x11()
Mbox(Pset, ylim = c(-1,1), col=colors,main="RLE")
#相对标准差箱线图（Normalized Unscaled Standard Errors(NUSE)）
#NUSE是一种比RLE更为敏感 的质量检测手段。如果你在RLE图当中对某组芯片的质量表示怀疑，那当你使用NUSE图时，这种怀疑就很容易被确定下来
#NUSE是某芯片基因标准差相对于全组标准差的比值
#我们期待全组芯片都是质量可靠的话，那么，它们的标准差会十分接近，于是它们的NUSE值就会都在1左右
#然而，如果有实验芯片质量有问题的话，它就会严重的偏离1，进而影响其它芯片的NUSE值偏向相反的方向。当然，还有一种非常极端的情况，那就是大部分芯片都有质量问题，但是它们的标准差却比较接近，反而会显得没有质量问题的芯片的NUSE值会明显偏离1，所以我们必须结合RLE及NUSE两个图来得出正确的结论。
boxplot(Pset, ylim= c(0.95, 1.20), col=colors, main="NUSE")
#信号强度分布曲线（直方图）
hist(Data[,1:6],col=colors)
legend("topright",rownames(pData(Data))[1:6],col=colors,lwd=1,inset=.05)
#信号强度分布箱线图
boxplot(Data,col=colors)
#绘制R-I(ratio-intensity)图，又称MA图
#X轴为A值[log2(Exp)+log2/(Control)]/2，代表点的整体信号强度
#Y轴为M值log2(Exp)－log2/(Control)表示点的两通道信号差
#可由MA图观察芯片数据是否存在强度依赖的系统偏移以及信号差异点的比例
#Media指示出中值偏离0的程度，IQR是四分位距（interquartile range），两者都是越接近0越好
MAplot(Data[,1:6],pairs=TRUE,plot.method="smoothScatter")
#标准化后的MVA图
Data.rma <- rma(Data)
MAplot(Data.rma[,1:6],pairs=TRUE,plot.method="smoothScatter")
#RNA降解图
#RNA通常都是从5’端开始降解的，所以我们可以想象5’端的探针荧光强度应该低于3’端的荧光强度
#斜率接近0说明降解较少或者全部被降解
#实际实验中，降解较少是不太可能的，所以当斜率过小时，很可能表明RNA降解严重
Data.deg <- AffyRNAdeg(Data[1:6])
plotAffyRNAdeg(Data.deg)
plotAffyRNAdeg(Data.deg,col=colors)
legend("topleft",rownames(pData(Data)),col=colors,lwd=1,inset=.05)
#PCA分析 PCA(principal component analysis)，主成份分析
#经过PCA分析之后，平行实验所提供的基因芯片数据应该聚扰在一起，而不同设计的实验所提供的基因芯片数据应该分离
#这可以帮助我们很快的识别出一组平行实验当中，有哪些数据是可靠的，而哪些数据可以被放弃。
#实际的操作可能更好理解。我们使用gcrma(rma的一种扩展)来对数据进行预处理，然后使用affycoretools库当中的plotPCA来进行PCA分析作图。我们注意到它的结果和NUSE的结果保持一致。
library(affycoretools)
eset=gcrma(Data)
eset1=rma(Data)
plotPCA(eset,groups =? as.numeric(as.factor(pData(eset)[,1])),
??????? groupnames =levels(as.factor(pData(eset)[,1])))
eset1=eset[,7:21]
plotPCA(eset1,groups =? as.numeric(as.factor(pData(eset1)[,4])),
??????? groupnames =levels(as.factor(pData(eset1)[,4])))
#QCReport能够生成6页pdf的质量报告
library(affyQCReport)
QCReport(Data,file="ExampleQC.pdf")
#---另外一种使用FitPLM生成权重，残差及NUSE图像的方法
library("affyPLM")
dataPLM = fitPLM(Data)
boxplot(dataPLM, main="NUSE", ylim = c(0.95, 1.22), outline = FALSE, col="lightblue", las=3,
??????? whisklty=0, staplelty=0)
Mbox(dataPLM, main="RLE", ylim = c(-0.4, 0.4),outline = FALSE, col="mistyrose", las=3,
???? whisklty=0, staplelty=0)
#----------------------------------------------------
#样本聚类
#原始数据读入，经AffyBatch目标转成ExpressionSet目标后，为提高后续分析（如差异表达基因的检测）的统计功效，往往需要进一步经过Detection Call Filter和IQR filter等过滤
#常规做法是先筛选出差异表达基因，然后只用差异表达基因进行聚类分析，本示例直接用了过滤后的数据集，聚类图的效果稍差一点
library(genefilter)
eset=rma(Data)
#不同的R语言包可能会出现相同的函数名称,dist2函数出现在两个包中，为避免冲突函数前面放上包的名称作为前缀
dd = genefilter::dist2(log2(exprs(eset)))
diag(dd) = 0
dd.row <- as.dendrogram(hclust(as.dist(dd)))
row.ord <- order.dendrogram(dd.row)
library("latticeExtra")
legend = list(top=list(fun=dendrogramGrob,args=list(x=dd.row, side="top")))
lp = levelplot(dd[row.ord, row.ord],scales=list(x=list(rot=90)), xlab="",
?????????????? ylab="", legend=legend)
lp
#很显然C3是显著异常样本，我们将C3剔除
badArray = match("C3.CEL", sampleNames(CLLbatch))
Data = Data[, -badArray]
##--------------数据分析---------
#以CLL数据为例
library("CLL")
data("CLLbatch")
data("disease")
rownames(disease) = disease$SampleID
#将CLLbatch的sampleNames去掉CEL后缀
sampleNames(CLLbatch) = sub("\\.CEL$", "",sampleNames(CLLbatch))
#将表型变量加入CLLbatch数据结构
mt = match(rownames(disease), sampleNames(CLLbatch))
vmd = data.frame(labelDescription = c("Sample ID",
????????????????????????????????????? "Disease status: progressive or stable disease"))
phenoData(CLLbatch) = new("AnnotatedDataFrame",
????????????????????????? data = disease[mt, ], varMetadata = vmd)
#其中一个Disease变量为NA，将去去除
CLLbatch = CLLbatch[, !is.na(CLLbatch$Disease)]
#经质量分析后发现CLL1为异常数据，将其剔除
badArray = match("CLL1", sampleNames(CLLbatch))
CLLB = CLLbatch[, -badArray]
#--数据预处理：芯片标准化（归一化），将A?yBatch数据转化为ExpressionSet结构
CLLrma = rma(CLLB)
e= exprs(CLLrma)
#--Ranking and ?ltering probe sets
#数据过滤nsFilter:非特异性的筛选方法(不依赖于样本的表型变量)
#详见：genefilter包中的nsFilter方法? 使用genefilter包对芯片数据进行筛选
#参见nsFilter {genefilter}
#remove.dupEntrez=FALSE
library("genefilter")
library("hgu95av2.db")
CLLf = nsFilter(CLLrma, remove.dupEntrez=FALSE,var.cutof =0.5)$eset
#Summary statistics and tests for ranking
#用rowttests检验（t检验或F检验）比较组建表达差异，dm指基因在各组见的表达差异值
CLLtt = rowttests(CLLf, "Disease")
names(CLLtt)
head(CLLtt)
#计算行平均值（所有样本中基因表达的平均值）
a = rowMeans(exprs(CLLf))
#绘图展示结果
plot(a,CLLtt$dm)
#rank(a)，以a的大小排序。用于分析dm是否与a有关
plot(rank(a),CLLtt$dm)
#当样本量较少时，用t检验比较组间基因表达差异并不妥当
#以下为改良的t检验： modi?ed t-statistic based on an empirical Bayes moderation approach
library("limma")
design = model.matrix(~CLLf$Disease)
CLLlim = lmFit(CLLf, design)
CLLeb = eBayes(CLLlim)
#注意：当组间样本大于10时，使用贝叶斯改良t检验的优势就不明显了（当然也没有什么缺点）。
#volcano plot:横轴表示dm，纵轴表示log10(P).横线上方的点表示P<0.01
lod = -log10(CLLtt$p.value)
plot(CLLtt$dm, lod, pch=".", xlab="log-ratio",
???? ylab=expression(-log[10]~p))
abline(h=2)
#points(CLLtt$dm[CLLtt$p.value<0.001], lod[lod>3], pch =18, col = "red")
#多重检验问题：limma包中的toptable函数提供了多种方法进行多重检验校正
#coef=2是指我们关注回归曲线的第二项回归系数，即我们感兴趣的变量系数
#coef=1是第一项回归系数，为常数项。n表示返回的基因数
#BH表示筛选方法
tab = topTable(CLLeb, coef=2, adjust.method="BH", n=10)
genenames = as.character(tab$ID)
#Annotation注释
#载入注释包
library("annotate")
#查看所需要的注释包
annotation(CLLf)
#载入所需要的注释包
library("hgu95av2.db")
#将芯片编码号与注释信息对应，并查看注释信息（获得Entiz基因ID、基因代号等）
ll = getEG(genenames, "hgu95av2")
sym = getSYMBOL(genenames, "hgu95av2")
#创建HTML页面查看
tab = data.frame(sym, signif(tab[,-1], 3))
htmlpage(list(ll), othernames=tab,
???????? filename="GeneList1.html",
???????? title="HTML report", table.center=TRUE,
???????? table.head=c("Entrez ID",colnames(tab)))
browseURL("GeneList1.html")
#使用anna?y包创建更详细的注释信息
library("KEGG.db")
library("hgu95av2.db")
library("annaffy")
atab = aafTableAnn(genenames, "hgu95av2.db", aaf.handler())
saveHTML(atab, file="GeneList2.html")
#--------------表达差异分析----------------
library("Biobase")
library("genefilter")
library("ALL")
data("ALL")
#选择亚组
ALL$mol.biol
ALL_bcrneg$mol.biol
bcell = grep("^B", as.character(ALL$BT))
moltyp = which(as.character(ALL$mol.biol) %in% c("NEG", "BCR/ABL"))
ALL_bcrneg = ALL[, intersect(bcell, moltyp)]
ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol)
#非特异性筛选
#有部分基因表达量在各个芯片上非常接近，我们需要将其剔除
library("genefilter")
#rowSds计算出每行（每个基因）表达量的标准差
sds = rowSds(exprs(ALL_bcrneg))
# shorth() calculates the midpoint of the shorth(the shortest interval containing half of the data)
#在大多数情况下用shorth()可以确定筛选临界值
sh = shorth(sds)
hist(sds, breaks=50, col="mistyrose", xlab="standard deviation")
abline(v=sh, col="blue", lwd=3, lty=2)
#剔除表达量的标准差在sh以下的基因
ALLsfilt = ALL_bcrneg[sds>=sh, ]
dim(exprs(ALLsfilt))
#筛选策略2：用nsFilter筛选，见上文
#筛选策略3：剔除所有芯片上均为低表达或不表达的基因
#差异表达分析(t检验)
table(ALLsfilt$mol.biol)
tt = rowttests(ALLsfilt, "mol.biol")
names(tt)
hist(tt$p.value, breaks=50, col="mistyrose", xlab="p-value",main="Retained")
#将非特异筛选掉的部分
ALLsrest = ALL_bcrneg[sds<sh, ]
ttrest = rowttests(ALLsrest, "mol.biol")
hist(ttrest$p.value, breaks=50, col="lightblue",xlab="p-value", main="Removed")
#hist图表明有明显差异基因大部分没有因非特异筛选而丢失
#多重检验的校正
#mt.rawp2adjp函数用于校正t检验的p值，BH表示校正方法
library("multtest")
mt = mt.rawp2adjp(tt$p.value, proc="AB")
g = featureNames(ALLsfilt)[mt$index[1:10]]
library("hgu95av2.db")
links(hgu95av2SYMBOL[g])
mb = ALLsfilt$mol.biol
y = exprs(ALLsfilt)[g[2],]
ord = order(mb)
plot(y[ord], pch=c(1,16)[mb[ord]],col=c("black", "red")[mb[ord]],
???? main=g[2], ylab=expression(log[2]~intensity),xlab="samples")
#实例
library("ALL")
data("ALL")
#获得B细胞白血病的亚组指针
bcell = grep("^B", as.character(ALL$BT))
#获得NEG或者BCR/ABL型白血病的亚组指针
moltyp = which(as.character(ALL$mol.biol) %in% c("NEG", "BCR/ABL"))
#获得同时符合上述两个条件的亚组
ALL_bcrneg = ALL[,intersect(bcell,moltyp)]
#清除空白的levels
ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol)
#查看标准差大小是否与基因表达水平相关
#如果标准差与表达水平相关，使用非特异性筛选将标准差小基因的删除可能会导致选择偏倚
#我们将基因表达水平（x轴）与基因表达水平的标准差（y轴）绘图
library("vsn")
meanSdPlot(ALL_bcrneg)
#可见基因表达水平（x轴）与基因表达水平的标准差（y轴）仅存在非常弱的相关性；
#因此,下一步我们可以将标准差小的基因剔除
#本例使用筛选的sds较大：设置为0.8，通常不需要设置这么大
sds = esApply(ALL, 1, sd)
sel = (sds > quantile(sds, 0.8))
ALLset1 = ALL_bcrneg[sel,]
#差异表达计算
library("genefilter")
tt = rowttests(ALLset1, "mol.biol")
plot(tt$dm, -log10(tt$p.value), pch=".",xlab = expression(mean~log[2]~fold~change),
???? ylab = expression(-log[10](p)))
#多重比较
library("multtest")
cl = as.numeric(ALLset1$mol.biol=="BCR/ABL")
resT = mt.maxT(exprs(ALLset1), classlabel=cl, B=1000)
ord = order(resT$index) ## the original gene order
rawp = resT$rawp[ord] ## permutation p-values
hist(rawp, breaks=50, col="#B2DF8A")
sum(resT$adjp<0.05)
res = mt.rawp2adjp(rawp, proc = "BH")
sum(res$adjp[,"BH"]<0.05)
?
