100个GEO基因表达芯片或转录组数据处理之GSE126848（003）

写在前边

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

数据信息检索

可以看到GSE126848是转录组高通量测序数据，因此可以使用GEOquery包下载数据临床信息，并且手动下载表达矩阵并整理

使用GEOquery包下载数据

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

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

geo_accession <- "GSE126848"
gset <- GEOquery::getGEO(geo_accession, destdir = "./", AnnotGPL = F, getGPL = F)
eSet <- gset[[1]]
gpl <- eSet@annotation

处理表型数据

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

pdata <- pData(eSet)

geo_accession	description	disease:ch1	gender:ch1	tissue:ch1
GSM3615293	2683	NAFLD	Male	Liver
GSM3615294	2685	NAFLD	Male	Liver
GSM3615295	2687	NAFLD	Male	Liver
GSM3615296	2689	NAFLD	Female	Liver
GSM3615297	2691	NAFLD	Female	Liver
GSM3615298	2693	NAFLD	Male	Liver

pdata %<>%
    dplyr::mutate(
        Sample = geo_accession,
        Group = case_when(`diagnosis:ch1` == "HC" ~ "Control", `diagnosis:ch1` == "NASH" ~ "Case", TRUE ~ NA),
        Age = `age (y):ch1`,
        Sex = str_to_title(`gender:ch1`),
        Stage = `fibrosis (stage):ch1`
    ) %>%
    dplyr::filter(!is.na(Group)) %>%
    dplyr::select(Sample, Group, Age, Sex)
fwrite(pdata, file = str_glue("{geo_accession}_pdata.csv"))

处理表达谱数据

原始数据为Count值，需要标准化为TPM，并且基因名是Ensembl ID转换为Symbol基因名，可以使用到我自己写的几个函数genekit、bioquest；有需要可以联系我的公众号@恩喜玛生物，加入交流群

import pandas as pd
import genekit as gk
import bioquest as bq

fdata = pd.read_csv("GSE126848_Gene_counts_raw.txt.gz",sep='\t',index_col=0)
pdata = pd.read_csv("GSE126848_pdata.csv",index_col=0)
pdata.drop(columns=["Sample2"]).to_csv("GSE126848_pdata.csv")

fdata与pdata样本名统一，这里使用了Python的字符串格式化方法

fdata = fdata.loc[:,["{0:0>4}".format(x) for x in pdata.Sample2]]
fdata.columns = pdata.index.to_list()

保存一份原始Count数据信息

fdata.to_csv("GSE126848_count.csv.gz")

Count 转 TPM

fdata = gk.countto(fdata, towhat='tpm', geneid='Ensembl', species='Human')

Ensembl ID转换为Symbol基因名

fdata=gk.geneIDconverter(
    frame=fdata,
    from_id='Ensembl',
    to_id='Symbol',
    keep_from=False,
    gene_type=False,
    )

去重复

根据每个基因表达量的中位数去除重复的基因

fdata=bq.tl.unique_exprs(fdata)

保存TPM基因表达量数据

fdata.to_csv("GSE126848_tpm.csv.gz")