100个GEO基因表达芯片或转录组数据处理之GSE26899(008)
写在前边
虽然现在是高通量测序的时代,但是GEO、ArrayExpress等数据库储存并公开大量的基因表达芯片数据,还是会有大量的需求去处理芯片数据,并且建模或验证自己所研究基因的表达情况,芯片数据的处理也可能是大部分刚学生信的道友入门R语言数据处理的第一次实战,因此准备更新100个基因表达芯片或转录组高通量数据的处理。
数据信息检索
可以看到GSE26899是基因表达芯片数据,因此可以使用GEOquery包下载数据
使用GEOquery包下载数据
using(tidyverse, GEOquery, magrittr, data.table, AnnoProbe, clusterProfiler, org.Hs.eg.db, org.Mm.eg.db)
注:using
是我写的函数,作用是一次性加载多个R包,不用写双引号,并且不在屏幕上打印包的加载信息,可以参考之前的推文《如何优雅地管理R包》using的定义;函数名字using是在模仿Julia语言中的包加载函数
因为文件太大,在R内下载失败,可通过图片中的方法下载文件,并且把下载后的文件放在destdir = "./"即当前工作目录,GEOquery::getGEO便能跳过下载,直接使用本地的文件。
geo_accession <- "GSE26899"
gset <- GEOquery::getGEO(geo_accession, destdir = "./", AnnotGPL = F, getGPL = F)
eSet <- gset[[1]]
gpl <- eSet@annotation
处理表型数据
这部分是很关键的,可以筛选一下分组表型信息,只保留自己需要的样本,在这里只保留tissue:ch1中 normal gastric tissue和gastric cancer tissue的样本,作为后续分析的样本(根据自己的研究目的筛选符合要求的样本)
pdata <- pData(eSet)
geo_accession | patient:ch1 | tissue:ch1 | age:ch1 | ajcc stage:ch1 | gender:ch1 | lauren classification:ch1 | location:ch1 |
---|---|---|---|---|---|---|---|
GSM662374 | KG004C | Gastric tumor tissue | 41 | 4 | M | Diffuse | Body |
GSM662375 | KG008N | Gastric Surrounding normal tissue | |||||
GSM662376 | KG009C | Gastric tumor tissue | 55 | 1 | M | Intestinal | Fundus |
GSM662377 | KG012C | Gastric tumor tissue | 54 | M | |||
GSM662378 | KG013C | Gastric tumor tissue | 45 | 3 | F | Intestinal | Antrum |
GSM662379 | KG014N | Gastric Surrounding normal tissue |
pdata %<>%
dplyr::mutate(
Sample = geo_accession,
Patient=`patient:ch1`,
Age=`age:ch1`,
Stage=`ajcc stage:ch1`,
Sex=`gender:ch1`,
Histology=`lauren classification:ch1`,
Location=`location:ch1`,
Group = dplyr::case_when(
stringr::str_detect(`tissue:ch1`, "normal") ~ "Nomal",
stringr::str_detect(`tissue:ch1`, "tumor tissue") ~ "Cancer",
TRUE ~ NA),
) %>%
dplyr::filter(!is.na(Group) %>%
dplyr::select(Sample,Patient, Group,Histology, Age,Stage, Location, Sex)
处理表达谱数据
数据大小不大于50不需要取log
exprs_mtx <- exprs(eSet)
range(exprs_mtx, na.rm = TRUE)
# 5.11 15.227
probe_exprs <- as.data.table(exprs_mtx, keep.rownames = "ProbeID")
探针与基因Symbol对应关系
从AnnoProbe包中获取探针与GeneID对应关系
probe2symbol <- AnnoProbe::idmap(gpl = gpl, type = "soft", mirror = "tencent", destdir = od) %>%
dplyr::rename(ProbeID = ID, GeneID = symbol) # pipe", "bioc", "soft"
ID转换
把表达矩阵中的探针名转换为基因名;transid是我写的一个R函数,有需要可以联系我的公众号@恩喜玛生物,加入交流群
fdata <- transid(probe2symbol, probe_exprs)
保存数据
common_samples <- base::intersect(colnames(fdata),pdata$Sample)
fdata %<>% select(all_of(c("GeneID",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"))