100个GEO基因表达芯片或转录组数据处理025 GSE248467

写在前边

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

数据信息检索

可以看到GSE248467是 高通量测序数据，因此可以使用GEOquery包下临床信息，直接在页面下载 NCBI 生产的转录组表达矩阵 counts 和 tpm 两种进一步处理

使用GEOquery包下载 临床 数据

BiocManager::install('ScienceAdvances/Canton')
Canton::using(using, tidyverse, GEOquery, magrittr, data.table, AnnoProbe, clusterProfiler, org.Hs.eg.db, org.Mm.eg.db,ggdendro,ComplexHeatmap)

注：using作用是一次性加载多个R包，不用写双引号，并且不在屏幕上打印包的加载信息

处理表型数据

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

geo_accession <- "GSE248467"
eSet <- GEOquery::getGEO(geo_accession, AnnotGPL = F, getGPL = F)

pdata=Biobase::pData(eSet[[1]]) %>% 
    dplyr::mutate(Sample = geo_accession,Group=stringr::str_split(pdata$title,"_",simplify=T) %>% .[,2]) %>% 
    dplyr::select(Sample,Group) %T>% 
    fwrite(file = stringr::str_glue("{geo_accession}_pdata.csv"))

处理表达谱数据

annot <- fread('Human.GRCh38.p13.annot.tsv.gz') %>% 
    dplyr::select(GeneID, Symbol )

counts=fread(glue::glue('{geo_accession}_norm_counts_TPM_GRCh38.p13_NCBI.tsv.gz')) %>% 
    merge(annot,.,by='GeneID') %>% dplyr::select(-GeneID) %>% unique_exprs %T>% 
    fwrite(file = stringr::str_glue("{geo_accession}_counts.csv.gz"))

tpm=fread(glue::glue('{geo_accession}_norm_counts_TPM_GRCh38.p13_NCBI.tsv.gz')) %>% 
    merge(annot,.,by='GeneID') %>% dplyr::select(-GeneID) %>% unique_exprs %T>% 
    fwrite(file = stringr::str_glue("{geo_accession}_tpm.csv.gz"))

数据质控

qcplot为自定义函数，作用是绘制用于质控判断的图，如PCA、top20基因热图、树状图，PCA图可以三组区别还是比较明显的

qcplot(
  data=tpm %>% column_to_rownames('Symbol'), 
  pdata=pdata$Group, 
  w=12,h=12
)