零底子入门转录组分析——数据处理（GEO数据库——高通量测序数据）
零底子入门转录组分析——数据处理（GEO数据库——高通量测序数据）GEO数据库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技术信息中心NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据，也就是说只要是目前已经发表的论文，论文中涉及到的基因表达检测的数据都可以通过这个数据库中找到。
而且GEO网站这个网站作为各种高通量实验数据的公共存储库。这些数据包括基于单通道和双通道微阵列的实验，检测mRNA，基因组DNA和蛋白质丰度，以及非阵列技术，如基因表达系列分析（SAGE），质谱蛋白质组学数据和高通量测序数据。可以按照文献数据集编号等浩繁形式进行检索。但是在这篇教程中仅先容如何从GEO网站上根据数据集编号下载所必要的GEO数据集，而且下载后在R中对数据集进一步处理成后续分析所要的形式。
本项目以妊娠期糖尿病GSE154414数据集(高通量测序数据)作为展示
选用的数据库是GEO。
实验分组：疾病组，对照组。
我这里使用的R版本是4.2.2
要用到的R包：tidyverse，GEOquery
1. 数据集获取
起首进入GEO网站官网（如下图所示），在检索位置输入数据集编号，点击箭头指向的位置进一步运行搜索。
https://img-blog.csdnimg.cn/direct/2d0a258a50d34a68957ebc96bc10a463.png
搜索之后会弹出如下界面：起首必要查抄物种范例（Homo sapiens），之后查看数据集的范例是否是高通量测序/芯片数据，我这里是高通量测序数据（Expression profiling by high throughput sequencing），页面往下拉。
https://img-blog.csdnimg.cn/direct/fb622c63c0b14a5a988ffa1617557776.png
如下图所示：包含了该数据集对应的解释文件GPL20301，而且还列出来了数据会合包含的样本。
注：但是解释信息可以不消过多的关注，由于后续分析用不到，样本数目可以大致瞅一眼
https://img-blog.csdnimg.cn/direct/44855e4fc99d4e87ac28f3cda9534c68.png
到此对于该数据集已经有了初步了解（实际上就是看是不是高通量测序数据），如果是高通量测序数据就按照下面的操纵进行，如果是芯片数据，可以参考之前的教程零底子入门转录组分析——数据处理（GEO数据库——芯片数据）
2. 数据处理（Rstudio）
rm(list = ls()) # 删除工作空间中所有的对象
setwd('/XX/XX/XX') # 设置工作路径
if(!dir.exists('./00_rawdata')){
dir.create('./00_rawdata')
} # 判断该工作路径下是否存在名为00_rawdata的文件夹，如果不存在则创建，如果存在则pass
setwd('./00_rawdata/') # 设置路径到刚才新建的00_rawdata下
加载包：
library(GEOquery)
library(tidyverse)
标注一下必要下载的数据集编号，而且在当前00_rawdata文件夹下创建一个名为GSE154414的文件夹（这是为了方便管理，如果不单独创建文件夹，数据集许多的话，就会显得很乱）
GEO_data <- 'GSE154414'
if(!dir.exists(paste0(GEO_data))){
dir.create(paste0(GEO_data))
}
setwd(paste0(GEO_data))
获取数据集：
[*]第一行指令是ncbi存放数据的地址
[*]第二行是拼接网址为GSE154414原始count存放地址（细心的小同伴或许已经观察到里面acc=XXX，file=XXX这里就是我们的GSE154414，这里实在用到了一点爬虫的思路！）
[*]第三行指令是获取GSE154414的原始count数据，等候下载完成
# load raw count of data_set
urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?format=file&type=rnaseq_counts"
path <- paste(urld, "acc=GSE154414", "file=GSE154414_raw_counts_GRCh38.p13_NCBI.tsv.gz", sep="&");
tbl <- as.matrix(data.table::fread(path, header=T, colClasses="integer"), rownames="GeneID")
https://img-blog.csdnimg.cn/direct/93caa8138d994c6bbf8e86b9ab6b97b1.png
下载完成后，优先转换成数据框格式并保存（这是由于有的数据集非常大，下载要很久，有时间碰到网欠好的时间还容易下载报错，各种失败，以是要及时保存）
expr <- as.data.frame(tbl)
write.csv(expr, paste0(GEO_data, '_raw_count.csv'))
可以先看一下下载好的这个名为expr 的数据长啥样，如下图所示，活动基因的ID，列为样本名，表中的那些整数就是raw_count。
https://img-blog.csdnimg.cn/direct/5c82ae0eff2e40b6a9ca0681223af1c6.png
接下来就是将这个表达矩阵的行名转换成我们寻常看到的基因symbol，要否则谁知道这些是个什么东西(>_>)
下载基因解释文件（这一步是用官方的解释文件）同样也是用的上述相似的代码获取，如下所示：
# load gene annotations
apath <- paste(urld, "type=rnaseq_counts", "file=Human.GRCh38.p13.annot.tsv.gz", sep="&")
annot <- data.table::fread(apath, header=T, quote="", stringsAsFactors=F, data.table=F)
rownames(annot) <- annot$GeneID
write.csv(annot, 'gene_annoation.csv')
同样打开下载好的这个名为annot的基因解释文件，看看都是些什么，如下图所示：
[*]第一列就是与表达矩阵行名一样的基因ID
[*]第二列就是基因的symbol（也就是咱们寻常见的基因名）
[*]第三列就是关于基因的简要描述
[*]第五列是基因范例（有的基因是编码基因，有的是非编码基因）
[*]反面还有基因对应的EnsemblGeneID，基因长度，位于染色体的位置，以及是否是正链/负链等信息（可以按需从这个解释表中获取基因的相应信息）
https://img-blog.csdnimg.cn/direct/79a9c0cfa4df419aa22c6fc3340fbb37.png
留意：这个解释文件适用于基本全部GEO高通量数据（因此只需下载一次即可，做好保存，碰到其他GEO高通量数据可以直接调用，不消重新下载）
接下来将这个基因解释表调整成后续必要用的格式（我这里第一步筛选了一下编码基因，如果不必要筛选的小同伴可以把这一步解释掉）
table(annot$GeneType) # 查看GeneType都有哪些类型
annot <- annot # 选择编码基因
colnames(annot) # 查看列名
probe2symbol <- dplyr::select(annot, 'GeneID', 'Symbol') # 选择基因ID和基因symbol这两列
probe2symbol <- filter(probe2symbol, `Symbol` != '') # 去掉Symbol中为空的行
names(probe2symbol) <- c('ID', 'symbol') # 给GeneID和Symbol两列名重命名为ID和symbol
probe2symbol <- na.omit(probe2symbol) # 去除为NA的行
终极处理好的probe2symbol 如下图所示，第一列是基因ID，第二列为基因symbol
https://img-blog.csdnimg.cn/direct/b48f3b4f314546c68dd46d282df09f41.png
接下来要将前面的expr表达矩阵和处理好的probe2symbol 联系到一起（说白了就是给expr表达矩阵的行名解释成我们寻常看到的基因名，而且加个去重处理）
expr$ID <- rownames(expr) # 给expr定义一个名为ID的新列，每一行的值为expr的行名
dat <- expr # 将expr赋值给dat
dat$ID <- as.character(dat$ID) # 将dat的ID列转成字符串型
probe2symbol$ID <- as.character(probe2symbol$ID)
dat <- dat %>%
merge(probe2symbol, by='ID')%>%
dplyr::select(-ID)%>% ## 去除多余信息
dplyr::select(symbol, everything())%>% ## 重新排列
mutate(rowMean = rowMeans(.))%>% ## 求出平均数
arrange(desc(rowMean))%>% ## 把表达量的平均值从大到小排序
distinct(symbol, .keep_all = T)%>% ## symbol留下第一个
dplyr::select(-rowMean)%>% ## 反向选择去除rowMean这一列
tibble::column_to_rownames(colnames(.)) ## 把第一列变成行名并删除
终极处理好的dat 如下图所示，列名为样本名，行名为平经常见的基因名（可以看出来这个和没处理前expr的最大区别就是：基因ID转成了熟悉的基因名）。
https://img-blog.csdnimg.cn/direct/677d5bbdf6e649c48b205af54a444dcf.png
接下来获取样本的临床信息表，其中getGEO函数是从GEO数据库中获取数据集的相关信息（包括样本的临床信息和表达矩阵等，但是由于这个数据集是高通量测序数据，因此表达矩阵不在这里）。
gset <- getGEO(GEO_data ,
destdir = '.',
GSEMatrix = T,
getGPL = F)
如下图所示gset是一个列表，phenoData里面就含有样本的临床信息。
https://img-blog.csdnimg.cn/direct/865ea24c2cdb4aa9a91cafd49d568491.png
通过pData函数获取样本临床信息
a <- gset[]
pd <- pData(a)
pd文件如下图所示每一行是一个样本，每一列为样本的不同临床信息，包括但不局限于：样本名，肿瘤分期，肿瘤分级，生存信息，样本范例等。
https://img-blog.csdnimg.cn/direct/14e7c18253a54534acadd488b16e7dc5.png
注：这个样本信息表中每一列命名不是规范的，如果想要获取本身必要的哪列临床信息（例如：样本范例——是否是癌，生存信息等），必要一一去查看对应的列中信息都是什么
在这里我们必要的是样本范例信息（用于做分组），这个信息存储在characteristics_ch1.2这列中
table(pd$characteristics_ch1.2)
table之后可以看到这列中包含两种范例，每种范例有4个样本：
[*]patient diagnosis: gestational diabetes mellitus（患者诊断:妊娠期糖尿病）
[*]patient diagnosis: healthy control（患者诊断:康健对照）
https://img-blog.csdnimg.cn/direct/6d1710023ac146b0814a5e9ad9f52831.png
接下来就可以根据这列信息做成我们本身想要的分组形式（即：一列是样本ID，另一列是分组环境），好比：将patient diagnosis: gestational diabetes mellitus作为疾病组，patient diagnosis: healthy control作为对照组
group <- data.frame(sample = pd$geo_accession, group = pd$characteristics_ch1.2) # 创建一个数据框，包含两列：sample和group
table(group$group)
group$group <- ifelse(group$group == 'patient diagnosis: healthy control', 'control', 'disease') # 判断，如果group列中值为patient diagnosis: healthy control，则重命名为control，否则为disease
table(group$group)
group <- group
group文件如下图所示：第一列是样本ID，第二列是分组范例。
https://img-blog.csdnimg.cn/direct/b0bce2630957475195a2ade7174d33a5.png
到这一步基本表达矩阵和样本分组信息表已经处理好了，通过 XX %in% XX筛选样本（这一步是由于有的数据会合不光有癌与癌旁，还有大概会有别的范例的样本，在前面做完筛选后和处理好的表达矩阵样本数目就不一样了，因此要保持一致）后就可以保存了。
colnames(dat)
group$sample
dat <- dat[, colnames(dat) %in% group$sample] # 筛选group中存在的样本
dat <- na.omit(dat)
group <- group # 筛选dat中存在的样本
write.csv(dat, file = paste0('dat.', GEO_data, '_count.csv'))
write.csv(group, file = paste0('group.', GEO_data, '.csv'), row.names = F)
到这一步数据集原始表达矩阵（raw_count）就已经处理完了，这个表达矩阵可以接下来通过DESeq2进行差异分析。
注：raw_count仅适用于做差异分析，如果是做其他分析，则必要对raw_count进行标准化处理，转成FPKM后才气用于其他范例的分析（例如：WGCNA，生存分析等）
接下来就是对raw_count做标准化处理，将其转换成FPKM形式。
3. 数据标准化（Rstudio）
关于FPKM的解释，可以参考RPKM、FPKM 和 TPM，解释清楚这篇文章（简单说就是由于存在测序深度不同/基因长度不同的环境，以是原始count并不能真实反映基因表达量，有的基因长度较长，匹配到的reads会非常多，因此会根据基因长度做个校正）
起首要从前面annot文件中获取基因的长度
colnames(annot)
gene_length <- dplyr::select(annot, Symbol, Length)
定义一个名为countToFpkm 的计算函数：
countToFpkm <- function(counts, effLen){
N <- sum(counts)
exp( log(counts) + log(1e9) - log(effLen) - log(N) )
}
通过apply函数和之前自定义的countToFpkm 函数对表达矩阵进行FPKM转换
dat <- dat
expr_fpkm <- apply(dat, 2, countToFpkm, effLen = gene_length$Length)
expr_fpkm就是经过FPKM转换的数据集（但是目前还不能直接用，是由于除了差异分析之外，基本全部的分析点用到的数据集都是log2（）处理的，这是由于纯fpkm值有的会很大，不同分析点在计算的时间大概会由于个别值特别大会导致较大的误差，因此要通过log2处理。）
接下来通过GEO官方的判定公式，如果logC的值为TRUE则对expr_fpkm数据集进行log2处理，如果为FALSE则不处理
qx <- as.numeric(quantile(expr_fpkm, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx > 100) ||
(qx-qx > 50 && qx > 0) ||
(qx > 0 && qx < 1 && qx > 1 && qx < 2)
if (LogC) {
expr_fpkm <- NaN
expr_fpkm <- log2(expr_fpkm)
print("log2 transform finished")
}else{
print("log2 transform not needed")
}
注：这里必要特别留意一点的是这个判定函数并不是万能的，有些数据集在这一步判定的时间为FALSE，但是里面的fpkm值都很大，像这种特别环境的就不能用判定了，可以直接对其做log2处理
末了将处理好的数据集保存即可
expr_fpkm <- data.frame(expr_fpkm)
write.csv(expr_fpkm, file = paste0('dat.', GEO_data, '_fpkm.csv'))
补充内容：
[*]为什么要进行log2处理：一样寻常数据集的表达值范围都是非常大的，从0到25000不等。将这些数据绘制密度分布后，一样寻常出现右偏，即大部分信号都是在左侧，右侧拖个长长的尾巴（偏态分布），倒霉于研究，而经过log2转化后，数据更加会合，更加接近正态分布，更方便套用正态分布那一套理论进行后续分析。
以上就是GEO高通量测序数据下载及处理的全部过程，如果有什么必要补充或不懂的地方，大家可以私聊我大概在下方评论。
如果以为本教程对你有所帮助，点赞关注不迷路！！！
与教程配套的原始数据+代码+处理好的数据见零底子入门转录组分析——数据处理（GEO数据库——高通量测序数据）配套资源
注：配套资源只要改个路径就能运行，本人已检测过可以跑通，请放心食用，食用过程碰到标题，可先自行百度，实在解决不了可以私信
[*]目录部分跳转链接：零底子入门生信数据分析——导读
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。
页:
[1]