Appendix. Plot with R

0) Prepare

本章我们介绍如何使用 R 进行数据可视化,我们将提供两种方案:一是在自己电脑使用 Rstudio 来画图(基于 .Rmd 文件),优点是使用方便,交互性强;一是在 Docker 容器中用命令行的方式来画图,优点是无需额外的安装和配置。

0a) 方案一: 在自己电脑上用 Rstudio 画图

  1. 在自己电脑安装软件和配置

    下载 相关资源中的 R, Rstudio 软件及 lulab-plot-master.zip

    1. 安装 R。
    2. 安装 RStudio。
    3. 下载并解压 lulab-plot-master.zip, 双击其中的 lulab-plot.Rproj
    4. 安装需要的package:

  2. 打开 .Rmd 文件

    用Rstudio打开all.Rmd文件, 即可阅读教程,并执行相关代码。

    tips: 如果你更喜欢每个文件仅包含一节的内容(一种 plot 类型),可以先打开 index.Rmd,安装需要的 packages,然后依次打开每一节对应的 .Rmd 文件(动画展示了第1、2小节对应的 1.box-plots.Rmd2.violin-plots.Rmd

0b) 方案二: 在 Docker 中使用 R 来画图

如果你在使用方案一时遇到了问题,也可以用我们提供的 Docker(里面已经预装好了 R 语言和需要的 packages)。

0b.1) 在容器中使用R

首先进入容器:

docker exec -it bioinfo_tsinghua bash

本章的操作均在 /home/test/plot/ 下进行:

cd /home/test/plot/

进入容器后,用以下命令进入 R 语言环境:

R

现在就可以运行 R 代码了,这里我们展示了计算 1, 2, ..., 10 的平均数。

mean(1:10)

在实际画图时,依次将下文给出的 R 代码复制到 Terminal 中运行。

运行完毕之后,用以下命令退出(按完 Enter 后,按 n 和 Enter),返回到容器:

q()

0b.2) load data, install packages, etc

  1. Prepare output directory

    dir.create('output')
    
  2. Load the data

    # Read the input files
    # “header=T” means that the data has a title, and sep="\t" is used as the separator
    data <-read.table("input/box_plots_mtcars.txt",header=T,sep="\t")
    # The function c(,,) means create the vector type data 
    df <- data[, c("mpg", "cyl", "wt")]
    
    df2 <-read.table("input/histogram_plots.txt",header=T,sep="\t")
    
    df3 <- read.table("input/volcano_plots.txt", header=T)
    
    df4 <- read.table("input/manhattan_plots_gwasResults.txt",header=T,sep="\t")
    
    df5 <-read.table("input/heatmaps.txt",header=T,sep="\t")
    
    # Covert data into matrix format
    # nrow(df5) and ncol(df5) return the number of rows and columns of matrix df5 respectively.
    dm <- data.matrix(df5[1:nrow(df5),2:ncol(df5)])
    
    # Get the row names
    row.names(dm) <- df5[,1]
    
    df6 <- read.table("input/ballon_plots_GO.txt", header=T, sep="\t")
    
    df7 <- read.table("input/box_plots_David_GO.txt",header=T,sep="\t")
    df7 <- df7[1:10,]
    
  3. Install R packages

    Docker 中已经装好所需要的 R 包,如果你是在自己电脑上运行,则需要安装 ggplot2, qqman, gplots, pheatmap, scales, reshape2, RColorBrewer 和 plotrix(使用 install.packages(), 如 install.packages('ggplot2'))。

  4. library R packages

    library(ggplot2)
    library(qqman)
    library(gplots)
    library(pheatmap)
    library(scales)
    library(reshape2)
    library(RColorBrewer)
    library(plyr)
    library(plotrix)
    

0b.3) Save the plot and view it

If you want to save the plot, please use either pdf() + dev.off() or ggsave().
The second one is specific for the ggplot2 package (i.e., if the code for plot starts with ggplot, then you can use the second one).

Let's see an example:

  1. pdf() + dev.off()

    # Begin to plot
    # Output as pdf
    pdf("output/1.1.Basic_boxplot.pdf", height = 3, width = 3)
    # Mapping the X and Y 
    # Components are constructed by using "+"
    ggplot(df, aes(x=cyl, y=mpg))+ 
    # draw the boxplot and fill it with gray
      geom_boxplot(fill="gray")+
    # Use the labs function to set the title and modify x and y
      labs(title="Plot of mpg per cyl",x="Cyl", y = "Mpg")+
    # Set the theme style
      theme_classic()
    
    # Save the plot
    dev.off()
    
  2. ggsave()

    # Begin to plot
    p <- ggplot(df, aes(x=cyl, y=mpg)) + 
      geom_boxplot(fill="gray")+
      labs(title="Plot of mpg per cyl",x="Cyl", y = "Mpg")+
      theme_classic()
    # Sava as pdf
    ggsave("output/1.1.Basic_boxplot.pdf", plot=p, height = 3,    width = 3)
    

If you want to view the produced file, you need to copy the file to /home/test/share, then open the bioinfo_tsinghua_share folder on the Desktop of host machine.

the following code is executed in Terminal, i.e., you need to quit R.

cp output/1.1.Basic_boxplot.pdf /home/test/share/

Here we only show one plot, in real use, you should replace the code for plot and change output file name to do more plots.

0b.4) 以上 3 步的动画

为了更清楚地展示方案二,我们制作了一个完整的动画:

动画中演示的是虚拟机中的操作步骤,我们首先用浏览器打开这一章,然后将一些代码复制到 Terminal 中去运行,最后查看生成的 plot。(On linux, you can use Ctrl + Insert to paste text in the clipboard to the terminal.)

正如你所看到的,方案二使用起来还是比较不方便的,所以如果没有特别的原因,我们还是推荐优先考虑方案一。

For the following sections, you can find all code in /home/test/plot/Rscripts/ or here (a file per chapter), and demo output in /home/test/plot/success/output/.

1) Box plots

  1. Basic box plot

    df$cyl <- as.factor(df$cyl)
    head(df)
    
    ###                    mpg cyl    wt
    ### Mazda RX4         21.0   6 2.620
    ### Mazda RX4 Wag     21.0   6 2.875
    ### Datsun 710        22.8   4 2.320
    ### Hornet 4 Drive    21.4   6 3.215
    ### Hornet Sportabout 18.7   8 3.440
    ### Valiant           18.1   6 3.460
    
    ggplot(df, aes(x=cyl, y=mpg)) + 
      geom_boxplot(fill="gray")+
      labs(title="Plot of mpg per cyl",x="Cyl", y = "Mpg")+
      theme_classic()
    

  2. Change continuous color by groups

    ggplot(df, aes(x=cyl, y=mpg, fill=cyl)) + 
      geom_boxplot()+
      labs(title="Plot of mpg per cyl",x="Cyl", y = "Mpg") +
      scale_fill_brewer(palette="Blues") + 
      theme_bw()
    

  3. Box plot for GO results

    df7$Term <- sapply(strsplit(as.vector(df7$Term),'~'),'[',2)
    head(df7)
    
    #          Category                                                         Term Count       X.      PValue
    #1 GOTERM_BP_DIRECT                               chemical synaptic transmission     6 4.651163 0.003873106
    #2 GOTERM_BP_DIRECT                                                cell motility     3 2.325581 0.007016402
    #3 GOTERM_BP_DIRECT negative regulation of intrinsic apoptotic signaling pathway     3 2.325581 0.011455205
    #4 GOTERM_BP_DIRECT                protein N-linked glycosylation via asparagine     3 2.325581 0.014940498
    #5 GOTERM_BP_DIRECT            positive regulation of androgen receptor activity     2 1.550388 0.017976476
    #6 GOTERM_BP_DIRECT                               photoreceptor cell maintenance     3 2.325581 0.024198625
    #                                                                                                                   Genes
    #1 ENSMUSG00000032360, ENSMUSG00000020882, ENSMUSG00000000766, ENSMUSG00000020745, ENSMUSG00000029763, ENSMUSG00000066392
    #2                                                             ENSMUSG00000022665, ENSMUSG00000043850, ENSMUSG00000031078
    #3                                                             ENSMUSG00000095567, ENSMUSG00000036199, ENSMUSG00000030421
    #4                                                             ENSMUSG00000031232, ENSMUSG00000028277, ENSMUSG00000024172
    #5                                                                                 ENSMUSG00000038722, ENSMUSG00000028964
    #6                                                             ENSMUSG00000037493, ENSMUSG00000043850, ENSMUSG00000020212
    #  List.Total Pop.Hits Pop.Total Fold.Enrichment Bonferroni Benjamini       FDR
    #1        110      172     18082        5.734249  0.8975036 0.8975036  5.554012
    #2        110       21     18082       23.483117  0.9839676 0.8733810  9.848665
    #3        110       27     18082       18.264646  0.9988443 0.8950571 15.604073
    #4        110       31     18082       15.907918  0.9998546 0.8901964 19.881092
    #5        110        3     18082      109.587879  0.9999763 0.8811197 23.441198
    #6        110       40     18082       12.328636  0.9999994 0.9089683 30.281607
    
    ggplot(df7) + geom_bar(stat="identity", width=0.6, aes(Term,Fold.Enrichment, fill=-1*log10(PValue)),colour="#1d2a33") + 
      coord_flip() +
      scale_fill_gradient(low="#e8f3f7",high="#236eba")+
      labs(fill=expression(-log10_Pvalue), x="GO Terms",y="foldEnrichment", title="GO Biological Process") +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5))  +
      theme(axis.title.x =element_text(size=16), 
            axis.title.y=element_text(size=14)) +
      theme(axis.text.y = element_text(size = 10,face="bold"),
            axis.text.x = element_text(size = 12,face="bold"))
    

    ggplot(df7) + geom_bar(stat="identity", width=0.6, aes(Term,Fold.Enrichment, fill=-1*log10(PValue)),colour="#1d2a33") + 
      coord_flip() +
      scale_fill_gradient(low="#feff2b",high="#fe0100")+
      labs(fill=expression(-log10_Pvalue), x="GO Terms",y="foldEnrichment", title="GO Biological Process") +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5))  +
      theme(axis.title.x =element_text(size=16), 
            axis.title.y=element_text(size=14)) +
      theme(axis.text.y = element_text(size = 10,face="bold"),
            axis.text.x = element_text(size = 12,face="bold"))
    

Reference: http://www.sthda.com/english/wiki/ggplot2-box-plot-quick-start-guide-r-software-and-data-visualization

2) Violin plots

  1. Basic violin plot

    df$cyl <- as.factor(df$cyl)
    head(df)
    
    ###                    mpg cyl    wt
    ### Mazda RX4         21.0   6 2.620
    ### Mazda RX4 Wag     21.0   6 2.875
    ### Datsun 710        22.8   4 2.320
    ### Hornet 4 Drive    21.4   6 3.215
    ### Hornet Sportabout 18.7   8 3.440
    ### Valiant           18.1   6 3.460
    
    ggplot(df, aes(x=cyl, y=mpg)) +
        geom_violin(trim=FALSE) +
        labs(title="Plot of mpg per cyl", x="Cyl", y = "Mpg")
    

  2. Add summary statistics on a violin plot

    • Add median and quartile

      ggplot(df, aes(x=cyl, y=mpg)) + 
        geom_violin(trim=FALSE) +
        labs(title="Plot of mpg per cyl", x="Cyl", y = "Mpg") +
        stat_summary(fun.y=mean, geom="point", shape=23, size=2, color="red")
      

      or

      ggplot(df, aes(x=cyl, y=mpg)) + 
        geom_violin(trim=FALSE) +
        labs(title="Plot of mpg per cyl", x="Cyl", y = "Mpg") +
        geom_boxplot(width=0.1)
      

    • Add mean and standard deviation

      ggplot(df, aes(x=cyl, y=mpg)) + 
        geom_violin(trim=FALSE) +
        labs(title="Plot of mpg per cyl", x="Cyl", y = "Mpg") +
        stat_summary(fun.data="mean_sdl", fun.args = list(mult = 1), geom="crossbar", width=0.1 )
      

      or

      ggplot(df, aes(x=cyl, y=mpg)) + 
        geom_violin(trim=FALSE) +
        labs(title="Plot of mpg per cyl", x="Cyl", y = "Mpg") +
        stat_summary(fun.data=mean_sdl, fun.args = list(mult = 1), geom="pointrange", color="red")
      

  3. Change violin plot fill colors

    ggplot(df, aes(x=cyl, y=mpg, fill=cyl)) + 
      geom_violin(trim=FALSE) +
      geom_boxplot(width=0.1, fill="white") +
      labs(title="Plot of mpg per cyl", x="Cyl", y = "Mpg") +
      scale_fill_brewer(palette="Blues") + 
      theme_classic()
    

Reference: http://www.sthda.com/english/wiki/ggplot2-violin-plot-quick-start-guide-r-software-and-data-visualization

3) Histogram plots

  1. Basic histogram plot

    head(df2)
    
    ###   sex weight
    ### 1   F     49
    ### 2   F     56
    ### 3   F     60
    ### 4   F     43
    ### 5   F     57
    ### 6   F     58
    
    ggplot(df2, aes(x=weight)) + geom_histogram(binwidth=1)
    

  2. Add mean line on a histogram plot

    ggplot(df2, aes(x=weight)) + 
      geom_histogram(binwidth=1, color="black", fill="white") +
      geom_vline(aes(xintercept=mean(weight)),color="black", linetype="dashed", size=0.5)
    

  3. Change histogram plot fill colors

    ##Use the plyr package to calculate the average weight of each group :
    mu <- ddply(df2, "sex", summarise, grp.mean=mean(weight))
    head(mu)
    
    ###   sex grp.mean
    ### 1   F    54.70
    ### 2   M    65.36
    
    ##draw the plot
    ggplot(df2, aes(x=weight, color=sex)) +
      geom_histogram(binwidth=1, fill="white", position="dodge")+
      geom_vline(data=mu, aes(xintercept=grp.mean, color=sex), linetype="dashed") +
      scale_color_brewer(palette="Paired") + 
      theme_classic()+
      theme(legend.position="top")
    

Reference: http://www.sthda.com/english/wiki/ggplot2-histogram-plot-quick-start-guide-r-software-and-data-visualization

4) Density plots

  1. Basic density

    head(df2)
    
    ###   sex weight
    ### 1   F     49
    ### 2   F     56
    ### 3   F     60
    ### 4   F     43
    ### 5   F     57
    ### 6   F     58
    
    ggplot(df2, aes(x=weight)) + 
      geom_density()
    

  2. Add mean line on a density plot

    ggplot(df2, aes(x=weight)) +
      geom_density() +
      geom_vline(aes(xintercept=mean(weight)), color="black", linetype="dashed", size=0.5)
    

  3. Change density plot fill colors

    ##Use the plyr package plyr to calculate the average weight of each group :
    mu <- ddply(df2, "sex", summarise, grp.mean=mean(weight))
    head(mu)
    
    ###   sex grp.mean
    ### 1   F    54.70
    ### 2   M    65.36
    

    draw the plot

    • Change fill colors

      ggplot(df2, aes(x=weight, fill=sex)) +
        geom_density(alpha=0.7)+
        geom_vline(data=mu, aes(xintercept=grp.mean, color=sex), linetype="dashed")+
        labs(title="Weight density curve",x="Weight(kg)", y = "Density") + 
        scale_color_brewer(palette="Paired") +
        scale_fill_brewer(palette="Blues") +
        theme_classic()
      

    • Change line colors

      ggplot(df2, aes(x=weight, color=sex)) +
        geom_density()+
        geom_vline(data=mu, aes(xintercept=grp.mean, color=sex), linetype="dashed")+
        labs(title="Weight density curve",x="Weight(kg)", y = "Density") + 
        scale_color_brewer(palette="Paired") +
        theme_classic()
      

    • Combine histogram and density plots

      ggplot(df2, aes(x=weight, color=sex, fill=sex)) + 
        geom_histogram(binwidth=1, aes(y=..density..), alpha=0.5, position="identity") +
        geom_density(alpha=.2) +
        labs(title="Weight density curve",x="Weight(kg)", y = "Density") + 
        scale_color_brewer(palette="Paired") +
        scale_fill_brewer(palette="Blues") +
        theme_classic()
      

Reference: http://www.sthda.com/english/wiki/ggplot2-density-plot-quick-start-guide-r-software-and-data-visualization

5) Dot plots

  1. Basic dot plots

    df$cyl <- as.factor(df$cyl)
    head(df)
    
    ###                    mpg cyl    wt
    ### Mazda RX4         21.0   6 2.620
    ### Mazda RX4 Wag     21.0   6 2.875
    ### Datsun 710        22.8   4 2.320
    ### Hornet 4 Drive    21.4   6 3.215
    ### Hornet Sportabout 18.7   8 3.440
    ### Valiant           18.1   6 3.460
    
    ggplot(df, aes(x=cyl, y=mpg)) + 
      geom_dotplot(binaxis='y', stackdir='center', binwidth=1)
    

  2. Add mean and standard deviation

    ggplot(df, aes(x=cyl, y=mpg)) + 
      geom_dotplot(binaxis='y', stackdir='center', binwidth=1) + 
      stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), geom="crossbar", width=0.5)
    

    or

    ggplot(df, aes(x=cyl, y=mpg)) + 
      geom_dotplot(binaxis='y', stackdir='center', binwidth=1) + 
      stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), geom="pointrange", color="red")
    

  3. Change dot colors

    ggplot(df, aes(x=cyl, y=mpg, fill=cyl, shape=cyl)) + 
      geom_dotplot(binaxis='y', stackdir='center', binwidth=1, dotsize=0.8) + 
      labs(title="Plot of mpg per cyl",x="Cyl", y = "Mpg") +
      #stat_summary(fun.data="mean_sdl", fun.args = list(mult=1), geom="crossbar", width=0.5) +
      scale_fill_brewer(palette="Blues") +
      #scale_color_brewer(palette="Blues") +
      theme_classic()
    

  4. Change dot colors, shapes and align types

    ggplot(df, aes(x=cyl, y=mpg, color=cyl, shape=cyl)) + 
      geom_jitter(position=position_jitter(0.1), cex=2)+
      labs(title="Plot of mpg per cyl",x="Cyl", y = "Mpg") + 
      scale_color_brewer(palette="Blues") + 
      theme_classic()
    

Reference: http://www.sthda.com/english/wiki/ggplot2-dot-plot-quick-start-guide-r-software-and-data-visualization

6) Scatter plots

  1. Basic scatter plots

    df$cyl <- as.factor(df$cyl)
    head(df)
    
    ###                    mpg cyl    wt
    ### Mazda RX4         21.0   6 2.620
    ### Mazda RX4 Wag     21.0   6 2.875
    ### Datsun 710        22.8   4 2.320
    ### Hornet 4 Drive    21.4   6 3.215
    ### Hornet Sportabout 18.7   8 3.440
    ### Valiant           18.1   6 3.460
    
    ggplot(df, aes(x=wt, y=mpg)) + 
      geom_point(size=1.5)
    

  2. Add regression lines and change the point colors, shapes and sizes

    ggplot(df, aes(x=wt, y=mpg, color=cyl, shape=cyl)) +
      geom_point(size=1.5) + 
      geom_smooth(method=lm, se=FALSE, fullrange=TRUE) +
      theme_classic()
    

Reference: http://www.sthda.com/english/wiki/ggplot2-scatter-plots-quick-start-guide-r-software-and-data-visualization

7) Volcano plots

head(df3)
###      Gene log2FoldChange    pvalue      padj
### 1    DOK6         0.5100 1.861e-08 0.0003053
### 2    TBX5        -2.1290 5.655e-08 0.0004191
### 3 SLC32A1         0.9003 7.664e-08 0.0004191
### 4  IFITM1        -1.6870 3.735e-06 0.0068090
### 5   NUP93         0.3659 3.373e-06 0.0068090
### 6 EMILIN2         1.5340 2.976e-06 0.0068090
df3$threshold <- as.factor(ifelse(df3$padj < 0.05 & abs(df3$log2FoldChange) >=1,ifelse(df3$log2FoldChange > 1 ,'Up','Down'),'Not'))
ggplot(data=df3, aes(x=log2FoldChange, y =-log10(padj), color=threshold,fill=threshold)) +
  scale_color_manual(values=c("blue", "grey","red"))+
  geom_point(size=1) +
  xlim(c(-3, 3)) +
  theme_bw(base_size = 12, base_family = "Times") +
  geom_vline(xintercept=c(-1,1),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="log2FoldChange",y="-log10 (adjusted p-value)",title="Volcano plot of DEG", face="bold")

8) Manhattan plots

head(df4)
###   SNP CHR BP         P
### 1 rs1   1  1 0.9148060
### 2 rs2   1  2 0.9370754
### 3 rs3   1  3 0.2861395
### 4 rs4   1  4 0.8304476
### 5 rs5   1  5 0.6417455
### 6 rs6   1  6 0.5190959
manhattan(df4, main = "GWAS results", ylim = c(0, 8),
          cex = 0.5, cex.axis=0.8, col=c("dodgerblue4","deepskyblue"),
          #suggestiveline = F, genomewideline = F, #remove the suggestive and genome-wide significance lines
          chrlabs = as.character(c(1:22)))

9) Heatmaps

  1. Draw the heatmap with the gplots package, heatmap.2() function

    head(dm)
    
    ###       Control1      Tumor2 Control3     Tumor4 Control5     Tumor1
    ### Gene1 3.646058 -0.98990248 2.210404 -0.2063050 2.859744  1.3304284
    ### Gene2 4.271172 -1.16217765 2.734119 -2.4782173 3.752013  0.0255639
    ### Gene3 3.530448  1.11451101 1.635485 -0.4241215 3.701427  1.2263312
    ### Gene4 3.061122 -1.18791027 4.331229  0.8733314 2.349352  0.4825479
    ### Gene5 1.956817  0.25431042 1.984438  1.2713845 1.685917  1.4554739
    ### Gene6 2.000919  0.06015972 4.480901  0.9780682 3.063475 -0.4222994
    ###       Control2     Tumor3 Control4     Tumor5
    ### Gene1 2.690376  0.6135943 2.470413  0.5158246
    ### Gene2 4.471795  1.6516242 2.735508 -0.5837784
    ### Gene3 3.588787 -0.6349656 1.999844  0.1417349
    ### Gene4 1.854433 -1.2237684 1.154377 -0.9301261
    ### Gene5 2.445830  0.3316909 2.715163  0.1866400
    ### Gene6 3.585366  1.0689000 2.563422  1.3465830
    
    ##to draw high expression value in red, we use colorRampPalette instead of redblue in heatmap.2
    ##colorRampPalette is a function in the RColorBrewer package
    cr <- colorRampPalette(c("blue","white","red"))
    heatmap.2(dm,
              scale="row", #scale the rows, scale each gene's expression value
              key=T, keysize=1.1, 
              cexCol=0.9,cexRow=0.8,
              col=cr(1000),
              ColSideColors=c(rep(c("blue","red"),5)),
              density.info="none",trace="none",
              #dendrogram='none', #if you want to remove dendrogram 
              Colv = T,Rowv = T) #clusters by both row and col
    

  2. Draw the heatmap with the pheatmap package, pheatmap function

    ##add column and row annotations
    annotation_col = data.frame(CellType = factor(rep(c("Control", "Tumor"), 5)), Time = 1:5)
    rownames(annotation_col) = colnames(dm)
    annotation_row = data.frame(GeneClass = factor(rep(c("Path1", "Path2", "Path3"), c(10, 4, 6))))
    rownames(annotation_row) = paste("Gene", 1:20, sep = "")
    ##set colors of each group
    ann_colors = list(Time = c("white", "springgreen4"), 
                      CellType = c(Control = "#7FBC41", Tumor = "#DE77AE"),
                      GeneClass = c(Path1 = "#807DBA", Path2 = "#9E9AC8", Path3 = "#BCBDDC"))
    ##draw the heatmap
    pheatmap(dm, 
             cutree_col = 2, cutree_row = 3, #break up the heatmap by clusters you define
             cluster_rows=TRUE, show_rownames=TRUE, cluster_cols=TRUE, #by default, pheatmap clusters by both row and col
             annotation_col = annotation_col, annotation_row = annotation_row,annotation_colors = ann_colors)
    

  3. Draw the heatmap with the ggplot2 package

    ##9.3.1.cluster by row and col
    ##cluster and re-order rows
    rowclust = hclust(dist(dm))
    reordered = dm[rowclust$order,]
    ##cluster and re-order columns
    colclust = hclust(dist(t(dm)))
    ##9.3.2.scale each row value in [0,1]
    dm.reordered = reordered[, colclust$order]
    dm.reordered=apply(dm.reordered,1,rescale) #rescale is a function in the scales package
    dm.reordered=t(dm.reordered) #transposed matrix
    ##9.3.3.save col and row names before changing the matrix format
    col_name=colnames(dm.reordered) 
    row_name=rownames(dm.reordered) 
    ##9.3.4.change data format for geom_title 
    colnames(dm.reordered)=1:ncol(dm.reordered)
    rownames(dm.reordered)=1:nrow(dm.reordered)
    dm.reordered=melt(dm.reordered) #melt is a function in the reshape2 package
    head(dm.reordered)
    ##9.3.5.draw the heatmap
    ggplot(dm.reordered, aes(Var2, Var1)) + 
      geom_tile(aes(fill = value), color = "white") + 
      scale_fill_gradient(low = "white", high = "steelblue") +
      theme_grey(base_size = 10) + 
      labs(x = "", y = "") + 
      scale_x_continuous(expand = c(0, 0),labels=col_name,breaks=1:length(col_name)) + 
      scale_y_continuous(expand = c(0, 0),labels=row_name,breaks=1:length(row_name))
    

10) Ballon plots

  1. basic ballon plots

    head(df6)
    
    ###                    Biological.process Fold.enrichment X.log10.Pvalue. col
    ### 1    Small molecule metabolic process             1.0              16   1
    ### 2   Single-organism catabolic process             1.5              12   1
    ### 3           Oxoacid metabolic process             2.0              23   1
    ### 4 Small molecule biosynthetic process             2.5               6   1
    ### 5   Carboxylic acid metabolic process             2.7              24   1
    ### 6      Organic acid metabolic process             2.7              25   1
    
    ggplot(df6, aes(x=Fold.enrichment, y=Biological.process)) +
      geom_point(aes(size = X.log10.Pvalue.)) +
      scale_x_continuous(limits=c(0,7),breaks=0:7) +
      scale_size(breaks=c(1,5,10,15,20,25)) +
      theme_light() +
      theme(panel.border=element_rect(fill='transparent', color='black', size=1),
            plot.title = element_text(color="black", size=14, hjust=0.5, face="bold", lineheight=1),
            axis.title.x = element_text(color="black", size=12, face="bold"),
            axis.title.y = element_text(color="black", size=12, vjust=1.5, face="bold"),
            axis.text.x = element_text(size=12,color="black",face="bold"),
            axis.text.y = element_text(size=12,color="black",face="bold"),
            legend.text = element_text(color="black", size=10, hjust=-2),
            legend.position="bottom") +
      labs(x="Fold Enrichment",y="Biological Process",size="-log10(Pvalue)", title="GO Enrichment",face="bold")
    

  2. change the dot colors

    ggplot(df6, aes(x=col, y=Biological.process,color=X.log10.Pvalue.)) +
      geom_point(aes(size = Fold.enrichment)) +
      scale_x_discrete(limits=c("1")) +
      scale_size(breaks=c(1,2,4,6)) +
      scale_color_gradient(low="#fcbba1", high="#a50f15") +
      theme_classic() +
      theme(panel.border=element_rect(fill='transparent', color='black', size=1),
            plot.title = element_text(color="black", size=14, hjust=0.5, face="bold", lineheight=1),
            axis.title.x = element_blank(),
            axis.title.y = element_text(color="black", size=12, face="bold"),
            axis.text.x = element_blank(),
            axis.ticks = element_blank(),
            axis.text.y = element_text(size=12,color="black",face="bold"),
            legend.text = element_text(color="black", size=10)) +
      labs(y="Biological Process",size="Fold Enrichment", color="-Log10(Pvalue)",title="GO Enrichment",face="bold")
    

11) Vennpie plots

The vennpie plot is the combination of a venn diagram and a pie chart.

##11.1.data input (number of reads mapped to each category)
total=100
rRNA=5
mtRNA=7
intergenic=48 
introns=12
exons=30
upstream=3
downstream=6
not_near_genes=40

rest=total-rRNA-mtRNA
genic=rest-intergenic
introns_and_exons=introns+exons-genic


##11.2 draw the plot
## parameter for pie chart
iniR=0.2 # initial radius
colors=list(NO='white',total='black',mtRNA='#e5f5e0',rRNA='#a1d99b',
            genic='#3182bd',intergenic='#fec44f',introns='#fc9272',
            exons='#9ecae1',upstream='#ffeda0',downstream='#fee0d2',
            not_near_genes='#d95f0e')

## from outer circle to inner circle
##0 circle: blank
pie(1, radius=iniR, init.angle=90, col=c('white'), border = NA, labels='')
##4 circle: show genic:exons and intergenic:downstream
floating.pie(0,0,
             c(exons, genic-exons+not_near_genes, downstream, mtRNA+rRNA+intergenic-not_near_genes-downstream),
             radius=5*iniR, 
             startpos=pi/2, 
             col=as.character(colors[c('exons','NO','downstream','NO')]),
             border=NA)
##3 circle: show genic:introns and intergenic:not_near_genes | upstream
floating.pie(0,0,
             c(genic-introns, introns, not_near_genes, intergenic-upstream-not_near_genes, upstream, mtRNA+rRNA),
             radius=4*iniR,
             startpos=pi/2, 
             col=as.character(colors[c('NO','introns','not_near_genes','NO','upstream','NO')]),
             border=NA)
##2 circle: divide the rest into genic and intergenic
floating.pie(0,0,
             c(genic, intergenic, mtRNA+rRNA),
             radius=3*iniR, 
             startpos=pi/2, 
             col=as.character(colors[c('genic','intergenic','NO')]),
             border=NA)
##1 circle: for rRNA+mtRNA+rest
floating.pie(0,0, 
             c(rest, rRNA,mtRNA), 
             radius=2*iniR, 
             startpos=pi/2, 
             col=as.character(colors[c('NO','rRNA','mtRNA')]), 
             border = NA)
legend(0, 6*iniR, gsub("_"," ",names(colors)[-1]), 
       col=as.character(colors[-1]), 
       pch=19, bty='n', ncol=2)

### or, in one column with reads count and %
##names=gsub("_"," ",names(colors)[-1])
##values = sapply(names(colors)[-1], get)
##percent=format(100*values/total, digits=2, trim=T)
##values = format(values, big.mark=",", scientific=FALSE, trim=T)
##cl=as.character(colors[-1])
##pchs=rep(19, length(cl)); pchs[1]=1;
##legend(0, 5*iniR, paste(names," (",values,", ", percent,"%)", sep=""), 
##       col=cl, pch=pchs,bty='n', ncol=1, cex=0.6)

Reference: http://onetipperday.sterding.com/2014/09/vennpier-combination-of-venn-diagram.html

12) Learn more

  1. Guide to Great Beautiful Graphics in R

    http://www.sthda.com/english/wiki/ggplot2-essentials

  2. Top 50 ggplot2 Visualizations - The Master List (With Full R Code)

    http://r-statistics.co/Top50-Ggplot2-Visualizations-MasterList-R-Code.html

  3. Color Scheme Suggestion

    http://colorbrewer2.org
    Rcolor.pdf

  4. Plots Gitbook of Xiaochen Xi

results matching ""

    No results matching ""