热图、韦恩图、GO富集分析图（有了转录组数据不知道该怎么写文章，看我就对了！）

你好啊，我是Christine，这里是周一数据挖掘专栏，欢迎和我们一起学习~

上个月我们曾经学习过一篇细胞内PD-L1调控RNA稳定性的文章：TCGA正常血液样本中PD-L1与BRCA1和NBS1的表达量相关性（没找到原数据，但重复出了结果）(PS,目前为止我还是没能找到TCGA中"Blood Derived Normal"样本的RNA-seq数据 ╯︿╰ )，当时文中Fig4的测序数据还没有公开，现在GEO上已经可以下载了，不过作者只提供了处理好的excel表格数据和RNA-seq raw-data，对处理RNA-seq raw data感兴趣的欢迎点击“阅读原文”去B站学习Jimmy老师的视频，但是需要有linux服务器。

本周我们就系统性的学习一下本文的4张图（几乎全部结论）的画法吧！

A. 分别用PD-L1抗体或IgG进行RNA免疫共沉淀(RIP)，热图展示DNA损伤相关基因的表达量

B. 分别用2种shRNA沉默PD-L1基因，热图同样为DNA损伤相关基因的表达量

C. 图A与图B的重叠基因及重叠的DNA损伤基因，也就是受PD-L1调控的RNAs

D. RIP-seq及RNA-seq重叠基因的GO富集分析

热图

GSE128613和GSE128742分别为RIP和RNA-seq数据，GEO网站上提供了xlsx格式的数据

#### RIP peak-score
library(readxl)
RIP <- read_excel("./GSE128613_processed_RIP_data.xlsx")
# DNA damage基因，在excel的第3个sheet中
RIP_gene <- read_excel("./GSE128613_processed_RIP_data.xlsx",sheet = 3,col_names = F)
setdiff(RIP_gene$...2,RIP$`Gene name`) 
#基因名有重叠现象,"USP1,UBP1" "BRIP1,BACH1" "POLQ,POLH"
#逗号前后是不同的基因，查了一下HGNC号，作者用的应该是前一个
choose_matrix <- as.matrix(RIP[RIP$`Gene name` %in% c(RIP_gene$...2,"USP1","BRIP1","POLQ"),3:8]) 
colnames(choose_matrix) <- c("IgG_1","IgG_2","IgG_3","PDL1_1","PDL1_2","PDL1_3")
choose_matrix <- log2(choose_matrix + 1)
library(pheatmap)
pheatmap(choose_matrix, scale = "row",
         cluster_rows = F, cluster_cols = F,
         show_colnames = T,angle_col=45)

#### RNA-seq
RNA_seq <- read_excel("./GSE128742_processed_RNA-seq_data.xls")
# DNA damage基因
RNA_seq_gene <- read_excel("./GSE128742_processed_RNA-seq_data.xls",sheet = 3,col_names = F)
setdiff(RNA_seq_gene$...2,RNA_seq$`gene name`) 
#基因名也有重叠,"USP1,UBP1" "POLQ,POLH" "MED1,MBD4"

choose_matrix <- as.matrix(RNA_seq[RNA_seq$`gene name` %in% c(RNA_seq_gene$...2,"USP1","POLQ","MED1"),2:7]) 
colnames(choose_matrix) <- c("Ctrl shRNA(1)","Ctrl shRNA(2)",
                             "shPD-L1-1(1)","shPD-L1-1(2)",
                             "shPD-L1-2(1)","shPD-L1-2(2)")
choose_matrix <- log2(choose_matrix + 1)
pheatmap(choose_matrix, scale = "row",
         cluster_rows = F, cluster_cols = F,
         show_colnames = T)

heatmap.png

韦恩图

## PD-L1敲除后下调的基因与RIP-seq的重叠基因
library("VennDiagram")
VENN.LIST=list(RNA_seq = RNA_seq$`gene name`,
               RIP = RIP$`Gene name`)
venn.plot <- venn.diagram(VENN.LIST, NULL,
                          fill=c("RoyalBlue", "Salmon"), 
                          alpha=c(0.8,0.8), 
                          cex = 2,cat.cex=1.5)
grid.draw(venn.plot) 

# PD-L1敲除后下调的DNA damage基因与RIP-seq的DNA damage基因的重叠情况
VENN.LIST=list(RNA_seq = RNA_seq_gene$...2,
               RIP = RIP_gene$...2)
venn.plot <- venn.diagram(VENN.LIST, NULL,
                          fill=c("RoyalBlue", "Salmon"), 
                          alpha=c(0.8,0.8), 
                          cex = 2,cat.cex=1.5)
grid.draw(venn.plot)

GO富集分析图

图中展示的是受PD-L1调控的基因的GO功能富集分析，也就是：与PD-L1抗体结合的基因和PD-L1敲除后下调的基因的交集。

GO富集分析的R包和网页工具很多，我习惯用clusterProfiler包，因为它提供了很多画图函数，比较方便。

genes <- intersect(RIP$`Gene name`,RNA_seq$`gene name`) #1756个基因
library(clusterProfiler)
library(org.Hs.eg.db)
res <- enrichGO(genes, 
                OrgDb = org.Hs.eg.db, 
                keyType = "SYMBOL", 
                ont = "BP",
                pvalueCutoff = 0.05
                )
head(res@result,n=20) # 查看结果
# 画图（这个R包提供了多种绘图函数，感兴趣的可以看看说明书自行尝试）
barplot(res, showCategory=15) #横轴显示的是富集的基因数目

因为用的方法不同，和原文中富集的terms有一定出入，但总体结果是类似，都在细胞代谢、转录、蛋白质修饰等途径中有富集。

作者是在Gene Ontology Consortium website上做的分析，而且文章的图看起来是作者手动选择了一些自己感兴趣的显著terms（所以要搞清楚自己课题的需求）

# 将需要的基因写入文件
write.table(genes,file = "./GO_genes.txt", quote = F, row.names = F, col.names = F)

# 将基因上传到http://geneontology.org/上，下载分析结果，读入R (记得删除文件最上面几行的注释)
df_go <- read.delim("./analysis.txt",header = T, stringsAsFactors = F)

# 挑选文章中用的GO-terms
df <- df_go[c(274,174,198,378,269,182,138,305,107,163,395,81,116,75,393),c(1,7)]
df$upload_1..raw.P.value. <- -log10(df$upload_1..raw.P.value.)
colnames(df) <- c("terms","val")

# 按照pvalue排序画barplot
library(forcats)
library(ggplot2)
ggplot(df, aes(x=fct_reorder(terms,val), y=val)) +
  geom_bar(stat="identity") +
  labs(x=" ",y="-log10(p-val)") +
  coord_flip() #翻转x,y轴

这张图看起来和文中的差不多了，虽然有些terms的p值顺序不大一致，但因为都是显著的，也不影响结论。

是不是很激动，自己的转录组数据马上就能用起来啦，如果你完全看不懂上面的图表和代码，那么你距离成文还差一个学习班！

■   ■   ■ 

全国巡讲约你

七月份我们不外出，只专注单细胞！

系统学习单细胞分析，报名生信技能树的线下培训，手慢无。

一年一度的生信技能树单细胞线下培训班火热招生