热图、韦恩图、GO富集分析图（有了转录组数据不知道该怎么写文章，看我就对了！）
上个月我们曾经学习过一篇细胞内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)
韦恩图
## 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-termsdf <- 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排序画barplotlibrary(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值顺序不大一致，但因为都是显著的，也不影响结论。
是不是很激动，自己的转录组数据马上就能用起来啦，如果你完全看不懂上面的图表和代码，那么你距离成文还差一个学习班！