安捷伦芯片原始数据处理
双通道芯片有时候实验设计挺复杂的，agilent的原始数据数据处理在中文互联网上也不算常见。
本文借助limma的帮助文档，完成一篇使用agilent双色表达分析，研究肺鳞癌早期肿瘤发生和免疫逃避机制的nature的原始数据处理和文章复现。
由于本文涉及了比较详细的双通道芯片上游分析及原理，因此预计将这篇文章分成上下两部分推送。
背景知识
双通道芯片RNA样本制备，与芯片杂交：
eg：GSE33479
零、 数据下载，读取
型号：Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (Probe Name version)
文章
可以看到有10*32个EA对照，后面做MAplot的时候可以高亮出来
基因表达谱的获取和分析
通俗讲就是：「Cy3」标记了「从未吸烟」的「正常人」的支气管组织
「Cy5」分别标记了SCC，原位癌，组织学异常低荧光，组织学正常荧光正常的122个「吸烟或曾经吸烟」样本。
GSE_number<-"GSE33479"
library(limma)
library(R.utils)
#批量解压数据，只能运行一次，再运行时候会报错（压缩包不能解压两遍）
files<-list.files('./RAWdata',pattern = "txt",full.names = T,recursive = F)
for(i in 1:length(files)){
gunzip(files[i], remove = T)
}
#列出解压完成后的./RAWdata文件夹下面有txt的文件，
txt<-list.files('./RAWdata',pattern = "txt",full.names = T,recursive = F)
#对于Agilent的gpr文件或者txt文件，可以使用limma的read.maimages来读取。
#文件source参数需根据图像分析软件进行选择：
#"generic", "agilent", "agilent.median", "agilent.mean", "arrayvision", "arrayvision.ARM",
#"arrayvision.MTM", "bluefuse", "genepix", "genepix.custom", "genepix.median", "imagene",
#"imagene9", "quantarray", "scanarrayexpress", "smd.old", "smd", "spot" or "spot.close.open"。
#对于单色芯片，注意将green.only设置成TRUE.
RG <- read.maimages(txt, source="agilent",green.only = F)
原始数据读进来是一个「RGList」对象，由limma包产生。
帮助文档写的很清楚，我翻译出来给大家看一下吧。
一、 RGList对象结构
Red, Green Intensity List - class
什么意思：我的理解是红色绿色荧光信号强度列表，包含双通道芯片原始的荧光信号强度的矩阵和探针/基因对应到哪个荧光点的注释信息。
Description
一个基于列表的S4类，用于存储一批荧光斑点微阵列的红色和绿色通道前景和背景强度。RGList对象通常由read.maimages函数创建。
Slots/List Components
RGList对象可以用创建S4类对象的函数new()来创建。需要包含以下组分：
「R」 | 「包含红色荧光(cy5)前景信号强度的实数矩阵。行对应点，列对应芯片。」 |
「G」 | 「包含绿色荧光(cy3)前景信号强度的实数矩阵。行对应点，列对应芯片。」 |
「Optional components include 可选组分」
「Rb」 | 「包含红色荧光(cy5)背景信号强度的实数矩阵。」 |
「Gb」 | 「包含绿色荧光(cy3)背景信号强度的实数矩阵。」 |
「weights」 | 「含有相对点质量权重的与R（红色前景信号强度）相同维数的数值矩阵。元素应该是非负的。」 |
「other」 | 「列表中包含的其他矩阵，维度同R和G一致」 |
「genes」 | 「包含探针信息的数据框，每个荧光点必须要对应一行，可以有任意列」 |
「targets」 | 「含有RNA样本信息的数据框，行对应芯片数量，可能有任意列。通常包含Cy3和Cy5，指定将哪个RNA杂交到每个芯片。」 |
「printer」 | 「包含用于在芯片上打印点的过程信息的列表。详细请查看limma帮助文档的printlayout」 |
有效的RGList对象可以包含其他可选组份，但是所有探针或芯片的信息都应该包含在上述组件中。
「我们来看看我们这个数据集读进来的对象存储了些什么东东：」
1. targets
#1. targets
head(RG$targets)
dim(RG$targets)
targets存了文件名，一个文件对应到一个芯片的RNA样本嘛。
可以看到我们读进来了122个文件，做了122个样品，稍微记一下，待会儿要考。
2. genes
#2. genes
RG$genes->RG_gene
View(RG_gene)
dim(RG$genes)
可以看到有探针序列，探针id基因名称和描述等信息
动动短期记忆记一下这里dim的行的结果为45015行，也就是有45015个点，列是根据芯片平台包含信息多少决定的，所以换个平台就不一定是这么多行了。
3. G,Gb,R,Rb
#3. G,Gb,R,Rb
head(RG$G)
RG$G->G
paste0("c",1:ncol(G))->colnames(G)
dim(RG$G);dim(RG$Gb);dim(RG$R);dim(RG$Rb)
可以看到是列名为文件名，行为基因的矩阵，里面每一个值都是该样本的对应荧光点的信号值，以下是简化了行名之后看到的结果。
回旋镖这么快就扎回来了，我们在genes数据框下看到有45015个点及对应的注释信息，这里我们看到绿色前景荧光，绿色背景荧光，红色前景荧光，红色背景荧光都是有45015行，对应每个点的信号强度，而122列则是对应targets有122个RNA样品
4. source
来源就是平台，芯片厂商安捷伦了
使用MAplot检查RGList中每个样本质量
我们将对照基因整理成controlStatus函数需要的格式，用于在MAplot检查时高亮出来
通常情况下要设置对照探针需要使用STF(Spot Types Files)，它通常用于区分对照探针和对应基因的常规探针，以及区分阳性对照与阴性对照、校准对照的比率等。
STF应该有一个「SpotType列」（必须要），给出不同点类型的名称。一个或多个其他列应具有与genelist中的列相同的名称，并且应包含足以识别斑点类型的模式或正则表达式。
以ApoAI数据的STF为例，图片出自limmauserguide：
在本例中，列ID和列Name在genelist中，并包含要「匹配的模式」。星号是通配符，可以表示任何内容。小心使用大写或小写，不要插入任何多余的空格。剩下的列提供了要与不同类型的点关联的颜色。此代码假定探测注释data.frame包括列ID和列Name。如果GenePix已用于图像分析，则通常如此，但其他图像分析软件可能使用其他列名。
GEO中没有STF文件，我们自己建一个数据框来放点类型相关的信息。因为在RGList对象中GeneName列包含了点类型的信息我们提出来看一下
RG$genes$GeneName->genename
sort(table(genename),decreasing = T) %>% head(.,100)
#可以看到有DarkCorner,NegativeControl,EQC,10*32的E1A,GE_BrightCorner,除此之外应该都是对应基因的探针
genename[str_detect(genename,"E1A_")] %>% length
#length()结果表明确实是320个
#创建可以被controlStatus函数识别的数据框
#这里的SpotType列会作为MAplot的图例，GeneName列与RG$genes中为symbol的列同名，内容为要检索匹配到的名称（可以想象成str_detect函数的检索），然后color列就是MAplot中点的颜色了
spottypes<-data.frame(SpotType=c("gene","EQC","DarkCorner","NegativeControl","GE_BrightCorner","E1A"),
GeneName=c("*","EQC","DarkCorner","NegativeControl","GE_BrightCorner","E1A*"),
color=c("#696969","red","yellow","darkgreen","blue","purple"))
#我们先画一个测试一下
plotMA(RG[,1])
#新建no_normalized_MAplot文件夹，用于存放未标准化的MAplot
dir.create("./no_normalized_MAplot")
#plotMA3by2是生成每张图3*2个芯片MAplot的函数，就是plotMA的自动拼图版本...
plotMA3by2(RG,path = "./no_normalize_MAplot")
#可以看到读入原始数据和画MAplot是限速步骤，因为122个样本量太大了，把读好的数据存成.Rdata方便日后取用
save(RG,GSE_number,file=paste0(GSE_number,"_agilent_RGList.Rdata"))
这张图中黄色和蓝色的corner点应该是芯片物理空间上的角落，不对应到基因。绿色是阴性对照，紫色E1A则是产品目录中人为设置的有指定倍数差异表达的外源RNA参照，理想情况下2的M值次幂应该等于设计时的差异表达倍数。灰色的就是基因，这张图对应的还没有进行背景校正和芯片间标准化。
二、背景矫正、芯片内标准化、芯片间标准化
先进行背景校正和芯片内标准化
rm(list = ls())
GSE_number<-"GSE33479"
load(paste0(GSE_number,"_agilent_RGList.Rdata"))
RG<-limma::backgroundCorrect(RG, method="normexp", offset=25)
#normalizeWithinArrays使用的默认标准化方法是打印探头的loess标准化（printtiploess）
#有一些不合适的情况是值得注意的。例如，Agilent芯片没有打印头组，所以应该使用全局勒斯标准化：
MA<-normalizeWithinArrays(RG, method="loess")
class(MA)
#也先画一个
plotMA(MA[,1])
dir.create("./normalized_MAplot")
#绘制芯片内标准化后的MAplot到工作目录下的normalized_MAplot里
plotMA3by2(MA,path = "./normalized_MAplot")
#检查一下M值得箱线图看箱体范围相差大吗，如果相差大则表明需要normalizeBetweenArrays
pdf("Mboxplot.pdf")
boxplot(MA$M)
dev.off()
#boxplot的结果显示需要进行芯片间标准化
MA <- normalizeBetweenArrays(MA,method="scale")
save(MA,GSE_number,file=paste0(GSE_number,"_agilent_MAList.Rdata"))
芯片内标准化之后大多数基因的M值都在0附近了（毕竟差异表达基因是少数，稳定表达的占多数），这个图的E1A对照的回归看着不是很好，
可以在normalizewithinarrays函数中使用参数“method="control"设置对照点的loess回归，但是使用前提是你需要知道稳定表达的点，且对照点足够多，详见limmauserguide，推荐使用默认的标准化方法就行。
注意本文没有使用默认的printtioloess标准化，而是使用全局的loess标准化是因为agilent芯片制造原理上是没有tip的。
又多一个新的对象，MAList
别慌，听我慢慢巴扯
M-value, A-value Expression List - class
还是解释下标题，MAList是一个列表存储芯片M值（M值=前景信号强度-背景信号强度）和A值（A值=前景信号强度+背景信号强度）的列表。
Description
一个简单的基于列表的类，用于存储一批芯片荧光点的M值和A值。MAList对象通常在标准化过程中由normalizewithinArrays或MA.RG函数创建。
Slots/List Components
MAList对象可以用创建S4类对象的函数new()来创建。需要包含以下组分：
「M」: | 「一个包含M值的数值矩阵，行与荧光点对应，列与芯片（做了几个RNA样本）对应」 |
「A」: | 「一个包含A值的数值矩阵，行列对应关系和M相同」 |
Optional components include:
可选组分包括：
「weights」: | 「一个和M有着相同维度的数值矩阵，包含每个点相对的质量。里面的元素需要是非负的。」 |
「other」: | 「补充其他信息矩阵，和M具有相同维度」 |
「genes」: | 「包含探针信息的数据框，每个一行都对应一个荧光点，列数不限」 |
「targets」: | 「包含RNA样品信息的数据框，每行对应M的每一列。通常包括Cy3和Cy5指定哪个RNA杂交到哪个芯片上了」 |
「printer」: | 「包含打印芯片的探头的信息」 |
有效的MAList对象可以包含其他可选组分，但是所有的探针或芯片信息都应该包含在上述组分中。
三、整理临床信息，提取分组信息
我们前文有写过文献中双通道的荧光都标记了什么RNA：
?将cDNA在体外转录成互补cRNA，并使用Cy5染料标记来自122个感兴趣样本的RNA, Cy3染料标记参考RNA。参考RNA从16名从未吸烟的人的正常支气管活检中等量提取(Agilent Technologies)。在PCR扩增和荧光标记后，根据芯片制造商(Agilent Technologies)的建议，cRNA被杂交到Whole Human Genome 4 x 44K双色芯片上。使用genesspring GX version 7.3.1软件进行进一步的标准化步骤:(1)每个斑点(按对照通道划分)，(2)每个芯片(标准化到芯片间的中位数表达值)和(3)每个基因(标准化到患者间的中位数表达值)。基因表达测量的报告展示了每个探针和基因的相对丰度，也就是说，与从未吸烟的人的正常活检相比，研究样本中红色和绿色强度(Cy5/Cy3)之间的比率。 ?
因为双通道芯片的实验设计都挺复杂的。为方便大家理解这里放个表格：
Cy5 | Cy3 | |
|---|---|---|
first step:「normal bronchial tissue」 | 「吸烟/以前吸烟者」：组织学正常且正常荧光或低荧光（stage 0-1）组织学增生且正常荧光或低荧光（stage 2） | 「从未吸烟者」的正常样本 |
second step:「low-grade lesions」 | 「吸烟/以前吸烟者」：化生（stage 3）轻度或中度异生（stage 4-5） | 「从未吸烟者」的正常样本 |
third step:「high-grade lesions」 | 「吸烟/以前吸烟者」：严重异生（stage 6）原位癌（stage 7） | 「从未吸烟者」的正常样本 |
fourth step:「SCC」 | 「吸烟/以前吸烟者」：肺鳞状细胞癌（stage 8） | 「从未吸烟者」的正常样本 |
我们从GEO数据库的series_matrix中提取临床信息
GSE_number<-"GSE33479"
load(paste0(GSE_number,"_agilent_MAList.Rdata"))
#1.3 根据临床信息获得实验设计获得分组信息，这部分自己做......
library(GEOquery)
library(stringr)
getGEO(filename = paste0(GSE_number,"_series_matrix.txt.gz"),getGPL = F)->geoObject
#提取临床信息
pData(geoObject)->pd
#这个数据集的实验设计中M值矩阵虽然不是表达矩阵，但是按照limmauserguide 16章
#对双通道芯片的处理，可以作为差异表达分析（lmFit,ebayes那套的输入）
MA$M->dat
head(dat)
#检查下pd的每行和exp的每列样本是否对应
identical(rownames(pd),str_split(str_split(colnames(dat),"_",simplify = T)[,1],"/",simplify = T)[,3])
#提取有用的临床信息
phe<-data.frame(pd$geo_accession,
pd$`phenotypical stage:ch1`,
pd$characteristics_ch1.2,
pd$characteristics_ch1.3,
pd$characteristics_ch1.4,
pd$characteristics_ch1.5,
pd$characteristics_ch1.6)
#按文献中的九个stage获取分组信息
group_list<-phe$pd..phenotypical.stage.ch1.
table(group_list)
#1.4 探针id转换
head(dat)
#行没有探针或基因信息，要从MAList的genes里面提取
head(MA$genes,20)
#两种方法1.是使用probeid列，转换成基因名称，之前有很多教程了这里不做介绍
MA$genes$ProbeName->probeid
#2.是直接使用GeneName列（如果MA$genes的数据框中包含该信息的话）
#这个数据集可以这么做的，但是需要手动去掉"EQC","DarkCorner",
#"NegativeControl","GE_BrightCorner","E1A"等对于实验设计无生物学意义的点，
#如果你充分理解这4.5w个探针谁该留谁该去的话，也可以这么做......
MA$genes$GeneName->genename
#这里只展示法1
dat->dat1
probeid->rownames(dat1)
library(AnnoProbe)
#获取平台GPL编号
GPL<-geoObject@annotation
#我们的AnnoProbe包，返回一列是probe id,一列是symbol的数据框
ids=idmap(GPL,'soft')
head(ids)
#检查我们MA$genes得到的M矩阵的列名是否都在ids$ID中
dat1=dat1[rownames(dat1) %in% ids$ID,]
#将ids的行顺序按probe_id整理至与dat1的行名一致
#反正结果是ids每行的probe_id与dat的rownames一致
ids=ids[match(rownames(dat1),ids$ID),]
head(ids)
#以下去重方法是保留相同基因名称中中位数最大的
ids$median=apply(dat1,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果
dat1=dat1[ids$ID,] #新的ids取探针id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat1)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat1[1:4,1:4]
save(dat1,group_list,phe,file = 'step1-output.Rdata')
这样我们保存的M矩阵dat1，分组信息和之后可能用到的临床信息为'step1-output.Rdata'
四、差异表达分析
差异表达分析之前需要画PCA图和样本相关性热图或方差前1000基因热图，对质量进行粗略检查
因为这部分代码之前分享很多，也比较基础，所以为不占用篇幅就不过多展示了。
#这里我们按照上文表格把小分组整理成大分组
table(group_list)
group_list[group_list=="normal normofluorescent"|group_list=="normal hypofluorescent"|group_list=="hyperplasia"]<-"normal_bronchial_tissue"
group_list[group_list=="metaplasia"|group_list=="mild dysplasia"|group_list=="moderate dysplasia"]<-"low_grade_lesions"
group_list[group_list=="severe dysplasia"|group_list=="carcinoma in situ"]<-"high_grade_lesions"
group_list[group_list=="squamous cell carcinoma"]<-"squamous_cell_carcinoma"
table(group_list)
#然后画个PCA图进行检查
#这里PCA图我使用了自己在写的一个R包easyArray，因为这部分代码之前分享很多，也比较基础，所以为不占用篇幅就不过多展示了，等成熟了也许之后会分享出来给大家使用
library(easyArray)
checkPCA(dat1,group_list)
#会发现结果不太理想，因为这个实验做的是肺鳞癌支气管的癌症发生过程，这四个分组本身对应癌症发生过程从早到晚一条时间线上的四个阶段。没有癌症vs正常那么大的差异也是比较正常的。
#将分组弄成因子，然后走limma差异表达分析
group_list<-factor(group_list,levels = c("normal_bronchial_tissue","low_grade_lesions","high_grade_lesions","squamous_cell_carcinoma"))
design <- model.matrix(~ 0+group_list)
colnames(design)<-levels(group_list)
library(limma)
fit<-lmFit(dat1,design)
#这里随便设置三个差异表达分析的比较，分别是SCC比正常样本，高水平损伤比正常，低水平损伤比正常
contrast.matrix <- makeContrasts("squamous_cell_carcinoma-normal_bronchial_tissue",
"high_grade_lesions-normal_bronchial_tissue",
"low_grade_lesions-normal_bronchial_tissue",
levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
#从fit2中提取差异表达分析的结果，coef是选择第几个比较，adjust是p值校正方法，
#默认为BH法，n是选取p.adjust从小到大排列的前几个基因
SCC_normal_DEG<-na.omit(topTable(fit2, coef=1, adjust="BH",n=Inf))
high_normal_DEG<-na.omit(topTable(fit2, coef=2, adjust="BH",n=Inf))
low_normal_DEG<-na.omit(topTable(fit2, coef=3, adjust="BH",n=Inf))
#定义一个函数，用于给limma生成的差异表达分析结果加上up,down，以及画图需要的-log10P.adjust
addchange<-function(deg){
logFC_t <- with(deg, mean(abs(deg$logFC)) + 2 * sd(abs(deg$logFC)))
deg$v = -log10(deg$adj.P.Val)
deg$g = ifelse(deg$adj.P.Val > 0.05, "stable",
ifelse(deg$logFC > logFC_t, "up",
ifelse(deg$logFC < -logFC_t, "down", "stable")))
deg$name = rownames(deg)
return(deg)
}
addchange(SCC_normal_DEG)->SCC_normal_DEG
addchange(high_normal_DEG)->high_normal_DEG
addchange(low_normal_DEG)->low_normal_DEG
#结果是三个有上下调和p.adjust和基因名的数据框
#这里还是使用了笔者自己写的easyArray包来进行火山图绘制，结果如下：
DEG_volcano(SCC_normal_DEG)
DEG_volcano(high_normal_DEG)
DEG_volcano(low_normal_DEG)
结果还行
我们这篇推文就到此结束了，内容主要是agilent双色芯片的原始数据读入和预处理，
下周周五更新这篇文章的下集，准备给大家复现这篇文章的WGCNA及富集分析的结果。