GEO在学习-代码

2019-08-05 本文已影响40人小梦游仙境

GEO再学习

GEO数据挖掘已经成为生信学者必备技能,我以为自己会跑代码了就是会了,其实呢,细细去领会每个细节,还是会发现很多盲区.听了晓娟老师的课,发现很多很多细节我并不甚清,一定要明白,才能融会贯通.同时,晓娟老师又把思路捋地非常清,练习实操几次,一定可以像别人家的孩子那样,做到数据挖掘,复现文章美图.同时,如果自己有了数据,也可以这样分析了!

GEO数据挖掘思路

1.下载数据
2.数据预处理
3.差异分析
4.差异分析后的可视化-热图、火山图、富集分析

1.下载数据

(1)用GSE号下载

从文章中搜索GSE号,然后去GEO网址:https://www.ncbi.nlm.nih.gov/geo/下载,以``GSE42872`为例

Platform (GPL) :芯片或测序平台的相关信息
Sample (GSM) :单个样本的数据
Series (GSE) :一系列GSM整合后的数据
DataSets (GDS):由GEO的工作人员整合的样品数据

如果上面这些发懵,记不住,那看下图

用GEOquery包的下载方法

平台与对应的基因注释包

代码如下

library(GEOquery)
f<-'GSE42872.Rdata'
####getGPL获得平台的注释信息，但下载速度会慢很多
####而且注释文件格式大多不如bioconductor包好用
if(!file.exists(f)){
  gset<-getGEO('GSE42872',destdir='.',
               AnnotGPL=F,
               getGPL=F)
  save(gset,file=f)
}
#数据提取
load('GSE42872.Rdata')

(2)下载原始数据-txt.gz

某些情况下,下载的series Matrix File(s)文件并不是整理好的,这时可能需要下载原始文件,就是RAW.tar.参考学徒数据挖掘干扰MYC-WWP1通路重新激活PTEN的抑癌活性——3步搞定GSEA分析一样可以整理成想要的表达矩阵信息.文章中的GSE号为gse126789

下载RAW data

下载后的文件是txt.gz结尾的

image-20190804040439655

下载数据后,处理的代码如下

rm(list = ls())  
options(stringsAsFactors = F)
# 下载GSE126789_RAW.tar
# 批量读取文件处理为表达矩阵
dir <- "/Users/mengmeng/Downloads/GSE126789_RAW/"#上一步下载的文件路径

files <- list.files(path = dir, pattern = "*wwp1.*.txt.gz", recursive = T) #找到wwp1的WT和敲除样本,也就是文件名,用的list.files函数仅仅是文件的名字,并不包括里面的数据
expr <- lapply(files,
               function(x){
                 # 只读取基因名和count
                 expr <- read.table(file = file.path(dir,x), sep = "\t", header = T, stringsAsFactors = F)[,c(1,7)]#取出来第一列和第七列,即基因名和count计数,矩阵中取不相邻的两列就可以[,c(1,7)]这么取,也可以cbind(b$Geneid,b$Count)取,但啰嗦
                 return(expr)#注意此时的files仅仅是文件名字,并不是解压后文件含有的矩阵信息,而这时的files是character,是字符型,要清楚的是,lapply不仅可以作用于列表list,同样可以应用在向量上.向量包括数值、字符、逻辑型向量
               }) #这个循环非常棒,要理解function(x)这个里面的x是什么.files是什么,x就是什么,files是文件名,那么x也是文件名,所以才有file = file.path(dir,x),这个file.path(dir,x)就相当于read.table后面括号中的文件名
df <- do.call(cbind, expr)#得到的每一个矩阵横向叠加(cbind就是列相加)
rownames(df) <- df[,1]
df <- df[,seq(2,ncol(df),by=2)] #去掉重复读取的基因名#seq函数:从第二列开始,直到最后一列,by间隔为2来取
colnames(df) <- substr(files,1,10) #列名为样本名,用substr提取字符串,仅仅保留样本名如GSM3613325
df <- df[apply(df, 1, sum)!=0,] #去掉在所有样本中count都为0的基因,取出来不等于0的矩阵
grouplist <- c(rep("KO",3),rep("WT",3)) #记录下分组信息,新建一个grouplist向量,即'KO'重复3次,'WT'重复3次
expr <- df
save(expr,grouplist,file = "./expr_group.Rdata")

(3)下载原始数据-CEL.gz

同上面一样,如果下载的RAW.tar数据是CEL.gz,那就去看jimmy老师的文章吧,用affy包读取affymetix的基因表达芯片数据-CEL格式数据

2.数据预处理

本次对表达矩阵进行差异分析时,用的是limma包.通常使用limma进行差异分析处理时，需要经过log2后的矩阵作为表达矩阵输入。根据log2FC的定义，这个数字表示变化倍数经过log2后的一个值，比如log2FC=1,则变化为2倍；log2FC=2,则变化为4倍。这是常用的一种表述方法。在使用limma函数计算时，如果输入的矩阵没有经过log2处理，则会把FC当成log2FC输入，这或许是因为limma默认输入的是log2后的表达式。这里有必要提到log的一个运算，即，可见对于已经log2后的数据，计算log2FC = log2(A/B)只需要直接使用log2A-log2B。所以如果给出的是一个未经log2的数值，函数也会直接相减以得到log2FC，这就导致计算出来的差异表达高达几百甚至上千。在判断counts是否需要重新标准化以及是否需要log2时，可以根据数值大小粗略估计。如果表达丰度的数值在50以内，通常是经过log2转化的。如果数字在几百几千，则是未经转化的。因为2的几十次方已经非常巨大，如果2的几百次方，则不符合实际情况。参考

1.https://blog.csdn.net/tuanzide5233/article/details/88542805,

2.https://blog.csdn.net/tuanzide5233/article/details/88542558

3.http://www.bio-info-trainee.com/1209.html

用代码查看矩阵是否需要log

qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
  (qx[6]-qx[1] > 50 && qx[2] > 0) ||
  (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
LogC
ex <- log2(ex+1)

用箱线图看数据的分布

boxplot(ex,las=2,cex.axis=0.6,main='data check')

image-20190804071930092

用PCA看分组,即检验现在的表达矩阵中的样本信息所对应的分组信息是否有以下的情况:
- 是否有离群样本
- 实验组和对照组是否正确(有没有标反)
- 有没有批次效应

主成分分析（PCA)是一种降维分析，将多个维度的数据简化降维，突出主要成分，所以叫做主成分分析。PCA是一种数据降维技巧，它能将大量相关变量转化为一组很少的不相关变量，这些无关变量称为主成分。例如，使用PCA可将30个相关（很可能冗余）的环境变量转化为5个无关的成分变量，并且尽可能地保留原始数据集的信息。

用一个通俗的例子来举例:

image-20190804071242809

这个是大英帝国四个地区的17种食物的消费数据（单位克/人/周），你能看出来四个地方的消费有何不同吗？或者说，四个地方哪几个更为接近？谁更为特殊？

如果光肉眼看，肯定是很费劲的，而且也不一定准？

从数据的结构上来看，以英国4个地区为自变量的话，这么这个数据实际上含有17个维度，每一个维度都不是完全一样，但是每一个维度的近似程度或者疏远程度是有差别的，也就是说，每一个维度对整体变异度的贡献不一样。

通俗点说，每种食物的消耗量的差别对“四个地区的饮食区别”这一结果的影响程度是不一样的。

那么，我们就需要把哪些对“四个地区的饮食区别”影响最大的食物给筛选出来，然后排序，并依据相关公式进行整合计算，计算出一个新的参数——也就是主成分。排名第一的主成分对整体变异度的贡献最大。

对上述数据进行PCA分析，结果如下：

image-20190804071350641

这是第一主成分的区分情况，我们再加上第二个维度：PC2

image-20190804071408923

PC1对四个地区的区分结果最好，PC2就看不出有什么区分结果了。当然，数据不同，PC1和PC2的贡献度也不同，有时候PC2对数据的区分力度也很大，要看具体的数据.

比如有8种样本，6个是近似病理类型的，另外两个是另一大类的疾病类型的样本。这时候对8个样本的蛋白质组学数据进行PCA分析，那么，鉴定到的数千个蛋白的丰度信息经过PCA的计算，其结果一定是某6个聚集在一起，另外2个聚集在一起。

假设样本中平均鉴定到3000个蛋白，那么也就意味着本数据有3000个维度，每个蛋白（维度）的相对定量信息的变化在这8个样本中的分布都是不一样的，有些蛋白在8个样本中变异很小，有些很大，有些居中。那么，应用PCA进行降维，将3000个维度降至2-3个维度，从而简化了数据，突出了主要矛盾。参考:

http://bbs.foodmate.net/thread-1071830-1-1.html

那么PC1和PC2是什么呢?上面举例中,比如3000个蛋白,3000个维度,PC1是多少个维度的组合呢?是哪些个维度的组合呢?下面这张R语言实战中的图片可以回答,并不是某一个维度,而是

image-20190804072619643

PCA的目标是用一组较少的不相关变量代替大量相关变量，同时尽可能保留初始变量的信息，这些推导所得的变量称为主成分，它们是观测变量的线性组合。如第一主成分为

image-20190804073137978

它是k个观测变量的加权组合，对初始变量集的方差解释性最大。第二主成分也是初始变量的线性组合，对方差的解释性排第二，同时与第一主成分正交（不相关）。后面每一个主成分都最大化它对方差的解释程度，同时与之前所有的主成分都正交.目标是能用较少的主成分来解释全部变量。

####PCA看分组情况
library("FactoMineR")
library("factoextra") 
df.pca <- PCA(t(ex), graph = FALSE)
fviz_pca_ind(df.pca,
             geom.ind = "point",
             col.ind = group_list, 
             addEllipses = TRUE, 
             legend.title = "Groups"
)
save(ex,pd,group_list,file='ex_pd.Rdata')

当处理组和未处理组分得很开时,说明实验分组

处理组和未处理组分得很开

3.差异分析

下面是GSE5282的数据集,整理好所得到的表达矩阵和分组矩阵.

image-20190804073913072

image-20190804074017364

表达矩阵中的样本名的顺序要和分组信息中的顺序一一对应上,非常重要

rm(list = ls()) 
options(stringsAsFactors = F)
load('ex_pd.Rdata')
ex <- ex[,1:6]
#####可以自己生成，group_list不可以有空格；
group_list <- group_list[1:6]
group_list<- gsub(' ','_',trimws(group_list))#
###构建实验设计矩阵
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(ex)
library(limma)#limma
####构建比较矩阵,设置用来对比的组别
contrast.matrix<-makeContrasts(contrasts=c('EGF_4h-Control_4h'),levels = design)
######limma三部曲,只需要归一化后的数据、实验矩阵、比较矩阵的输入
fit <- lmFit(ex,design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1)
####上面得到了p值等统计的结果，topTable对p值校验，对基因排序
tT <- topTable(fit2, adjust="fdr", number=nrow(fit2))
dim(tT)
####提取出我们需要的指标
tT <- subset(tT, select=c('adj.P.Val',"P.Value","logFC"))
save(tT,file='DEG.Rdata')

>**多重检验时，产生adjust.pvalue的需求**

接下来,进行id转换

image-20190804214503093

同上面的一样,在http://www.bio-info-trainee.com/1399.html这个网址内找到GPL1355对应的R包注释集,是rat2302,下载的时候加上'db'.

BiocManager::install("rat2302.db",ask = F,update = F)
library('rat2302.db')
load('DEG.Rdata')
#####
ls("package:rat2302.db")
p2se<- select(rat2302.db,keys = keys(rat2302.db),columns = c('SYMBOL','ENTREZID'))#select这个函数可以直接将包中的'SYMBOL','ENTREZID'提取出来,相比用to table 方便很多
tT$probe_id <- rownames(tT)
tT <- merge(tT,p2se,by.x='probe_id',by.y='PROBEID')#大小写不同,可以用by.x,by.y,对应上
save(tT,file='DEG_symbol.Rdata')

4.差异分析后的可视化-热图、火山图、富集分析

热图

load('ex_pd.Rdata')
load('DEG.Rdata')
ex <- ex[,1:6]
save(ex,pd,group_list,file='ex_pd.Rdata')
load('ex_pd.Rdata')
######pheatmap 这是取前100个gene进行展示，可根据自己的需求进行调整
choose_gene<-rownames(tT)[1:100]
#####用normalize后的数据进行展示
data_matrix <-ex[choose_gene,]
#####热图更详细的了解https://www.jianshu.com/p/0e5ce59fa79e
#####scale对数据处理，求得的结果为z-score，即数据与均值之间差几个标准差 
####t是对数据做转置处理，为了适应scale的要求
data_matrix=t(scale(t(data_matrix)))
#####查看scale处理后数据的范围
fivenum(data_matrix)
####目的是避免出现极大极小值影响可视化的效果
data_matrix[data_matrix>1.2]=1.2
data_matrix[data_matrix< -1.2] = -1.2
library(pheatmap)
pheatmap(data_matrix)
####调整下颜色，使正负值颜色的对比会更加鲜明
require(RColorBrewer)
bk = c(seq(-1.2,1.2, length=100))
annotation_col = data.frame(Group = rep(c("EGF_4h", "Control_4h"), c(3,3)))#annotation_col
rownames(annotation_col)<-colnames(ex)
####annotation_col和annotation_row的格式需为数据框
####breaks参数用于值匹配颜色
####看下，color和breaks的对应，进行理解
pheatmap(data_matrix,
         breaks=bk,#从bk开始
         show_rownames = F,#默认是展示
         annotation_col = annotation_col,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
         filename = 'test.png')
####可以调整的内容有是否聚类、是否分组、是否显示行和列的内容等

image-20190804220955224

test.png

火山图

load(file = "DEG_symbol.Rdata")
####主要用两个值，logFC和p.value
####可参考进行logFC_cutoff的计算
logFC_cutoff <- with(tT,mean(abs( logFC)) + 2*sd(abs( logFC)) )
###logFC_cutoff也可自行根据需求设置
####依据logFC和adj.P.val为火山图的颜色分组做准备
tT$change = ifelse(tT$adj.P.Val < 0.05 & abs(tT$logFC) > logFC_cutoff,
                         ifelse(tT$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
####记录上调、下调基因数目，作为title的内容
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(tT[tT$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(tT[tT$change =='DOWN',]))

rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load(file = "DEG_symbol.Rdata")
####主要用两个值，logFC和p.value
####可参考进行logFC_cutoff的计算
logFC_cutoff <- with(tT,mean(abs( logFC)) + 2*sd(abs( logFC)) )
###logFC_cutoff也可自行根据需求设置
####依据logFC和adj.P.val为火山图的颜色分组做准备
tT$change = ifelse(tT$adj.P.Val < 0.05 & abs(tT$logFC) > logFC_cutoff,
                         ifelse(tT$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
####记录上调、下调基因数目，作为title的内容
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(tT[tT$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(tT[tT$change =='DOWN',]))

library(ggplot2)
g = ggplot(data=tT, 
           aes(x=logFC, y=-log10(P.Value), color=change)) +
    geom_point(alpha=0.4, size=1.75) + #出点状图
    theme_set(theme_set(theme_bw(base_size=20)))+ 
    xlab("log2 fold change") + ylab("-log10 p-value") + 
    ggtitle( this_tile ) + 
    theme(plot.title = element_text(size=15,hjust = 0.5))+
    scale_colour_manual(values = c('blue','black','red'))+
   annotate('text',x=tT$logFC[tT$logFC>5],y=-log10(tT$P.Value[tT$logFC>5]),label=tT$SYMBOL[tT$logFC>5])
print(g)
ggsave(g,filename = 'volcano.png')
save(tT,file='DEG_change.Rdata')

火山图

load('DEG_change.Rdata')
tT <- na.omit(tT)
####富集前的基因处理
library(clusterProfiler)
gene_up= tT[tT$change == 'UP','ENTREZID'] 
gene_down=tT[tT$change == 'DOWN','ENTREZID'] 
gene_diff=c(gene_up,gene_down)
gene_all = tT$ENTREZID
###############################KEGG富集分析#####################################################
## 1.KEGG pathway analysis
#上调、下调、差异、所有基因
#clusterProfiler作kegg富集分析：
save(kk.down,file='down.Rdata')
kk.up<- enrichKEGG(gene         = gene_up,
                   organism     = 'rno',
                   universe     = gene_all,
                   pvalueCutoff = 1)
kk.down <- enrichKEGG(gene         =  gene_down,
                      organism     = 'rno',
                      universe     = gene_all,
                      pvalueCutoff = 1)
#######################下述ggplot代码用于画图########################################
down_kegg <- kk.down@result
down_kegg <- down_kegg[down_kegg$pvalue<0.05,]
library(ggplot2)
ggplot(down_kegg, aes(x=Description, y=Count,fill=pvalue)) + 
  geom_bar(stat="identity") + 
  scale_fill_gradient(low = 'blue',high='red')+
  scale_x_discrete(name ="pathway names") +
  scale_y_continuous(name ="Count") +
  coord_flip() + 
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  ggtitle("Pathway Enrichment")

image-20190804235105060

image-20190804234953521

KEGG注释

#首先是KEGG的注释
rm(list = ls())  ## 魔幻操作，一键清空~
options(stringsAsFactors = F)
load('DEG_change.Rdata')
tT <- na.omit(tT)
####富集前的基因处理
library(clusterProfiler)
gene_up= tT[tT$change == 'UP','ENTREZID'] 
gene_down=tT[tT$change == 'DOWN','ENTREZID'] 
gene_diff=c(gene_up,gene_down)
gene_all = tT$ENTREZID
#######KEGG富集分析
## 1.KEGG pathway analysis
#上调、下调、差异、所有基因
#clusterProfiler作kegg富集分析：
save(kk.down,file='down.Rdata')
kk.up<- enrichKEGG(gene         = gene_up,
                   organism     = 'rno',
                   universe     = gene_all,
                   pvalueCutoff = 1)
kk.down <- enrichKEGG(gene         =  gene_down,
                      organism     = 'rno',
                      universe     = gene_all,
                      pvalueCutoff = 1)
#######################下述ggplot代码用于画图########################################
down_kegg <- kk.down@result
down_kegg <- down_kegg[down_kegg$pvalue<0.05,]
#down_kegg
up_kegg<- kk.up@result
up_kegg <- up_kegg[up_kegg$pvalue<0.05,]
library(ggplot2)
ggplot(down_kegg, aes(x=Description, y=Count,fill=pvalue)) + 
  geom_bar(stat="identity") + 
  scale_fill_gradient(low = 'blue',high='red')+
  scale_x_discrete(name ="pathway names") +
  scale_y_continuous(name ="Count") +
  coord_flip() + 
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  ggtitle("Pathway Enrichment")
#up_kegg
library(ggplot2)
ggplot(up_kegg, aes(x=Description, y=Count,fill=pvalue)) + 
  geom_bar(stat="identity") + 
  scale_fill_gradient(low = 'blue',high='red')+
  scale_x_discrete(name ="pathway names") +
  scale_y_continuous(name ="Count") +
  coord_flip() + 
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  ggtitle("Pathway Enrichment")

image-20190805001223395

down_kegg

image-20190805001139289

up_kegg

富集程度看颜色深浅,而柱状图长短代表Count数,最后出图就是上面富集分析(enrichKEGG函数)可视化后的展示)

down_kegg-dotplot

GO注释

主要包括这三个方面的注释

image-20190805002532500

#代码如下,可画图展示
library(org.Rn.eg.db)
ego_CC <- enrichGO(gene = gene_diff,
                   OrgDb= org.Rn.eg.db,
                   universe = gene_all,
                   ont = "CC",
                   pAdjustMethod = "BH")
#生物过程
ego_BP <- enrichGO(gene = gene_diff,
                   OrgDb= org.Rn.eg.db,
                   universe = gene_all,
                   ont = "BP",
                   pAdjustMethod = "BH")
#分子功能：
ego_MF <- enrichGO(gene = gene_diff,
                   OrgDb= org.Rn.eg.db,
                   universe = gene_all,
                   ont = "MF",
                   pAdjustMethod = "BH")
#可视化展示-Barplot
barplot(ego_CC, showCategory=12,title="bar_CC")
barplot(ego_BP, showCategory=12,title="bar_BP")
barplot(ego_MF, showCategory=12,title="bar_MF")
#可视化展示-Dotplot
dotplot(ego_CC,showCategory = 12,title="dot_CC")
dotplot(ego_BP,showCategory = 12,title="dot_BP")
dotplot(ego_MF,showCategory = 12,title="dot_MF")
save(kk.up,kk.down,ego_BP,ego_CC,ego_MF,file='annotation.Rdata')