Annotation and Visualisation of RNA-seq results
本文将介绍RNA-seq差异分析结果可视化的实战代码。
书接上文（用R也可以完成的RNA-Seq下游分析-3）
在差异分析后我们经常要可视化我们的分析结果，例如我们可以使用火山图将差异表达的基因画出来。
Annotaion
在本次分析中，我们都是采用基因的ID来进行的，但如果想便于我们自己查阅的话，单纯的ID是比较麻烦的，我们还得一个个ID找它对应的基因名。因此，我们可以在差异分析结果的表格中添加上基因的symbol，genename等信息。
由于这次分析的物种是小鼠，我们采用以下的注释包：
> library(org.Mm.eg.db)
# 同时载入之前分析的数据
> load("DE.Rdata")
让我们查看一下该注释包存有的数据
> keytypes(org.Mm.eg.db)
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID"
[7] "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" "GO" "GOALL"
[13] "IPI" "MGI" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM"
[19] "PMID" "PROSITE" "REFSEQ" "SYMBOL" "UNIGENE" "UNIPROT"
在以上的数据中，我们需要的"ENTREZID" 、"GENENAME"、"SYMBOL" 都在其中。
接下来，我们只要对应已有结果中的ENTREZID提取相应的GENENAME和SYMBOL即可。
# 将差异分析结果存储到results变量当中
> results <- as.data.frame(topTags(lrt.BvsL,n = Inf))
> ann <- select(org.Mm.eg.db,keys=rownames(results),
+ columns=c("ENTREZID","SYMBOL","GENENAME"))
'select()' returned 1:1 mapping between keys and columns
> head(ann)
ENTREZID SYMBOL GENENAME
1 110308 Krt5 keratin 5
2 50916 Irx4 Iroquois homeobox 4
3 12293 Cacna2d1 calcium channel, voltage-dependent, alpha2/delta subunit 1
4 56069 Il17b interleukin 17B
5 24117 Wif1 Wnt inhibitory factor 1
6 12818 Col14a1 collagen, type XIV, alpha 1
利用cbind函数将新的注释信息与原表格合并之后即可
> results.annotated <- cbind(results, ann)
> head(results.annotated)
logFC logCPM LR PValue FDR ENTREZID SYMBOL
110308 -8.940579 10.264297 24.89789 6.044844e-07 0.0004377961 110308 Krt5
50916 -8.636503 5.749781 24.80037 6.358512e-07 0.0004377961 50916 Irx4
12293 -8.362247 6.794788 24.68526 6.749827e-07 0.0004377961 12293 Cacna2d1
56069 -8.419433 6.124377 24.41532 7.764861e-07 0.0004377961 56069 Il17b
24117 -9.290691 6.757163 24.32506 8.137331e-07 0.0004377961 24117 Wif1
12818 -8.216790 8.172247 24.24233 8.494462e-07 0.0004377961 12818 Col14a1
GENENAME
110308 keratin 5
50916 Iroquois homeobox 4
12293 calcium channel, voltage-dependent, alpha2/delta subunit 1
56069 interleukin 17B
24117 Wnt inhibitory factor 1
12818 collagen, type XIV, alpha 1
此外，除了基因名的信息之外，还可以利用有基因组位点信息的包，如TxDb.Mmusculus.UCSC.mm10.knownGene(小鼠的)和GenomicRanges来注释基因座信息。
Visualization
相信大家对差异分析的可视化结果---火山图，都不陌生。该图呈现出了差异表达基因的基本分布情况，还可以通过p值的设定，让我们清晰地观察显著差异表达基因。
先看看最简单的图
> signif <- -log10(results.annotated$FDR)
> plot(results.annotated$logFC,signif,pch=16)
naive volcano plot
当然，各位R学习者肯定不会止步于画了这个图就算了。让我们用ggplot画出更漂亮的图吧！
> library(ggplot2)
> library(ggtheme)
# 首先添加上threshold
> results.annotated$threshold <- as.factor(
+ ifelse(results.annotated$FDR < 0.05 & abs(results.annotated$logFC) >=log2(2),
+ ifelse(results.annotated$logFC > log2(2) ,'Up','Down'),'Not'))
> plt <- ggplot(data=results.annotated, aes(x=logFC, y =-log10(FDR), colour=threshold,fill=threshold)) +
+ scale_color_manual(values=c("blue", "grey","red"))+
+ geom_point(alpha=0.4, size=1.2) +
+ theme_bw(base_size = 12, base_family = "Times") +
+ geom_vline(xintercept=c(-0.5,0.5),lty=4,col="grey",lwd=0.6)+
+ geom_hline(yintercept = -log10(0.05),lty=4,col="grey",lwd=0.6)+
+ theme(legend.position="right",
+ panel.grid=element_blank(),
+ legend.title = element_blank(),
+ legend.text= element_text(face="bold", color="black",family = "Times", size=8),
+ plot.title = element_text(hjust = 0.5),
+ axis.text.x = element_text(face="bold", color="black", size=12),
+ axis.text.y = element_text(face="bold", color="black", size=12),
+ axis.title.x = element_text(face="bold", color="black", size=12),
+ axis.title.y = element_text(face="bold",color="black", size=12)) +
+ labs( x="log2 (Fold Change)",y="-log10 (p-value)")
> plt
beatuiful volcano plot
除了火山图，还有各种对差异分析结果可视化的图就不在此展开了，有兴趣的朋友可以到原网站上继续了解。
暂完。
