零基础入门转录组下游分析——数据处理惩罚（GEO数据库——芯片数据） ...

零基础入门转录组下游分析——数据处理惩罚（GEO数据库——芯片数据）


  

---

---

GEO数据库全称GENE EXPRESSION OMNIBUS，是由美国国立生物技能信息中央NCBI创建并维护的基因表达数据库。它创建于2000年，收录了世界各国研究机构提交的高通量基因表达数据，也就是说只要是目前已经发表的论文，论文中涉及到的基因表达检测的数据都可以通过这个数据库中找到。
并且GEO网站这个网站作为各种高通量实行数据的公共存储库。这些数据包括基于单通道和双通道微阵列的实行，检测mRNA，基因组DNA和卵白质丰度，以及非阵列技能，如基因表达系列分析（SAGE），质谱卵白质组学数据和高通量测序数据。可以按照文献数据集编号等众多情势举行检索。但是在这篇教程中仅先容如何从GEO网站上根据数据集编号下载所必要的GEO数据集，并且下载后在R中对数据集进一步处理惩罚成后续分析所要的情势。

---

---

1. 本项目以非小细胞肺癌GSE37745数据集（肺腺癌）作为展示
3. 选用的数据库是GEO。
5. 实验分组：疾病组，对照组。
7. 我这里使用的R版本是4.2.2
9. 要用到的R包：tidyverse，GEOquery

复制代码

1. 数据集获取

起首进入GEO网站官网（如下图所示），在检索位置输入数据集编号，点击箭头指向的位置进一步运行搜索。

搜索之后会弹出如下界面：起首必要查抄物种类型（Homo sapiens），之后查察数据集的类型是否是高通量测序/芯片数据，我这里是芯片数据（Expression profiling by array），页面往下拉。

如下图所示：包含了使用该数据集的一些文献，以及该数据集对应的注释文件，并且还列出来了数据集中包含的样本。在这里我们比力关注的是注释文件GPL570。

可以点击GPL570查察注释文件信息，页面拉到最后，如下图所示，这个注释文件里包含了芯片的ID及其对应的基因symbol，这也是我们后面必要用到的内容。（并且我们可以看到基因symbol中有的基因后面有三条反斜杠，那是他的别名，后续必要删除）

注意：GEO数据集如果是芯片数据则不必要手动下载，如果是高通量测序的数据则必要手动下载表达矩阵。
2. 数据处理惩罚（Rstudio）

获取以上信息后即可直接进入R用代码自动下载

1. rm(list = ls()) # 删除工作空间中所有的对象
2. setwd('/XX/XX/XX') # 设置工作路径
3. if(!dir.exists('./00_rawdata')){
4.   dir.create('./00_rawdata')
5. } # 判断该工作路径下是否存在名为00_rawdata的文件夹，如果不存在则创建，如果存在则pass
6. setwd('./00_rawdata/') # 设置路径到刚才新建的00_rawdata下

复制代码

加载包：

1. library(GEOquery)
2. library(tidyverse)

复制代码

标注一下必要下载的数据集编号，及前面查察的注释文件编号，并且在当前00_rawdata文件夹下创建一个名为GSE37745的文件夹（这是为了方便管理，如果不但独创建文件夹，数据集很多的话，就会显得很乱）

1. GEO_data <- 'GSE37745'
2. gene_annotation <- 'GPL570'
4. if(!dir.exists(paste0(GEO_data))){
5.   dir.create(paste0(GEO_data))
6. }
7. setwd(paste0(GEO_data))

复制代码

获取数据集：

1. gset <- getGEO(GEO_data , # 前面创建GEO_data对象
2.                destdir = '.',
3.                GSEMatrix = T,
4.                getGPL = F) # getGEO函数自动从官网中获取对应的数据集
5. expr <- as.data.frame(exprs(gset[[1]])) # 将获取的数据集的表达矩阵提取出来，如果没有表达矩阵，说明是高通量数据需要到网站自己下载

复制代码

起首来看一下获取到的表达矩阵长啥样：行名为芯片的ID（后面必要比对注释文件改成对应的基因symbol），列名为样本ID。

获取注释文件GPL570：

1. gpl <- getGEO(gene_annotation, destdir = '.')
2. gpl <- Table(gpl)
3. colnames(gpl)
4. [1] "ID"                               "GB_ACC"                          
5. [3] "SPOT_ID"                          "Species Scientific Name"         
6. [5] "Annotation Date"                  "Sequence Type"                  
7. [7] "Sequence Source"                  "Target Description"              
8. [9] "Representative Public ID"         "Gene Title"                     
9. [11] "Gene Symbol"                      "ENTREZ_GENE_ID"                  
10. [13] "RefSeq Transcript ID"             "Gene Ontology Biological Process"
11. [15] "Gene Ontology Cellular Component" "Gene Ontology Molecular Function"

复制代码

这里面我们必要用到的就是ID和Gene Symbol，进一步处理惩罚ID和基因ID，前面我们也说了注释文件的基因symbol里有的有别称，必要剔除，用到的就是separate函数，可以按照sep = 'XXX’的情势分割。

1. probe2symbol <- dplyr::select(gpl, 'ID', 'Gene Symbol')
2. probe2symbol <- filter(probe2symbol, `Gene Symbol` != '')
3. probe2symbol <- separate(probe2symbol, `Gene Symbol`, into = c('symbol', 'drop'), sep = '///')
4. probe2symbol <- dplyr::select(probe2symbol, -drop)
5. names(probe2symbol) <- c('ID', 'symbol')

复制代码

可以看一下处理惩罚后的probe2symbol长啥样，如下图所示，第一列就是芯片ID，第二列就是基因symbol，做到这大家也就可以想象到后面的步调了，无非就是做个匹配，将表达矩阵中的ID对应成相应的基因symbol。

如下代码所示进一步处理惩罚表达矩阵：

1. dat <- expr
2. dat$ID <- rownames(dat)
3. dat$ID <- as.character(dat$ID)
4. probe2symbol$ID <- as.character(probe2symbol$ID)
6. dat <- dat %>%
7.   merge(probe2symbol, by='ID')%>% ## 合并注释文件和表达矩阵
8.   dplyr::select(-ID)%>%     ## 去除多余信息——芯片ID
9.   dplyr::select(symbol, everything())%>%     ## 重新排列
10.   mutate(rowMean = rowMeans(.[grep('GSM', names(.))]))%>%    ## 求出平均数
11.   arrange(desc(rowMean))%>%       ## 把表达量的平均值从大到小排序
12.   distinct(symbol, .keep_all = T)%>%      ## symbol留下第一个
13.   dplyr::select(-rowMean)%>%     ## 反向选择去除rowMean这一列
14.   tibble::column_to_rownames(colnames(.)[1])   ## 把第一列变成行名并删除

复制代码

处理惩罚后的dat文件，如下图所示：行名为基因symbol，列为样本名。到了这一步表达矩阵已经处理惩罚好了，接下来就是新的问题：样本是怎么分组的？

前面的getGEO函数不仅能获取数据集的表达矩阵，还可以获取杨信息，只必要通过pData函数即可展示样本信息

1. a <- gset[[1]]
2. pd <- pData(a) #

复制代码

如下图所示为pd对象的样本信息：有样本名，类型等等（可自行查察）

这里我们必要用到的是geo_accession和characteristics_ch1.2这两列（前一列就是样本ID，后一列就是分组信息），固然该数据集中有三种类型的疾病，分别是adeno（肺腺癌），large（大细胞），squamous（肺鳞癌）

1. table(pd$characteristics_ch1.2)
2. # histology: adeno    histology: large histology: squamous
3. # 106                  24                  66
4. group <- data.frame(sample = pd$geo_accession, group = pd$characteristics_ch1.2)
5. table(group$group)
6. group <- group[group$group != 'histology: large', ] # 筛选掉large（大细胞）
7. table(group$group)
8. group$group <- ifelse(group$group == 'histology: adeno', 'LUAD', 'LUSC')
9. table(group$group)
10. # LUAD LUSC
11. # 106   66

复制代码

做到这一步基本可以说是大功告成了，目前我们有了基因表达矩阵，尚有了样天职组信息，接下来将我们的分组样本信息和表达矩阵匹配一下，生存文件即可。

1. dat <- dat[, group$sample]
2. dat <- na.omit(dat)
3. write.csv(dat, file = paste0('dat.', GEO_data, '.csv'))
4. write.csv(group, file = paste0('group.', GEO_data, '.csv'), row.names = F)

复制代码

注：做到这一步，可以说是生信分析已经完成一半了，后续很多的分析都会基于前面处理惩罚的这两个文件，因此这最初的第一步也是最重要的一步，只有基础打好，后面分析才会可靠一些。

---

---

结语：
以上就是零基础入门转录组下游分析——数据处理惩罚（GEO数据库——芯片数据）的全部过程，如果有什么必要增补或不懂的地方，大家可以私聊我大概在下方批评。
如果觉得本教程对你有所资助，盼望广大学习者能够点赞，收藏，加关注
关于我们： 我们的团队是领航生信，如果大家想要系统学习常规SCI生信套路和流程，可以在微信中搜索领航生信，即可接洽到我们，感谢大家的支持
祝大家能够开心学习，轻松学习，在学习的路上少一些坎坷~~~

---

• 目录部分跳转链接：零基础入弟子信转录组数据分析——导读
  

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。