rm(list = ls())library(data.table)library(magrittr)expr <- fread('TCGA-LIHC.htseq_fpkm.tsv.gz',h=T,check.names = F)ann <- read.table('gencode.v22.annotation.gene.probeMap',h=T)expr2 <- merge(ann[,1:2],expr,by=1) %>% .[,-1]library(limma)expr2 <- avereps(expr2[,-1],ID=expr2$gene) %>% as.data.frame() expr3 <- as.data.frame(t(expr2))expr4 <- expr3[order(names(expr3))] 相關性分析常用的有pearson相關係數和spearman相關係數，其中pearson相關係數適用於連續性變量，且變量服從正態分布的情況，為參數性的相關係數；而spearman相關係數適用於連續性及分類型變量，為非參數性的相關係數。以下採用spearman相關。
gene <- as.numeric(expr4[,"A1BG"]) col <- colnames(expr4)data <- data.frame(col)for (i in 1:length(col)){ test <- cor.test(as.numeric(expr4[,i]),gene,type="spearman") data[i,2] <- test$estimate data[i,3] <- test$p.value }這裡的p值<0.05會有很多基因滿足，我們需要進一步根據篩選出來的基因個數來設置相關係數cor的閾值，這裡小王子設置成0.6，最終得到87個顯著相關基因進行下遊分析。
names(data) <- c("gene","correlation","pvalue")sig_data <- data[data$pvalue<0.05&abs(data$correlation)>0.6,]sig_data <- na.omit(sig_data)四、利用org.Hs.eg.db和clusterProfiler包進行GO和KEGG富集分析：這兩個包都不怎麼好安裝，大家如果沒有的話可以嘗試下小王子隱藏的代碼。最終獲得與我們所選目的基因顯著相關的基因在GO和KEGG資料庫中富集的通路。
######利用顯著相關基因進行GO和KEGG富集分析######library("org.Hs.eg.db")###沒有的話以下代碼安裝#if (!requireNamespace("BiocManager", quietly = TRUE))# install.packages("BiocManager")#BiocManager::install("org.Hs.eg.db")library("clusterProfiler")###沒有的話以下代碼安裝#if (!requireNamespace("BiocManager", quietly = TRUE))# install.packages("BiocManager")#BiocManager::install("clusterProfiler")gene_name=as.vector(sig_data$gene)geneID <- mget(gene_name, org.Hs.egSYMBOL2EG, ifnotfound=NA)%>%as.character()data=cbind(Gene=sig_data,entrezID=geneID)%>%as.data.frame()go <- enrichGO(gene=data$entrezID, OrgDb = org.Hs.eg.db, ont='ALL',pvalueCutoff = 0.05)%>%as.data.frame()kegg <- enrichKEGG(gene = data$entrezID,organism ="human",pvalueCutoff = 0.05)%>%as.data.frame()由於GO富集到了275條通路，這裡根據p.adjust選取前20條進行展示；而KEGG一共富集到9條密切相關通路，這裡完全展示。
go1 <- go[order(go$p.adjust)[1:20],] go1$p.adjust <- -log10(go1$p.adjust)go1$GeneRatio <- ifelse(go1$ONTOLOGY=="BP",go1$Count/4541*100,go1$Count/4709*100)library(ggplot2)ggplot(go1,aes(Description,Count,fill=p.adjust))+ geom_bar(stat = "identity")+ coord_flip()+ scale_fill_gradient(low="blue",high="red")+ labs(x="",fill="-log10(adj.P)")+ ggtitle("Top20 GO enrichment of our genes")+ theme_bw()+ theme(axis.title = element_text(size = 16), axis.text = element_text(size = 12), legend.title = element_text(size = 12), plot.title = element_text(size = 18,face = "bold",hjust = 0.5))######繪製KEGG水平柱狀圖######kegg$p.adjust <- -log10(kegg$p.adjust)##將p.adjust進行-log10換算library(ggplot2)ggplot(kegg,aes(Description,Count,fill=p.adjust))+ geom_bar(stat = "identity")+ coord_flip()+ scale_fill_gradient(low="blue",high="red")+ labs(x="",fill="-log10(adj.P)")+ ggtitle("KEGG enrichment of our genes")+ theme_bw()+ theme(axis.title = element_text(size = 16), axis.text = element_text(size = 12), legend.title = element_text(size = 12), plot.title = element_text(size = 18,face = "bold",hjust = 0.5))#題目設置：face = "bold"加粗字體；hjust = 0.5位置居中##kegg為輸入數據框，將Description作x，Count作y，繪圖後再利用coord_flip()進行x、y軸轉換，實現水平柱狀圖的繪製；##scale_fill_gradient(low="blue",high="red")設置漸變色（越紅P值越顯著）；ggtitle("")設置題目；theme設置主題背景、字體等。ggplot(go1,aes(GeneRatio,reorder(Description,GeneRatio),color=p.adjust,size=Count))+ geom_point(shape=19)+ scale_color_gradient(low = "green",high = "red")+ scale_size_continuous(range = c(2,9))+ labs(x="GeneRatio(%)",y="",color="-log10(adj.P)",size="Gene Count")+ ggtitle("Top20 GO enrichment of our genes")+ theme_bw()+ theme(axis.title = element_text(size = 16), axis.text = element_text(size = 12), legend.title = element_text(size = 12), plot.title = element_text(size = 18,face = "bold",hjust = 0.5))最後幾張結果展示部分，KEGG結果不太美觀，當然與我們選取的基因有很大關係，不過你也可以設置篩選條件時選擇cor=0.5作為閾值，會篩選出657個顯著相關基因，或許結果又不一樣了呢？多嘗試幾次，總有理想結果。或許還有人有這樣的疑惑：正相關和負相關的基因能否放一起進行富集分析呢？這裡小王子收集到的資料採用放在一起分析，映射到差異分析上調和下調放一起富集分析或許大家更好理解一些。
1.R-可視化基礎(4)—— GO/KEGG柱狀&氣泡圖--《醫學僧的科研日記》
2.R語言：單基因批量相關性分析的妙用--《愛科學》