100个GEO基因表达芯片或转录组数据处理GSE60542（012）

最新推荐文章于 2024-09-13 11:18:53 发布

探因科技

最新推荐文章于 2024-09-13 11:18:53 发布

阅读量421

点赞数 3

文章标签：后端

本文链接：https://blog.csdn.net/weixin_44493991/article/details/140897450

版权

写在前边

虽然现在是高通量测序的时代，但是GEO、ArrayExpress等数据库储存并公开大量的基因表达芯片数据，还是会有大量的需求去处理芯片数据，并且建模或验证自己所研究基因的表达情况，芯片数据的处理也可能是大部分刚学生信的道友入门R语言数据处理的第一次实战，因此准备更新100个基因表达芯片或转录组高通量数据的处理。

数据信息检索

可以看到GSE60542是基因表达芯片数据，因此可以使用GEOquery包下! alt

使用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 <- "GSE60542"
eSet <- getGEO(filename=stringr::str_glue('{geo_accession}_series_matrix.txt.gz'), AnnotGPL = F, getGPL = F)
gpl <- eSet@annotation

处理表型数据

这部分是很关键的，可以筛选一下分组表型信息，只保留自己需要的样本，作为后续分析的样本（根据自己的研究目的筛选符合要求的样本）

title	geo_accession	age.at.operation:ch1	histological..subtype:ch1	Sex:ch1	Stage:ch1	tnm:ch1	tumor.in.left.right.lobe:ch1	vascular.invasion:ch1
23-N,Normal thyroid,N0	GSM1481838	52	Trabecular	F	III	pT3N0Mx	right	no
1-PT-N1,Papillary thyroid carcinoma,N1	GSM1481839	35	Papillary	M	I	pT3N1bMx	left	no
9-LNM,Lymph node metastasis,N1	GSM1481840	43	Papillary	M	I	pT1N1bMx	right	no
24-PT-N0,Papillary thyroid carcinoma,N0	GSM1481841	29	Trabecular	F	I	pT1N0Mx	left	no
23a-N,Normal thyroid,N0	GSM1481842	52	Trabecular	F	III	pT3N0Mx	right	no

pdata <- pData(eSet)

pdata %<>%
    dplyr::mutate(
        Sample = geo_accession,
        SampleType=ifelse(str_detect(title,"Normal"),"Normal","Tumor"),
        Age=`age.at.operation:ch1`,
        HistologicalType=`histological..subtype:ch1`,
        Sex=`Sex:ch1`,
        Stage=`Stage:ch1`,
        TNM_Stage=`tnm:ch1`,
        TumorSite=`tumor.in.left.right.lobe:ch1`,
        VascularInvasion=`vascular.invasion:ch1`
    ) %>%
    dplyr::select(Sample,SampleType,Age,HistologicalType,Sex,Stage,TNM_Stage,TumorSite,VascularInvasion) %>% 
    dplyr::filter(!str_detect(TNM_Stage,"M1$")) # 删除转移的肿瘤样本

处理表达谱数据

数据大小不大于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对应关系

probe2symbol <- AnnoProbe::idmap(gpl = gpl, type = "soft", mirror = "tencent", destdir = od) %>%
  dplyr::rename(ProbeID = ID, Feature = symbol) # pipe", "bioc", "soft"

处理一个探针对应多个基因的情况，保留全部基因

probe2symbol %<>% tidyr::separate_longer_delim(cols = Feature, delim = " /// ") %>% 
    dplyr::filter(!is.na(Feature),Feature!="")

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"))

本文由 mdnice 多平台发布