Annotation and Visualisation of RNA-seq results

本文将介绍RNA-seq差异分析结果可视化的实战代码。

原文地址：https://bioinformatics-core-shared-training.github.io/RNAseq-R/rna-seq-annotation-visualisation.nb.html

书接上文（用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

除了火山图，还有各种对差异分析结果可视化的图就不在此展开了，有兴趣的朋友可以到原网站上继续了解。

暂完。