Annotation and Visualisation of RNA-seq results
本文將介紹RNA-seq差異分析結果可視化的實戰(zhàn)代碼。
書接上文（用R也可以完成的RNA-Seq下游分析-3）
在差異分析后我們經(jīng)常要可視化我們的分析結果培慌，例如我們可以使用火山圖將差異表達的基因畫出來。
Annotaion
在本次分析中，我們都是采用基因的ID來進行的，但如果想便于我們自己查閱的話蔼紧，單純的ID是比較麻煩的拨拓，我們還得一個個ID找它對應的基因名。因此服猪，我們可以在差異分析結果的表格中添加上基因的symbol供填，genename等信息。
由于這次分析的物種是小鼠罢猪，我們采用以下的注釋包：
> library(org.Mm.eg.db)
# 同時載入之前分析的數(shù)據(jù)
> load("DE.Rdata")
讓我們查看一下該注釋包存有的數(shù)據(jù)
> 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"
在以上的數(shù)據(jù)中近她，我們需要的"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函數(shù)將新的注釋信息與原表格合并之后即可
> 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
相信大家對差異分析的可視化結果---火山圖，都不陌生泻骤。該圖呈現(xiàn)出了差異表達基因的基本分布情況漆羔，還可以通過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
除了火山圖鸟顺，還有各種對差異分析結果可視化的圖就不在此展開了惦蚊，有興趣的朋友可以到原網(wǎng)站上繼續(xù)了解。
暫完讯嫂。
