在进行目的芯片的差异分析之前，有两大拦路虎，一只是混杂因素；另一只，就是今天的问题，目的芯片的Agilent平台的芯片。因为是完全的新手，没有人带入门，只能抓虾。
网上的例子，包括金典的教材都是以Affy平台作为例子。
这里在GitHub上找了比人做过的例子，先运行一遍，弄懂别人代码的含义，再进行自己的芯片分析，我能怎么说，这个例子很坑爹。我直接运行自己的例子GSE84846

我们看一下RStudio的界面：

在红线出，Console右边有Terminal，这里可以进行一些有意思的操作，我在此处创建了一个文件夹，可以调给R使用，并且进入了Python：

设置路径，载入包

##### GSE84846   OSCC  加载必要的软件包和设置路径#####
source("https://bioconductor.org/biocLite.R")
# biocLite()
biocLite("sva")
getwd()
setwd("D:/RWork/Rworking")
library(GEOquery)
library(Biobase)
library(limma)

#GSE name of the study
GSEName = "GSE84846"

样本分组和筛选

# 字符串LNM来自GEO2R，是样本的分组信息，也用来筛除不需要的样本，它们对应“X"
#颈部淋巴结转移情况
LNM <- paste0("1011000000200X221010221100X002XX21X0200XX100202001",
               "00X21001X21XXX00021X02102022022221XXX2202X2202X00")
               
lnm <- c()

#对样本进行分组，有两个分组变量，分别是淋巴结转移情况LNM，和临床分期ClinStage，临床分期我去掉了。
for (i in 1:nchar(LNM)) { 
  lnm[i] <- substr(LNM, i, i)
 # clinstage[i] <- substr(ClinStage, i, i)
  }
# 删除了LNM标为 "X"的分组信息
sel <- which(lnm != "X")   #sel 全称selected
lnm <- lnm[sel]
#clinstage <- clinstage[sel]

获取原始数据

#### 用getGEOSUppFiles()函数获取原始数据#####
# 以下代码如果是再次运行则不需要运行
#rawdata = getGEOSuppFiles(GSEName)
getwd()

GSEpath = paste("./", GSEName, sep = "")
setwd(GSEpath)
#dir.create("CEL")
CELpath = "./CEL"
#untar(list.files()[2], exdir = CELpath)
getwd()

#下载平台注释文件
#GPL6480 <- getGEO("GPL6480", destdir = ".")

GPL6480 <- getGEO(filename = "GPL6480.soft")
colnames(Table(GPL6480))
head(Table(GPL6480))
write.csv(Table(GPL6480)[,c(1, 6, 7)], "GPL6480.csv")

读取数据

#获取原始数据，以下代码可行，并且删除了不需要的样本的txt格式文件
dat = list.files(CELpath, pattern = "txt")
dat <- dat[sel]#删除了LNM中标为 "X"的样本
dat
dat=read.maimages(files= dat, path = "./CEL", source="agilent")

前期处理

####背景校正和标准化处理####
raw_BGcorrected = backgroundCorrect(dat, "normexp", normexp.method = "rma", offset=50) 
raw_WithinArrays = normalizeWithinArrays(raw_BGcorrected, method = "robustspline")
raw_BGandNormalized = normalizeBetweenArrays(raw_WithinArrays, method="quantile")
#标准化之后就获得了表达矩阵
raw_aver = avereps(raw_BGandNormalized,ID=raw_BGandNormalized$genes$ProbeName)
class(raw_aver)

查看变化

背景校正后	芯片内校正后
芯片间校正后	均一化之后

创建设计矩阵

#### 创建设计矩阵 ####
lnm = sub("^1", "Node1", lnm)
lnm = sub("^0", "Node0", lnm)
lnm = sub("^2", "Node2", lnm)

lnm = as.factor(lnm)
#clinstage = as.factor(clinstage)
for (i in 1:length(raw_aver$targets$FileName)){
  raw_aver$targets$FileName[i] = substr(raw_aver$targets$FileName[i], 1, 10)
 
}
rownames(raw_aver$targets) = raw_aver$targets$FileName


pdata = as.data.frame(raw_aver$targets)
pdata$NodeStaging = lnm
design = model.matrix(~0 + NodeStaging, data = pdata)
colnames(design)[1] = "Intercept"
design
contr = makeContrasts(Node1vs0 = NodeStagingNode1 - NodeStagingNode0, 
                      Node2vs0 = NodeStagingNode2 - NodeStagingNode0,
                      Node2vs1 =NodeStagingNode2 - NodeStagingNode1, 
                      levels = design)

准备好的三个文件

#### 准备好的表达矩阵、分组矩阵和平台 ####
eset<- raw_aver #表达矩阵
summary(eset)
contr #实验设计矩阵
GPL = Table(GPL6480)[,c(1, 6, 7)]#平台注释文件
GPL

差异分析

差异分析

fit = lmFit(eset, design)
fit = eBayes(fit)
Node1 = topTable(fit, number=10000, coef = 1)
Node2 = topTable(fit, number=10000, coef = 2)
Node2vs1 = topTable(fit, number=10000, coef = 3)
head(Node1)
head(Node2)
Node2
colnames(Node1)
Node1[, c(7, 8, 11, 12, 13, 14, 15, 16)]
Node2[, c(7, 8, 11, 12, 13, 14, 15, 16)]
Node3[, c(7, 8, 11, 12, 13, 14, 15, 16)]
colnames(Node2vs1)

保存结果

Step1DownGenes = subset(Node1, Node1$logFC < -2 & Node1$adj.P.Val < 0.01)

Step2UpGenes = subset(Node2, Node2$logFC > 2 & Node2$adj.P.Val < 0.01)
Step2DownGenes = subset(Node2, Node2$logFC < -2 & Node2$adj.P.Val < 0.01)

Step2vs1UpGenes = subset(Node1, Node1$logFC > 2 & Node1$adj.P.Val < 0.01)
Step2vs1DownGenes = subset(Node1, Node1$logFC > 2 & Node1$adj.P.Val < 0.01)
getwd()
# Step1DownGenes
# 
write.table(Step1UpGenes[, c(7, 8, 11, 12, 13, 14, 15, 16)],
            file = "Node1UpGenes.txt", sep = "\t")
write.table(Step2UpGenes[, c(7, 8, 11, 12, 13, 14, 15, 16)],
            file = "Node2UpGenes.txt", sep = "\t")
write.table(Step2vs1UpGenes[, c(7, 8, 11, 12, 13, 14, 15, 16)],
            file = "Node2vs1UpGenes.txt", sep = "\t")

write.table(Step1DownGenes[, c(7, 8, 11, 12, 13, 14, 15, 16)],
            file = "Node1DownGenes.txt", sep = "\t")
write.table(Step2DownGenes[, c(7, 8, 11, 12, 13, 14, 15, 16)],
            file = "Node2DownGenes.txt", sep = "\t")
write.table(Step2vs1DownGenes[, c(7, 8, 11, 12, 13, 14, 15, 16)],
            file = "Node2vs1DownGenes.txt", sep = "\t")

绘制一个火山图