R语言轻松实现生物信息分析——GEO数据库之绘图

在上周，小编和其他同事一起学习了GEO数据库的基础知识，主要的是基因的注释，基因的合并，预处理以及基因差异性分析，GO分析，我们一起来学习一下绘图的部分吧。大家是不是也很好奇呢，先给大家看看效果图吧。

怎么样？是不是一下就有格调啦，下面我们来看看如何实现吧。

BiocManager::install("limma")
library(limma)
BiocManager::install("pheatmap")
library(pheatmap)
BiocManager::install("ggplot2")
library(ggplot2)

加载包之后，我们将需要的数据合并，并读取此数据。

rt<-read.csv("normalize.csv")
head(rt)
dim(rt)

# 建立分组，分为normal和abnormal两组，总共是60个样本。

group <- c(rep("normal",25),rep("abnormal",25),rep("normal",5),rep("abnormal",5))
dim(group)
#group <- factor(pData(eset.rma)$treatment)
group <- factor(group,levels = c("normal","abnormal"),ordered = F)

分组完成后构建分组比较矩阵

design <- model.matrix(~ group)
design

之后对矩阵命名进行比较

colnames(design) <- levels(group)
dim(design)

构建对比模型，比较两个实验条件下表达数据，

contrast.matrix <- makeContrasts(contrasts = "normal - abnormal",
levels=design)
contrast.matrix

进行线性模型的拟合，

exprSet_new<-rt[,3:ncol(rt)]
fit <- lmFit(exprSet_new, design)

根据对比模型进行差值计算

fit1 <- contrasts.fit(fit, contrast.matrix)

进行贝叶斯检验

fit2 <- eBayes(fit1)

对于"abnormal - normal"比较组（coef=1），表达变化超过2倍（lfc参数）的基因数为299个，而变化超过2倍、p<0.01的基因总数为224个

生成所有基因的检验结果报告(此时无筛选条件)

dif_gene0 <- topTable(fit2, coef="normal - abnormal", adjust.method="fdr", number = nrow(fit2))
nrow(dif_gene0)

加入筛选条件lfc=2

dif_gene1 <- topTable(fit2, coef="normal - abnormal", adjust.method="fdr", lfc=2, number = nrow(fit2))

nrow(dif_gene1)

dif_gene2 <- topTable(fit2, coef="normal - abnormal", adjust.method="fdr", p.value=0.05, lfc=1.5, number = nrow(fit2))

#加入筛选条件lfc=2，p.value=0.01
nrow(dif_gene2)
class(dif_gene2)
names(dif_gene2)

记得一定要保存全部基因结果，

dif_gene0$gene_names <- rownames(dif_gene0)
dim(dif_gene0)
write.table(dif_gene0,file="dif_gene0.txt",sep="\t",quote=F,row.names = F)

保存显著性差异基因，

dif_gene2$significant <- as.factor(dif_gene2$P.Value<0.01)
dif_gene2$gene_names <- rownames(dif_gene2)
dim(dif_gene2)
write.table(dif_gene2,file="dif_gene2.txt",sep="\t",quote=F,row.names = F)

上面就是我们必要的处理，之后我们才能进行差异基因热图的绘制。

第一步首先是设定差异基因阈值，减少差异基因用于提取表达矩阵。

dif_g<- topTable(fit2, coef="normal - abnormal", adjust.method="fdr", p.value=0.05, lfc=1, number = nrow(fit2))
nrow(dif_g)

现在我们提出前部分数据用作热图绘制，

heatdata <- exprSet_new[rownames(dif_g),]
head(heatdata)
dim(heatdata)

记住要制作一个分组信息用于注释，这是非常必要的，因为你进行GO分析的时候会用到，不然GO分析会报错。

annotation_col <- data.frame(group)
rownames(annotation_col) <- colnames(heatdata)

对了，告诉大家一个关于注释的调整方法，如果注释出界，可以通过调整格子比例和字体修正，代码如下：

pheatmap(heatdata,cluster_rows = TRUE,cluster_cols = TRUE,annotation_col =annotation_col, annotation_legend=TRUE,show_rownames = F,show_colnames = F,scale = "row", color =colorRampPalette(c("blue", "white","red"))(100))

如果大家对这段代码有兴趣，可以自己去加深一些了解。数据整理好了，我们来看看火山图吧。绘制火山图的代码如下：

dif_gene_volcano<-read.table("dif_gene0.txt",header=T)
names(dif_gene_volcano)
View(dif_gene_volcano)
dif_gene_volcano$change = as.factor(ifelse(dif_gene_volcano$adj.P.Val < 0.05 &
abs(dif_gene_volcano$logFC) > 2,
ifelse(dif_gene_volcano$logFC < -2,'UP','DOWN'),'NOT'))

这里强调一点儿，不同的阈值，筛选到的差异基因数量不一样。

table(dif_gene_volcano$change)
this_title = paste0('Cutoff for logFC is ',round(2,2),'\nThe number of up gene is',nrow(dif_gene_volcano[dif_gene_volcano$change == 'UP',]),'\nThe
number of down gene is ',nrow(dif_gene_volcano[dif_gene_volcano$change == 'DOWN',]))
g<-ggplot(data=dif_gene_volcano, aes(x=logFC, y =-log10(adj.P.Val),color=change)) +geom_point(alpha = 0.4, size = 1.75) +
theme_set(theme_set(theme_bw(base_size = 10))) +
xlab("log2 fold change") + ylab('-log10 p-value') +
ggtitle(this_title) +theme(plot.title = element_text(size = 10,hjust = 0.5)) +
scale_color_manual(values = c('blue','black','red'))
print(g)

这样就可以出来火山图啦。

下面我们来看一下GO分析的图怎么绘制，首先我们要将差异性基因的行号匹配成基因名称，

id <- read.csv("dif_gene2.csv",header = TRUE)
id[8]
id1 <- read.csv("normalize.csv",header = TRUE,stringsAsFactors = FALSE)
res <- merge.data.frame(id,id1,by.x = 'gene_names',by.y = 'number')
head(res)
gene_names<-res[,9]
write.csv(gene_names,"gene_names.csv")

第一步就完成了，之后我们加载R包，

BiocManager::install("org.Hs.eg.db")
library(DOSE)
library(GO.db)
library(org.Hs.eg.db)#人类注释包
library(topGO)
library(GSEABase)
library(clusterProfiler)

利用bitr函数将基因名称转换为ENTREZID号；物种是人org.Hs.eg.db；#可能会有部分基因对应不到ENTREZID

gene=read.csv("gene_names.csv",header = F)
gene=as.character(gene[,1])
gene.df <- bitr(gene, fromType ="SYMBOL",toType = c("ENSEMBL","ENTREZID","UNIPROT"), OrgDb = org.Hs.eg.db)head(gene.df)
id = gene.df[,3]

fromType是指你的数据ID类型是属于哪一类的，toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种，Orgdb是指对应的注释包是哪个。

现在我们进行最后一步GO分析，

ego_CC <- enrichGO(gene = id,
OrgDb=org.Hs.eg.db,
ont = "CC",#ont代表GO的3大类别，BP, CC, MF
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
write.csv(as.data.frame(ego_CC),row.names = F, file = "ego_CC.csv")
barplot(ego_CC,drop = TRUE,title = "enrichment_CC",showCategory = 12)

这样我们就做出来柱状图啦，以上就是所有的内容啦，期待大家的下期相约哟。