热图、韦恩图、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)
韦恩图
## 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值顺序不大一致，但因为都是显著的，也不影响结论。
是不是很激动，自己的转录组数据马上就能用起来啦，如果你完全看不懂上面的图表和代码，那么你距离成文还差一个学习班！
■ ■ ■
全国巡讲约你
第1-10站北上广深杭，西安，郑州， 吉林，武汉和成都（全部结束）
七月份我们不外出，只专注单细胞！
系统学习单细胞分析，报名生信技能树的线下培训，手慢无