三阴性乳腺癌表达矩阵探索之数据下载及理解

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！（视频观看方式见文末）

下面是《GEO数据挖掘课程》的配套笔记（第二篇）

了解数据挖掘

• 公共数据库：（数据来源）
• GEO和TCGA
• 国际三大数据中心：NCBI, ENSEMBL, UCSC
• 数据挖掘的概念
• 从大的数据背景中通过各种统计学方法得到数量大小合适的基因集找到的感兴趣的基因集
• 通过各种统计学方法来注释并解释这个基因集的意义

实战：

对文献解读的第三篇文章==Identification of Key Genes and Pathways in Triple-Negative Breast Cancer by Integrated Bioinformatics Analysis== 的分析过程进行重复

• 第一步：下载数据集

  GEO数据库基本介绍：

  一篇文章可以有一个或者多个GSE数据集，一个GSE里面可以有一个或者多个 GSM样本 ,多个研究的GSM样本介意根据研究目的整合为一个 GDS , 不过GDS本身用得很少，而且每个数据集都有自己对应的芯片平台，就是GPL

  GEO Platform：GPL

  GEO Sample: GSM

  GEO Series: GSE

  GEO Dataset: GDS

  GEO数据库，根据数据存放的标签GSE号进行查询

  

  找到GSE号

  进如GEO数据库

  

  进入GEO并搜索数据集

  点击目标查询进入目标数据集网页

  

  下载数据的详细介绍

  

  探针注释平台的位置

  

  表达矩阵下载位置

  表达矩阵下载的方式：

  

  表达矩阵下载的方法一

  

  表达矩阵下载方式二

1. 使用GEOqueryR 程序包从GEO数据库下载

   ==Note==:使用下面的代码下载的文件都会保存到本地，destdir参数指定数据存放的位置。此外，比较重要的三个参数为GSEMatrix=TRUE,AnnotGPL=FALSE, getGPL=TRUE

   #加载程序包
library(GEOquery)

   #根据GDS下载soft文件
   gds <- getGEO('GD858', destdir='.')

   #根据GPL号下载芯片设计信息
   gpl96 <- getGEO("GPL96", destdir=".")

   #根据GSE号下载series_matrix.txt.gz
gse1009 <- getGEO("GSE1009",dstdir=".")


1. 下载原始芯片表达数据（CEL）
1. 直接下载matrix文件，点击'Series Matrix File(s)’进入到矩阵存放位置，直接点击下载

• 第二步：开始分析

1. 新建一个R.project GSE76275.Rproj

2. 在新的project下分别创建每个流程的分析

3. 总共分step0-step5

   step0-install.R : 安装需要用到的程序包

   Notes: R版本高于3.5 使用BiocManager， 低于3.5用BiocInsraller

   rm(list = ls()) 
#清空当前工作空间变量  
options()$repos  
#查看当前工作空间默认的下载包路径
options()$BioC_mirror 
#查看使用BioCManager下载包的默认路径
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") 
# 指定使用BioCManager下载的路径
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) 
# 指定使用install.packages下载包的路径
options()$repos 
options()$BioC_mirror

   # https://bioconductor.org/packages/release/bioc/html/GEOquery.html
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager") 
#判断是否存在BiocManger包，没有的话下载该包 
#判断是否存在这些包，不存在的话安装这些包
if(!require("KEGG.db")) BiocManager::install("KEGG.db",ask = F,update = F)
if(!require("GSEABase")) BiocManager::install("GSEABase",ask = F,update = F)
if(!require("GSVA")) BiocManager::install("GSVA",ask = F,update = F)
if(!require("clusterProfiler")) BiocManager::install("clusterProfiler",ask = F,update = F)
if(!require("GEOquery")) BiocManager::install("GEOquery",ask = F,update = F)
if(!require("limma")) BiocManager::install("limma",ask = F,update = F)
if(!require("impute")) BiocManager::install("impute",ask = F,update = F)
if(!require("genefu")) BiocManager::install("genefu",ask = F,update = F)
if(!require("org.Hs.eg.db")) BiocManager::install("org.Hs.eg.db",ask = F,update = F)
if(!require("hgu133plus2.db")) BiocManager::install("hgu133plus2.db",ask = F,update = F)
if(!require("ConsensusClusterPlus")) BiocManager::install("ConsensusClusterPlus",ask = F,update = F)


   step1-download.R: 下载所需要的数据

   ##1.获取GEO数据
library(GEOquery)
f = "SE76275_eSet.Rdata" #如果文件不存在则进行下载
if(!
   file.exist(f)){ 
    gset <- getGEO("GSE76275", destdir=".", 
               AnnotGPL=T, #注释文件，可下可不下
               getGPL = T) #注释平台，可下可不下，可以改为F
    save(gset,file="GSE76275_eSet.Rdata") #保存到本地
     }
load("GSE76275_eSet.Rdata") #载入数据


   简单对下载的数据进行了解：

   ExpressionSet数据形式的组成：

   assayData

   phenoData

   featureData

   experimentData

   protocalData

   class(gset) #list类型
length(gset) #查看长度，只有一个元素
class(gset[[1]]) #取出第一个元素，并查看类型为"ExpressionSet"
?ExpressionSet #查看这个数据类型，getGEO函数的目的就是下载数据，而下载的数据最终以ExpressionSet的形式存在

   a<-gset[[1]] #取出该列表的第一个元素并赋值
   a@experimentData #访问不同的数据集
   a@assayData

   methods(class='ExpressionSet') #可用于查看该对象的操作函数
dat=exprs(a) #取出对象a中的表达矩阵
dim(dat) #检查维度，54675个探针，265个病人
dat[1:4,1:4] #查看前四行前四列


   了解实验设计：

   pd <- pData(a) #取出
pd$characteristics_ch1.1 #取出分组信息
ifelse(X==1,'X等于1'，'X不等于1') #首先判断X是否等于1，如果X等于1，返回'X等于1'的值，否则返回'X不等于1'的值

   group_list<-ifelse(pd$characteristics_ch1.1=='triple-negative status: not TN','noTNBC','TNBC') #获得分组信息

   save(dat,group_list, file='step1-output.RData')


   ==总结==：第一步结束之后，我们得到了连个数：256个样本在5万个探针上的表达量，256个样本的分组信息，此处的特征为是否是三阴性乳腺啊

   step2-check.R: 检查程序包以及数据

   (1) 检查PCA分组情况

   rm(list=ls()) #清空之前的一些值，避免混淆
## 学习PCA包的示例
install.packages(c('FactoMineR','factoextra'))
library(FactoMineR)
library(factoextra)
head(dat)
dat <- t(dat)  #交换行和列的位置
dat <- as.data.frame(dat)
dat <- cbind(dat,group_list)
#在进行PCA分析之前，变量group_list(index=54676)被移除
dat.pca <- PCA(dat[,-54676],graph = FALSE)
fviz_pca_ind(dat.pca,
             geom.ind = "point",#仅仅展示点
             col.ind=dat$group_list, #通过组群来进行染色
             polette = c("red","green"), #调色盘，几个组就设置几个颜色
             addEllipses = TRUE, #以椭圆方式集中
             legend.title="Group")

   

   PCA

   (2)检查热图,选取其中的一小部分基因绘制热图

   ##绘制热图
rm(list=ls())
load(file='step1-output.RData')
dat[1:4,1:4]
#为了更好的用热图体现基因的表达，需要挑选表达差异比较大的基因
#apply函数介绍
apply(mat,1,mean) #对矩阵mat的每一行求均值，第二个参数为1对行进行操作，为2则对列进行操作
apply(dat,1,sd)
system.time(apply(dat, 1, sd))

   cg <- names(tail(sort(apply(dat,1,sd)),1000)) #挑选sd值比较大的基因
   library(pheatmap)
   n <- t(scale(t(dat[cg,])))
   n[n>2]=2
   n[n<-2]=-2
   pheatmap(n, show_rownames = F,show_colnames = F) #scale,对样本进行均一化
   ac <-data.frame(g=group_list) #添加分组信息
   rownames(ac) = colnames(n)
   pheatmap(n, show_rownames = F,show_colnames = F, annotation_col = ac)

   

   heatmap_top1000_sd

   后面将详细介绍以下三个分析

   step3-DEG.R: 表达差异性分析

   step4-Annotation.R: 找出差异基因之后，看下这些差异基因属于哪些通路

   step5-wgcna.R

视频观看方式

我把3年前的收费视频课程：3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月，免费到B站:

• 这个课程超级棒，B站免费学习咯：https://m.bilibili.com/video/BV1dy4y1C7jz
• 配套代码在GitHub哈：https://github.com/jmzeng1314/GSE76275-TNBC
• TCGA数据库挖掘，代码在：https://github.com/jmzeng1314/TCGA_BRCA
• GTEx数据库挖掘，代码在：https://github.com/jmzeng1314/gtex_BRCA
• METABRIC数据库挖掘，代码在：https://github.com/jmzeng1314/METABRIC

然后马上就有了3千多学习量，而且有学员给出来了图文并茂版本万字笔记，让我非常感动！

扫描下面二维码马上就可以学习起来啦，笔记需要至少半个小时来阅读哦！