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