这章开始捡的是芯片数据分析
芯片简介
基因芯片就是一系列微小特征序列的（通常是DNA探针，也可能是蛋白质）的集合，它们可以被用于定性或者定量检查样品内特异分子的成份
可以对整个基因组上的某些基因marker进行定量分析
利用Bioconductor分析Affymetrix microarray数据
读取本地数据
library(BiocManager)
## 安装包
install("arrays")
install(c("affy", "limma"))
#载入包
library(affy) ##加载库文件
#读取数据
Data <- ReadAffy("ArabidopsisATH1-121502.CEL")
library(oligo)
xysFiles <- list.celfiles('myXYSs', full.names=TRUE)
rawData <- read.xysfiles(xysFiles)
读取在线数据
library(GEOquery)
gset <- getGEO("GSE46106", GSEMatrix =TRUE)
length(gset)
gset <- gset[[1]]
head(pData(gset)[,1:5])
QC
任何测序数据第一步都是QC，检测质量是否达标
library(affydata)
data(Dilution)
deg <- AffyRNAdeg(Dilution)
#生成总结
summaryAffyRNAdeg(deg)
#可视化
plotAffyRNAdeg(deg)
总结
斜率太大代表大部分降解了，斜率小代表降解小
可视化
标准化
采用gcrma包的函数：gcrma()
library(gcrma)
eset <- gcrma(Dilution)
我们看看标准化前后的表达量的差别：
boxplot(Dilution)
前
boxplot(eset)
后
MA图：
y <- (exprs(Dilution)[,c("20B","10A")]) ##载入Dilution中的20A,20B至y
ma.plot(rowMeans(log2(y)),log2(y[,1])-log2(y[,2]),cex=1, plot.method = "smoothScatter")
title("Pre-Norm Dilutions Dataset(array 20B v 20A)")
前
x <- exprs(eset)[, c("20B", "10A")]
ma.plot(rowMeans(log2(x)),log2(x[,1])-log2(x[,2]),cex=1, plot.method = "smoothScatter")
title("Post-Norm Dilutions Dataset(array 20B v 20A)")
后
若MA图，M的值偏于0，那么在差异表达分析获得的差异基因是比较准确的，否则，找到的差异基因用很大的假阳性
PCA
PCA是看聚类情况的，也就是我们处理的样品与对照组是否能够区分的开，生物学重复是否能够聚到一起
library(affycoretools)
plotPCA(eset)
差异表达
我们用limma来进行差异表达分析
首先设计好设计矩阵
实验设计
library(affydata)
data(Dilution)
library(limma)
library(gcrma)
eset <- gcrma(Dilution)
#设计矩阵
f <- factor(c("liver20", "liver20", "liver10", "liver10"))
design <- model.matrix(~ 0 + f)
colnames(design) <- levels(f)
contrast <- c("liver20-liver10")
cont.matrix <- makeContrasts(contrasts = contrast, levels=design)
df <- data.frame(Sample=c("A", "B", "C", "D", "E", "F"),
condition=c("WT", 'WT', "WT", "KD", "KD", "KD"),
pair = c(1, 2, 3, 1, 2, 3),
batch = c(1, 2, 1, 2, 1, 2))
des <- model.matrix(~0 + condition + pair + batch, data=df)
#差异表达
fit <- lmFit(eset, design)
fit1 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit1)
topTable(fit2, coef = contrast, n=5)
如果设计条件更多点：
df <- data.frame(Sample=c("A", "B", "C", "D", "E", "F", "G", "H"),
condition=c("WT", 'WT', "WT", "WT", "KD", "KD", "KD", "KD"),
time = c(0, 0, 3, 3, 0, 0, 3, 3), stringsAsFactors = FALSE)
df$f <- factor(paste(df$condition, df$time, sep="."))
des <- model.matrix(~0 + f, data=df)
colnames(des) <- levels(df$f)
#差异表达
fit <- lmFit(eset, design)
fit1 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit1)
topTable(fit2, coef = contrast, n=5)
注释
对于ExpressionSet对象，如果对象中含有featureData的话，当使用topTable等输出结果时，就会自动生成注释信息。所以，基因芯片的注释，多半可以在构建ExpressionSet对象时就准备好
我们上代码：
library(affydata)
data("Dilution")
library(annotate)
affydb<-annPkgName(Dilution@annotation,type="db")
require(affydb, character.only=TRUE, quietly=TRUE)
若没有featureData，则
library(annaffy)
library(limma)
library(gcrma)
eset <- gcrma(Dilution)
#检查是否
featureData
head(fData(eset))
#若没有，生成featureData
fD <- data.frame(symbol=as.character(aafSymbol(rownames(eset), affydb)),
gene =as.character(aafUniGene(rownames(eset), affydb)),
row.names=rownames(eset))
metadata <- data.frame(labelDescription=c("gene symbol", "UniGene"), row.names=colnames(fD))
adf <- new("AnnotatedDataFrame", data=fD, varMetadata=metadata)
featureData(eset) <- adf
#差异分析
fit <- lmFit(eset, design)
fit1 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit1)
topTable(fit2, coef = contrast, n=5)
