100个GEO基因表达芯片或转录组数据处理
写在前边
虽然现在是高通量测序的时代，但是GEO、ArrayExpress等数据库储存并公开大量的基因表达芯片数据，还是会有大量的需求去处理芯片数据，并且建模或验证自己所研究基因的表达情况，芯片数据的处理也可能是大部分刚学生信的道友入门R语言数据处理的第一次实战，因此准备更新100个基因表达芯片或转录组高通量数据的处理。
数据信息检索
可以看到GSE21257是基因表达芯片数据，因此可以使用GEOquery包处理数据
image.png
使用GEOquery包下载数据
using(tidyverse, GEOquery, magrittr, data.table, AnnoProbe, clusterProfiler, org.Hs.eg.db, org.Mm.eg.db)
注：using是我写的函数，有需要可以；using作用是一次性加载多个R包，不用写双引号，并且不在屏幕上打印包的加载信息
因为文件太大，在R内下载失败，可通过图片中的方法下载文件，GEOquery::getGEO直接读取本地的文件。
geo_accession <- "GSE21257"
eSet <- getGEO(filename='GSE21257_series_matrix.txt.gz', AnnotGPL = F, getGPL = F)
gpl <- eSet@annotation
image.png
image.png
处理表型数据
这部分是很关键的，可以筛选一下分组表型信息，只保留自己需要的样本，根据自己的研究目的筛选符合要求的样本，筛选临床信息。
pdata <- pData(eSet)
pdata %<>%
dplyr::mutate(
Sample = geo_accession,
Age = (str_remove(`age:ch1`,' months') %>% as.numeric()) / 12,
Gender = `gender:ch1`,
HistologicalSubtype=`histological subtype:ch1`,
Grade=`huvos grade:ch1`,
TumorLocation=`tumor location:ch1`,
OS.Status=ifelse(str_detect(`status:ch1`,'Alive'),0,1),
OS.Time=(str_split(`status:ch1`, ' ',simplify = T) %>% .[,3] %>% as.numeric()) * 12
) %>%
dplyr::select(Sample,OS.Status,OS.Time,Age,Gender,HistologicalSubtype,Grade,TumorLocation,everything())
处理表达谱数据
数据大小不大于50不需要取log
exprs_mtx <- exprs(eSet)
if(max(exprs_mtx, na.rm = TRUE)<50 | min(exprs_mtx, na.rm = TRUE)<0){
message("基因表达最大值小于50或者最小值小于0不需要log转化")
}else {
message("基因表达最大值大于50需要log转化")
exprs_mtx <- log2(exprs_mtx+1)
}
probe_exprs <- as.data.table(exprs_mtx, keep.rownames = "ProbeID")
探针与基因Symbol对应关系
从AnnoProbe包中获取探针与GeneID对应关系
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL10295[图片上传失败...(image-f40cc5-1720142997077)]
[图片上传失败...(image-71cb88-1720142997077)]
idmaps <- function(ann_file, ProbeID = "ID", Feature = "Symbol", skip = 229, pattern = "control") {
temp <- fread(ann_file, skip = skip, nThread = 8)
vars <- c(ProbeID, Feature)
temp <- temp[, ..vars]
data.table::setnames(temp, c("ProbeID", "Feature"))
temp <- temp[!is.null(Feature), ][!is.na(Feature), ][Feature != "", ][!stringr::str_detect(string = Feature, pattern = pattern), ]
return(as.data.frame(temp))
}
probe2symbol <- idmaps("GPL10295_family.soft.gz", Feature = "Symbol", skip = 413)
ID转换
把表达矩阵中的探针名转换为基因名；transid是我写的一个R函数，有需要可以
fdata <- transid(probe2symbol, probe_exprs)
保存数据
common_samples <- base::intersect(colnames(fdata),pdata$Sample)
fdata %<>% select(all_of(c("Feature",common_samples)))
fwrite(fdata, file = stringr::str_glue("{geo_accession}_{gpl}_fdata.csv.gz"))
pdata %<>% dplyr::filter(Sample %in% common_samples)
fwrite(pdata, file = stringr::str_glue("{geo_accession}_{gpl}_pdata.csv"))
