100个GEO基因表达芯片或转录组数据处理(15) GSE98895 GPL6947

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

写在前边

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

数据信息检索

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

使用GEOquery包下载数据

安装所需R包

BiocManager::install("lumi")
BiocManager::install("lumiHumanIDMapping")

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

注：using是我写的函数，有需要可以联系我，加入交流群；using作用是一次性加载多个R包，不用写双引号，并且不在屏幕上打印包的加载信息
因为文件太大，在R内下载失败，可通过图片中的方法下载文件，GEOquery::getGEO直接读取本地的文件。

geo_accession <- "GSE23317"
eSet <- getGEO(filename=stringr::str_glue('{geo_accession}_series_matrix.txt.gz'), AnnotGPL = F, getGPL = F)

处理表型数据

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

pdata <- pData(eSet)

pdata %<>%
    dplyr::mutate(
        Sample = geo_accession,
        Group = ifelse(`disease state:ch1`=='control','Control','MetabolicSyndrome')
    ) %>%
    dplyr::select(Sample,Group,everything())

处理表达谱数据

这里需要处理Illumina表达芯片的原始数据，因为GSE98895_series_matrix.txt.gz 中的基因表达值从-20到几千，无法用于后续分析，因此我们自己标准化处理.
下载原始数据
!


lumi代码主要两个步骤，读入数据和lumiExpresso（包括lumiB调整芯片背景；lumiT 对数据进行方差固定；lumiN 对方差固定后的数据进行标准化，lumiQ 评估数据质量。所有这些方法构成了一个预处理流程），可以参考lumi包发表的文章 PMID: 18467348



读入原始数据，并重命名行名，40个样本这里需要根据上面代码返回结果自己改动，前两行是探针id和对应的symbol，后面每两列对应一个样本，是荧光信号强度和pvalue
最后保存原始数据到tmp.txt文件中(参考生信技能树代码)

a <- fread('GSE98895_non-normalized.txt.gz')
a <- a[,1:82]
colnames(a) <- c("PROBE_ID", "Gene", paste(rep(rownames(pdata), each = 2),
    rep(c("AVG_Signal", "Detection Pval"), 40),
    sep = "."
))
fwrite(a, file = "tmp.txt", sep = "\t", quote = F)

探针信息

probe2symbol <- data.frame(ProbeID=a$PROBE_ID, Feature=a$Gene)

lumi处理原始数据并提取标准化后的表达矩阵

conflicts_prefer(data.table::second)
conflicts_prefer(dplyr::first)
x.lumi <- lumiR("tmp.txt", lib.mapping = "lumiHumanIDMapping")
lumi.N.Q <- lumiExpresso(x.lumi)
dat <- exprs(lumi.N.Q)

ID转换

第一次ID转换，把表达矩阵中的第一列ID转换成第二列PROBE_ID

probeid <- lumi.N.Q@featureData@data
probeid <- probeid[match(rownames(dat), rownames(probeid)),,drop=FALSE ]
identical(rownames(dat), rownames(probeid))
probeid$PROBE_ID -> rownames(dat)

	PROBE_ID
Ku8QhfS0n_hIOABXuE	ILMN_1762337
fqPEquJRRlSVSfL.8A	ILMN_2055271
ckiehnugOno9d7vf1Q	ILMN_1736007
x57Vw5B5Fbt5JUnQkI	ILMN_2383229
ritxUH.kuHlYqjozpE	ILMN_1806310
QpE5UiUgmJOJEkPXpc	ILMN_1779670

第二次ID转换，把ProbeID转换成基因名字，把表达矩阵中的探针名转换为基因名；transid是我写的一个R函数，有需要可以联系我，加入交流群

probe_exprs <- as.data.table(dat, keep.rownames = "ProbeID")
fdata <- transid(probe2symbol, probe_exprs)

ProbeID	Feature
ILMN_1762337	7A5
ILMN_2055271	A1BG
ILMN_1736007	A1BG
ILMN_2383229	A1CF
ILMN_1806310	A1CF
ILMN_1779670	A1CF

保存数据

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

本系列GEO处理数据过程文件包括代码、原始数据文件、中间数据文件、结果文件之后都会更新在以下百度网盘链接中

https://mp.weixin.qq.com/s/bbm49F8hC8kwIWdZM1lbnw

本文由mdnice多平台发布