3.1 蛋白互作网络分析

蛋白质的相互作用(PPI, Protein-protein interaction)是指蛋白质分子之间存在相互作用, 可以从生物化学、信号转导和遗传 网络的角度研究这种相关性。
PPI网络中的每个结点代表一种蛋白质,节点之间有连线代表它 们之间存在相互作用关系。

#安装STRINGdb 软件包
source("https://bioconductor.org/biocLite.R")
biocLite("STRINGdb")#用于连接STRING数据库
install.packages("rlist")
install.packages("STRINGdb")
library("STRINGdb")
library("rlist")
#设置工作目录
setwd("D:/train/TCGA/lab/PPI")
getwd()
#找到源文件
deg_file<-"/train/TCGA/lab/DEG/Gene/DE_genes.txt"

image.png

# 获取物种的分类编号
# get_STRING_species(version="10", species_name=NULL) 
# 9606 代表人类,NBCI的物种分类编号
string_db <- STRINGdb$new(version="10", species=9606,
                          score_threshold=700, input_directory= "D:/train/TCGA/lab/PPI")
# 读取差异表达的文件，获得差异表达基因列表
degs = read.table(deg_file,header=T,comment.char = "",check.names=F)
degs$gene <- rownames(degs)
head(degs)

image.png

# 查看有多少差异表达的基因需要分析 
cat("Total deg genes:", dim(degs)[1])

image.png

# 将基因的ID map 到string 数据库中， 不一定每个基因都能map上
deg_mapped <- string_db$map( degs, "gene", removeUnmappedRows = TRUE )

image.png

可以看到少了一百二十个左右，并不是所有基因都能映射的上

# 筛选出一部分结果，进行绘图,选择了前两百个，当然可以任意
hits <- deg_mapped$STRING_id[1:200]

options(par(mar=c(2.1, 0.1, 4.1, 2.1)))
# 绘图
network_file = "network_plot.png"
png(file=network_file, height = 900, width = 1600)
string_db$plot_network( hits,required_score = 400)
dev.off()

image.png

圆圈代表基因，线条代表基因间存在 相互作用，圆圈内部的结果代表蛋白 的结构。线头颜色代表证明蛋白之间 存在相互作用的不同证据

第二种方法：采用cytoscape 进行网络分析

# 将所有的结果输出到文件，后面采用cytoscape 进行网络分析
info <- string_db$get_interactions(deg_mapped$STRING_id)
write.table(info, file = "STRING_info.txt",sep="\t", row.names =F, quote = F)

image.png

在cytoscape软件中点击File→Import→Network→File
按下图操作

image.png
image.png

点击OK
结果呈现：

image.png
图里面就是完整的蛋白互作网络分析。

3.2 KEGG分析

在生物体内，不同基因相互协调行使其生物学功能，利用通路分析可以更进一步了解基因的参与代谢通路及具体的生物学功能。
Pathway富集分析以KEGG Pathway为单位，应用统计学检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway。

# 分析之前，先安装Bioconductor 和 clusterProfiler，org.Hs.eg.db
source("https://bioconductor.org/biocLite.R")
biocLite("org.Hs.eg.db")
biocLite("clusterProfiler")
biocLite("pathview")
############################################
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("clusterProfiler")
#因为最上面的语句我一直安装不成功，所以找到了这个办法，看起来成功了
#语句来源：http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("pathview")
################################################
library("clusterProfiler")
library("org.Hs.eg.db")
library("pathview")
#设置工作目录
setwd("D:/train/TCGA/lab/Annotation")
deg_file<-("/train/TCGA/lab/DEG/Gene/DE_genes.txt")
# 读取差异表达的文件，获得差异表达基因列表
degs = read.table(deg_file,header=T,comment.char = "",check.names=F)
DEG_list <- rownames(degs)
head(DEG_list)
# KEGG pathway 富集
gene_ids <- bitr(DEG_list, fromType="ENSEMBL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
#此处将DEG_LIST中的ID从ENSEMBL格式转化成ENTREZID格式
kegg <- enrichKEGG(gene         = gene_ids$ENTREZID,
                   organism     = 'hsa',
                   pvalueCutoff = 0.05)
#hsa代表的是人类
head(kegg)
write.table(kegg, file = "KEGG_enrichment_stat.txt",sep="\t", row.names =F, quote = F)

image.png

这个表目前还看不出什么有效信息，需要依重要性进行排序

# 点图
png(file = "kegg_dotplot.png")
dotplot(kegg,title="Enrichment KEGG_dot")
dev.off()

image.png

从图中信息可以看出这种癌症可能和碳代谢有很大关系

# 合并表达信息和基因信息
#知道哪些基因上调或下调
#logFC，Fold_Change为差异表达基因上调倍数，然后取Log，LogFC一般表达相差两倍以上是有意义的
#放宽至1.5或者1.2倍也可以接受
deg_info <- degs['logFC']
deg_info['ENSEMBL'] <- rownames(deg_info)
gene_ids_merge <- merge(gene_ids, deg_info,by='ENSEMBL')
# 去掉ENTREZID 重复的基因
index<-duplicated(gene_ids_merge$ENTREZID)
gene_ids_merge<-gene_ids_merge[!index,]
# 提取FC的值，在map上进行颜色标注
map_ids <- as.matrix(gene_ids_merge['logFC'])
rownames(map_ids) <- gene_ids_merge$ENTREZID

#  绘制富集的pathway 以hsa00280 为例
pathway_id <- "hsa00280"
map <- pathview(gene.data  = map_ids[,1],
                pathway.id = pathway_id,
                species    = "hsa", kegg.native = TRUE)

image.png
image.png

3.3GO富集分析

GO（Gene ontology）按照生物途径(Biology Process)，分子功能（Molecular Function）和细胞 定位(Cellular Location)对基因进行注释和分类。通 过对差异表达基因进行GO terms富集度统计学的分 析，计算出差异基因GO term的p-value和p-value的 FDR值（q-value），定位差异基因最可能相关的GO term。

#GO
#ont有三种 :BP(Biology Process),CL（Molecular Function）,MF(Cellular Location)
ego <- enrichGO(gene          = DEG_list,
                keyType       = "ENSEMBL",
                OrgDb         = org.Hs.eg.db,
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.05,
                readable      = TRUE)

write.table(ego, file = "GO_enrichment_stat.txt",sep="\t", row.names =F, quote = F)

看下文件内容

image.png

ego_results<-summary(ego)
ego_results

image.png
image.png
image.png

# 绘制相关图
pdf(file = "ego_barplot.pdf")
barplot(ego, showCategory=20, x = "GeneRatio")
#画前面20个
dev.off()

image.png

# 点图
pdf(file = "ego_dotplot.pdf")
dotplot(ego,showCategory=20)
dev.off()

image.png

# 绘制topGO图
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("topGO")
library("topGO")
pdf(file = "topgo.pdf")
plotGOgraph(ego)
dev.off()

image.png