从ArrayExpress数据库下载一个芯片做注释

---

前言

最近课题需要一组支撑数据，Characterization of a Large Panel of Patient-Derived Tumor
Xenografts Representing the Clinical Heterogeneity of Human Colorectal Cancer 里面的EMTAB-
991，是基因表达的矩阵。

---

# 首先加载必要的包
library(ArrayExpress)
library(oligo)

二、下载芯片

这次需要的芯片是"E-TAM-991”，来自Characterization of a Large Panel of Patient-Derived Tumor Xenografts Representing the Clinical Heterogeneity of Human Colorectal Cancer. Clin Cancer Res 18, 5314–5328 (2012). 也就是保存在arrayexpress数据库。直接搜E-MTAM-991，可能会搜不到，但是不要紧，https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-991/，直接跳转页面。

这些Files就是我们需要的，先下载下来吧。

1.用getAE()下载

dir.create("E-MTAB-991")  # 建一个文件夹，把E-TAM-911的raw data下载到里面。
dir <- "E-MTAB-991"
anno_AE <- getAE("E-MTAB-991", path = dir, type = "raw")

getAE() 里面两个参数比较重要：1. type: raw表示只下载raw文件，processed表示只下载process文件，full表示都下载。2. extract: 默认TRUE。如果选FALSE，下载的zip文件不会自动解压缩。

2.用wget 下载

这种就比较麻烦了，但是有时候速度快。

wget https://www.ebi.ac.uk/arrayexpress/files/E-MTAB-991/E-MTAB-991.sdrf.txt
wget https://www.ebi.ac.uk/arrayexpress/files/E-MTAB-991/E-MTAB-991.sdrf.txt
wget https://www.ebi.ac.uk/arrayexpress/files/E-MTAB-991/E-MTAB-991.raw.1.zip
wget https://www.ebi.ac.uk/arrayexpress/files/E-MTAB-991/E-MTAB-991.raw.2.zip
wget https://www.ebi.ac.uk/arrayexpress/files/E-MTAB-991/E-MTAB-991.raw.3.zip
wget https://www.ebi.ac.uk/arrayexpress/files/E-MTAB-991/A-AFFY-113.adf.txt

代码如下（示例）：

二、读取芯片数据

# 把样本对应的信息先读进去
raw_data_dir <- "E-MTAB-991"
sdrf_location <- file.path(raw_data_dir, "E-MTAB-991.sdrf.txt")
SDRF <- read.delim(sdrf_location)
rownames(SDRF) <- SDRF$Array.Data.File
SDRF <- AnnotatedDataFrame(SDRF)

然后注意到一个问题，我们下载的raw data是压缩包，需要手动解压缩。

library(R.utils) 
unzip("E-MTAB-991/E-MTAB-991.raw.1.zip",exdir = "E-MTAB-991",overwrite = F)
unzip("E-MTAB-991/E-MTAB-991.raw.2.zip",exdir = "E-MTAB-991",overwrite = F)
unzip("E-MTAB-991/E-MTAB-991.raw.3.zip",exdir = "E-MTAB-991",overwrite = F)

解压缩完成之后读入数据。

raw_data <- oligo::read.celfiles(filenames = file.path(raw_data_dir, 
                                                       SDRF$Array.Data.File),
                                 verbose = FALSE, phenoData = SDRF)
......                                 
Reading in : E-MTAB-991/CGCLG_030810_HT_HG-U133A_PHY03081001D05_1.CEL
Reading in : E-MTAB-991/CGCLG_030810_HT_HG-U133A_PHY03081001F11_1.CEL
Reading in : E-MTAB-991/CGCLG_061110_HT_HG-U133A_PHY06111001H04_1.CEL
Reading in : E-MTAB-991/CGCLG_030810_HT_HG-U133A_PHY03081001G11_1.CEL
Reading in : E-MTAB-991/CGCLG_030810_HT_HG-U133A_PHY03081001H11_1.CEL
......
Reading in : E-MTAB-991/CGCLG_030810_HT_HG-U133A_PHY03081001C12_1.CEL
Reading in : E-MTAB-991/CGCLG_030810_HT_HG-U133A_PHY03081001D12_1.CEL
Reading in : E-MTAB-991/CGCLG_061110_HT_HG-U133A_PHY06111001C05_1.CEL
Reading in : E-MTAB-991/CGCLG_061110_HT_HG-U133A_PHY06111001D05_1.CEL
Warning message:
In oligo::read.celfiles(filenames = file.path(raw_data_dir, SDRF$Array.Data.File),  :
  'channel' automatically added to varMetadata in phenoData.

读入以后做一个简单的验证

stopifnot(validObject(raw_data))

该处使用的url网络请求的数据。

pd<- Biobase::pData(raw_data) # pd里面是样本信息
# Robust Multi-array Average (RMA) 算法初处理
exp_raw<- oligo::rma(raw_data, normalize=TRUE, subset=NULL) 
exp_rma<- Biobase::exprs(exp_raw) # 我们的表达矩阵出现了
save(raw_data,file='raw.Rdata')
save(pd,file='pd.Rdata')
save(exp_rma,file='rma.Rdata')
write.csv(exp_raw,"exp_raw.csv")

三、注释芯片数据

下一步，就是做注释了。

name <- SDRF@data
table(name$Comment.Array.serial.number.)
Affy HG-U133A 
          115
# Affy HG-U133A 对应的包是hg133a.db 也就是 GPL571 
getGEO("GPL571", destdir="E-MTAB-991")

当然也可以通过https://www.ncbi.nlm.nih.gov/search/all/?term=GPL571直接下

rm(list = ls())
load("rma.Rdata")
gpl <- read.csv("E-MTAB-991/GPL571-17391.txt",comment.char = "#",sep = "\t")
namel <- gpl[,c("ID",colnames(gpl)[11])]
table(gpl$ID %in% row.names(exp_rma))
 TRUE 
22277 
mt <-merge(namel,exp_rma,by.x = "ID",by.y = "row.names") 
write.csv(mt,"E-MTAB-991_exp.csv")

---