用R分析芯片实战（上）

刘小泽写于18.9.10
实战上从数据下载到差异基因的获得、初步作图
实战下进行富集分析，使用数据库进行注释
在此感谢jimmy的limma包等教程

1 下载芯片数据

使用Fabbrini E于2015年发表在J Clin Invest上的文章：Metabolically normal obese people are protected from adverse effects following weight gain

文章研究了体重适度增加对代谢正常(MNO)和代谢异常(MAO)受试者脂肪组织基因表达的影响，选取11个代谢正常的人、7个代谢异常的人增重前后的皮下脂肪组织，得到36个数据

文章下载链接：https://pdfs.semanticscholar.org/6052/f035dfb296559e19314352506f048e13f02e.pdf

使用的数据集是GSE62832

实验平台是GPL6244， Affymetrix Human Gene 1.0 ST Array

• 方法一：打开网页https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62832 【可能会有网速限制】
  GEO数据GSE62832

• 方法二：使用GEOquery包

  if(T){
  source("http://bioconductor.org/biocLite.R")
  options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
  biocLite("GEOquery")
  
  library(GEOquery)
  eSet <- getGEO('GSE62832', destdir = '.', getGPL = F, AnnotGPL = F)#一般芯片的注释文件比较大，直接下载不容易成功，后两个选项可以避免下载GPL、Annotation
  save(eSet, file = 'GSE62832.eSet.Rdata')
  }
  

下载完，看一下eSet这个对象都包含什么

> eSet
$GSE62832_series_matrix.txt.gz
ExpressionSet (storageMode: lockedEnvironment)
assayData: 33297 features, 36 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: GSM1534034 GSM1534035 ... GSM1534069 (36 total) #36个样本
  varLabels: title geo_accession ... tissue:ch1 (35 total)
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: GPL6244 #查看对应芯片平台

class(eSet) #先查看数据种类【列表还是数据框】
str(eSet) #再看数据结构【可以看到第二行涉及到了表达矩阵】
#既然是列表，那么从列表中提取信息就是 eSet[[1]]
e <- exprs(eSet[[1]])

2. 将表达矩阵的探针ID转换为gene ID

• 首先需要知道GSE62832对应平台是GPL6244

• 然后需要知道GPL6244对应哪个注释包【阅读jimmy之前整理的平台信息https://www.jianshu.com/p/f6906ba703a0】

• 下载相应的注释包

  source("https://bioconductor.org/biocLite.R")
  options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
  biocLite("hugene10sttranscriptcluster.db")
  

• 开始探索、过滤、整合

  library(hugene10sttranscriptcluster.db) 
  ls("package:hugene10sttranscriptcluster.db") #查看所有包含的对象
  s <- toTable(hugene10sttranscriptclusterSYMBOL) #将SYMBOL对象转为数据框
  #【补充：这个注释包是别人上传的，我们只是拿来用。其中的表达矩阵中原来有3w多探针，作者过滤后只留下接近2w的探针。这仅仅是作者过滤后的结果，并非原始数据结果】
  

  简单的探索

  library(magrittr) 
  s$symbol %>% unique() %>% length() #看一下有多少基因【人类正常蛋白编码基因就2w左右】
  s$symbol %>% table() %>% sort() %>% tail() #基因使用探针最多的前6名【发现有两个基因用了10个探针】
  s$symbol %>% table() %>% sort() %>% table()
  

  过滤+整合

  dim(e) #过滤前
  #过滤【只保留和注释文件探针id相同的探针】
  efilt <- e[rownames(e)%in%s$probe_id,]
  dim(efilt)#过滤后
  #整合1【目的：保证一个基因对应一个探针；如果基因和探针一一对应很好说，但如果一个基因对应多个探针：每个探针取一行的均值-》对应同一基因的探针取表达量最大的探针-》按照基因名给他们建索引，因为是按照基因来过滤探针（不用s$probe_id构建索引的原因是，看清楚我们的目的是让注释包的一个基因对应我们自己表达矩阵的一个探针。如果用s$probe_id那么结果就成了让注释包的一个探针对应我们自己表达矩阵的一个探针，当然这样也运行不成功，因为自己表达矩阵的探针过滤后的数量和注释包的探针数量不相等，这样没法一一对应。但基因名数量是不变的，什么是索引？以不变应万变的就是索引）】
  maxp = by(efilt,s$symbol,function(x) rownames(x)[which.max(rowMeans(x))]) 
  uniprobes = as.character(maxp)
  efilt=efilt[rownames(efilt)%in%uniprobes,]
  #整合2【目的：将我们表达矩阵的行名换成刚才一对一的基因名，并且match这个函数保证了表达矩阵和注释包的顺序是一致的】
  rownames(efilt)=s[match(rownames(efilt),s$probe_id),2]
  

3. 对表达矩阵进行一些检验

表达矩阵不局限于GEO数据库的芯片分析，转录组及其他涉及基因、样本、分组关系的都会有一个表达量矩阵，就是一个基因在不同样本中（对照、处理；是否患病等）的表达差异。
拿到表达矩阵，在进行后续分析之前，首先要检测这个矩阵是不是合理的，比如看管家基因是否表达量突出、一致；样本分组是不是和实验设计一致，用PCA、hclust检验

3.1 检测一些管家基因表达量

boxplot(efilt[,1])#看看第一个样本中总体基因表达量分布，可以看到基本为5左右
efilt['ACTB',] #激动蛋白Beta-actin的基因名是ACTB，管家基因
efilt['GAPDH',] #也是管家基因

3.2 看表达矩阵的整体分布

先把表达矩阵=》tidy data【四列：基因名、样本、表达量、表型分组(看文献按MAO、MNO分组)】

library(reshape2)
pdata=pData(eSet[[1]]) #将样本表型信息从数据框中提取出来【取出来的是表型、样本的数据框】
group_list=as.character(pdata$`metabolic status:ch1`) 
m_efilt = melt(efilt) #先将原来矩阵“融化
colnames(m_efilt)=c('symbol','sample','value') #重新命名三列
m_efilt$group=rep(group_list,each=nrow(efilt)) 

以图为证

检查数据

4. 对表达矩阵进行差异分析

只需要提供表达矩阵efilt、分组信息group_list，就能使用limma进行分析
详细内容参考jimmy的http://www.bio-info-trainee.com/bioconductor_China/software/limma.html

suppressMessages(library(limma))
#limma需要三个矩阵：表达矩阵（efilt）、分组矩阵（design）、比较矩阵（contrast）
#先做一个分组矩阵～design，说明MAO是哪几个样本，MNO又是哪几个，其中1代表“是”
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(efilt)
design
#再做一个比较矩阵【一般是case比control】
contrast<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast

4.1 准备就绪，就可以开始差异分析

DEG <- function(efilt,design,contrast){
  ##step1
  fit <- lmFit(efilt,design)
  ##step2
  fit2 <- contrasts.fit(fit, contrast) 
  fit2 <- eBayes(fit2)  
  ##step3
  mtx = topTable(fit2, coef=1, n=Inf)
  deg_mtx = na.omit(mtx) 
  return(deg_mtx)
}
DEG_mtx <- DEG(efilt,design,contrast) #得到全部的差异基因矩阵

4.2 对一个小表达矩阵，如30、50个基因，可以用热图

top30_gene=head(rownames(DEG_mtx),30)
top30_matrix=efilt[top30_gene,] #得到top30的表达量矩阵
top30_matrix=t(scale(t(top30_matrix)))
#这里做个top30 gene heatmap

前30基因的热图

4.3 对一个大的表达矩阵，如全部的差异基因，可以用火山图

火山图实际上就是根据两列进行作图：logFC、pvalue

plot(DEG_mtx$logFC, -log10(DEG_mtx$P.Value)) #先画一个最简单的图（能说明原理）
#再根据jimmy的代码做一个。网络上也有很多火山图的代码可以参考

火山图

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到Bioplanet520@outlook.com

Welcome to our bioinfoplanet!