100个GEO基因表达芯片或转录组数据处理之GSE166193-GPL16686平台（013）

写在前边

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

数据信息检索

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

使用GEOquery包下载数据

using(tidyverse, GEOquery, magrittr, data.table, AnnoProbe, clusterProfiler, org.Hs.eg.db, org.Mm.eg.db)

注：using是我写的函数，作用是一次性加载多个R包，不用写双引号，并且不在屏幕上打印包的加载信息，可以参考之前的推文《如何优雅地管理R包》using的定义；函数名字using是在模仿Julia语言中的包加载函数

直接用GEOquery包下载数据集很多时候会网络中断导致下载失败，因此可以在网页直接下载GSEXXX_series_matrix.txt.gz文件，然后用getGEO函数读入。

geo_accession <- "GSE166193"
eSet <- getGEO(filename = str_glue("{geo_accession}_series_matrix.txt.gz"))
gpl <- eSet@annotation

处理表型数据

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

pdata <- pData(eSet)

pdata %<>%
    dplyr::mutate(
        Sample = geo_accession,
        Group = ifelse(str_detect(`treatment:ch1`,"without"),"Control","Treat")
    ) %>%
    dplyr::filter(!is.na(Group)) %>% 
    dplyr::select(Sample,Group,everything())

处理表达谱数据

exprs_mtx <- exprs(eSet)
if(max(exprs_mtx, na.rm = TRUE)<50){
    message("基因表达最大值小于50不需要log转化")
}else {
    message("基因表达最大值大于50需要log转化")
    exprs_mtx <- log2(exprs_mtx+1)
}
probe_exprs <- as.data.table(exprs_mtx, keep.rownames = "ProbeID")

探针与基因Symbol对应关系

下载GPL16686_family.soft注释文件https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL16686

soft=getGEO("GPL16686_family.soft")
a=soft@dataTable@table %>% 
    dplyr::mutate(ProbeID=ID) %>% 
    dplyr::select(ProbeID,GB_ACC) %>% 
    dplyr::filter(GB_ACC!="",GB_ACC!="-")

从https://ftp.ncbi.nih.gov/gene/DATA/下载GB_ACC与基因名的对应关系gene2accession.gz文件

b=fread("gene2accession.gz",select = c("#tax_id","Symbol","RNA_nucleotide_accession.version")) %>%
    dplyr::filter(`#tax_id`==9606) %>% 
    dplyr::mutate(Feature=Symbol,GB_ACC=str_remove(RNA_nucleotide_accession.version,"\\.\\d*")) %>% 
    dplyr::select(GB_ACC,Feature) %>% 
    dplyr::filter(GB_ACC!="",GB_ACC!="-",Feature!="",Feature!="-")

得到注释文件

probe2symbol=merge(a,b,by="GB_ACC") %>% select(ProbeID,Feature) %>% distinct(.keep_all = T)

ID转换

把表达矩阵中的探针名转换为基因名；transid是我写的一个R函数，有需要可以扫码加入交流群

fdata <- transid(probe2symbol, probe_exprs)

保存数据

common_samples <- base::intersect(colnames(fdata),pdata$Sample)
fdata %<>% select(Feature,all_of(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"))