Affy 包处理 affymetrix 表达谱芯片

约 858 字大约 3 分钟

2024-10-21

Questions

从Gene Expression Omnibus (GEO) 查询GSE65496，并下载GSE65496_family.soft和原始数据(.cel格式)
仔细研读soft格式文件
试着用Bioconductor的affy包下函数处理.cel文件 (可以网页搜索，或者上Bioconductor官网查询)，得到探针表达值
下载该数据对应的平台数据，通过文本处理，将探针和gene symbol对应，得到gene symbol的表达值
对表达谱数据做质量分析（质控）：①灰度图；②RNA降解分析
对芯片数据做如下分析：①Mas5或rma方法进行背景校正，标准化（归一化）；②提取表达值

Codes

# Rstudio Server URL: http://202.195.187.9:8787/
# 账号密码与Linux一致
# 切换R路径
.libPaths("/home/biotools/anaconda3/envs/R4.3.1/lib/R/")

# 工作路径
setwd(dir = "/public/workspace/stu22230111/12_R/01_BioChips/Homework_01/")

# 加载R包 (Biotools下已安装好相关R包，无需重复安装)
library("affy")
library("GEOquery")

###################################
## 问题1：下载GSE65496_family.soft和cel原始数据
# (虽然直接使用R包GEOquery比较方便，但下载速度比较慢，还是建议使用wget手动下载)

# 下载soft文件
# wget "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE65nnn/GSE65496/soft/GSE65496_family.soft.gz" -O "GSE65496.soft.gz"
getGEOfile("GSE65496", destdir="./RAW/")

# 下载cel原始数据
# wget "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE65496&format=file" -O "GSE65496_RAW.tar"
getGEOSuppFiles("GSE65496")

# 解包cel原始数据
untar("./GSE65496/GSE65496_RAW.tar", exdir = "./Data/")

###################################
## 问题2：研读soft文件
# NCBI官方文档：https://www.ncbi.nlm.nih.gov/geo/info/soft.html
soft_file = getGEO(filename = "./RAW/GSE65496.soft.gz")

###################################
## 问题3：用affy包处理.cel文件，得到探针表达值
# 读取CEL文件
cel_path = "./Data/"
cel_data = ReadAffy(celfile.path = cel_path)

##################################
# 画灰度图
# pdf记录矢量信息，保存和打开时很慢，需耐心等待一会儿
# 或者使用png (png("result.png", weight = 1920, height = 1080))
pdf("result.pdf")

# 两个一行，一共四行
par(mfrow = c(4,2))
for(i in 1:8){
  image(cel_data[,i])
}

dev.off()

##################################
# RNA降解分析
# 使用 RColorBrewer 中的调色板生成高对比度颜色
library(RColorBrewer)

colors = brewer.pal(n = length(sampleNames(cel_data)), name = "Dark2")

pdf("degradation.pdf", width = 10, height = 7)

RNAdeg = AffyRNAdeg(cel_data)

plotAffyRNAdeg(RNAdeg, cols = colors)

legend("topleft",
       sampleNames(cel_data), 
       col = colors, 
       lwd=1, 
       inset = 0.05,
       cex = 0.5)

dev.off()

#################################
# rma与mas5方法进行背景校正，标准化（归一化），得到探针表达值
# Bug: 需重新编译安装一下preprocessCore，否则无法使用rma
# 解决方法来源：https://github.com/Bioconductor/bioconductor_docker/issues/22
# BiocManager::install("preprocessCore", configure.args = c(preprocessCore = "--disable-threading"), force= TRUE, update=TRUE, type = "source")
# 记得重新打开R/Rstudio哦
# 这里我在 biotools 下的R里面已经重新编译安装好了

eset_rma = rma(cel_data)
eset_mas5 = mas5(cel_data)

# 提取表达值 (两种算法)
expr_rma = exprs(eset_rma)
expr_mas5 = exprs(eset_mas5)

# 保存表达值 (两种算法)
write.csv(expr_mas5, file = "expr_mas5.csv")
write.csv(expr_rma, file = "expr_rma.csv")

##################################
## 问题4：将探针和gene symbol对应，得到gene symbol的表达值
# 下载对应的平台数据(GPL570)
platform_data = getGEO("GPL570",destdir="./RAW/")

# 载入dplyr
library("dplyr")

## rma算法
# 提取ID与Gene Symbol的对应表格
gene_symbol_ID = Table(platform_data)[c(1,11)]
# 按照ID添加对应的Gene Symbol
expr_rma = as.data.frame(expr_rma)
expr_rma$ID = rownames(expr_rma)
# 内连接
processed_rma = inner_join(gene_symbol_ID, expr_rma, by="ID")
# 去除ID列
processed_rma = processed_rma[,-1]
write.table(processed_rma,"GeneSymbol_rma.csv",sep=",")

## mas5算法
# 提取ID与Gene Symbol的对应表格
gene_symbol_ID = Table(platform_data)[c(1,11)]
# 按照ID添加对应的Gene Symbol
expr_mas5 = as.data.frame(expr_mas5)
expr_mas5$ID = rownames(expr_mas5)
# 内连接
processed_mas5 = inner_join(gene_symbol_ID, expr_mas5, by="ID")
# 去除ID列
processed_mas5 = processed_mas5[,-1]
write.table(processed_mas5,"GeneSymbol_mas5.csv",sep=",")

Files

RMA结果文件下载

📄 GeneSymbol_rma.csv

MAS5结果文件下载

📄 GeneSymbol_mas5.csv