本文展示如何使用ggplot2包绘制测序reads深度分布图，即用条形图来展示每个位点的测序深度，以及用饼图来展示检出病原的分布。
-
作图数据说明
-
传递文件1“VM00986.hist.data”，为测序reads深度，第一列碱基位置，第二列序列数，内容如下：
VM00986.hist.data
-
传递文件2“VM00986.txt”，指定坐标轴及图标题，内容如下：
VM00986.txt
-
传递文件3“VM00986.pie.data”，第一列为检出病原名称，第二列为序列数，第三列为百分比，内容如下：
VM00986.pie.data
-
绘图脚本说明
### 使用ggplot2绘制测序深度的条形图、饼图
library(ggplot2)
library(gridExtra) ## gridExtra包常用于控制图形的复杂布局和绘制文本表，如把几张图拼成一组图
## 传递参数
args <- commandArgs(TRUE)
## 读取数据：args[1]为文件"VM00986.hist“,用于画条形图
data1 <- read.table(args[1],sep="\t")
x1 <- data1[,1] ## data1第一列：碱基位置
y1 <- data1[,2] ## data2第二列：序列数，若数值差异较大可执行下一步的"截断坐标轴"，若不需要直接跳过下一步。
## 截断坐标轴：将序列数大于600的值全部替换成600
n<-nrow(data1) ## 计算行数
y1<-c(1:n)
for(i in 1:n){
if (data1[i,2]>600)
data1[i,2]<-600
y1[i]<-data1[i,2]
}
## 截断坐标轴
## 读取数据：args[2]为文件"VM00986.txt“,用于画条形图
data2 <- read.table(args[2],header=FALSE,sep="\t")
a <- data2[1,] ## 第一行：x轴标题
b <- data2[2,] ## 第二行：图标题
c <- data2[3,] ## 第三行：y轴标题
## ggplot绘制条形图
p1 <- ggplot(data1,aes(x=x1)) +
geom_bar(aes(y=y1),stat="identity") +
labs(x=a,title=b)+ylab(c) + ## x轴标题、图标题、y轴标题
theme(plot.title = element_text(hjust=0.5), ## 图标题居中
title=element_text(size=9), ## 图标题字体大小
axis.title.x=element_text(size=10), ## x轴标题字体大小
axis.title.y=element_text(size=10)) + ## y轴标题字体大小
scale_x_continuous(breaks = c(0,50,100,150,200,250,300,350,400,450)) ## x轴刻度标签
## 读取数据：args[3]为文件"VM00986.pie“,用于画饼图
data3 <- read.table(args[3],header=TRUE,sep="\t")
## 设置图例标题d：“microbe(n= reads)"
d <- paste("microbe (n=",sum(data3[,2])," reads)",sep="")
## 设置系列x标题：比如“Others (n=419, 0.89%)”,round()求百分比，保留2位小数
x2 <- paste(data3[,1]," (n=",data3[,2],", ",round(data3[,3]*100,2),"%",")",sep="")
y2 <- data3[,2] ## 第二列：y值
## 设置主题样式：theme_minimal()简约主题，去除xy轴字体和标题，去除边框线和网格线
blank_theme <- theme_minimal()+theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
panel.border = element_blank(),
panel.grid=element_blank())
## ggplot绘制饼图
p2 <- ggplot(data3,aes(x=factor(1),fill=factor(x2,order=TRUE,levels=x2))) +
geom_bar(stat="identity",aes(y=y2), colour="black", size = 0.1) + ## stat="identity"：条形的高度表示数据的值，而数据的值是由aes()函数的y参数决定的；colour=边框线的颜色；size=粗细
coord_polar(theta="y",start=1.5) + ## 准换成饼图，start=起始位置
labs(fill=d) + ## 图例标题d
blank_theme +
scale_fill_manual(values=c("grey80","grey70","grey60","grey50","grey40")) #填充颜色：根据系列x的数量，更改颜色参数的个数
## 将两张图拼成一张：排列为ncol列、nrow行
p3 <- grid.arrange(p1,p2,ncol=1)
-
运行脚本：操作环境cmd (传递参数)
> D: ## 切换路径
> cd D:\工作资料\科研项目\赣州人民医院\ICU-钟振洲\结果图\测序深度图
> Rscript.exe depth.R VM00986.hist.data VM00986.txt VM00986.pie.data ## 输出文件“Rplot.pdf”，结果如下：
hist-pie.pdf
结果发现饼图中后四种病原的比例过小，可重新画一个放大版的饼图：
-
使用ggplot2绘制饼图
-
作图数据说明
-
元数据文件“pie.data”，第一列为检出病原名称，第二列为序列数，第三列为百分比，内容如下：
pie.data
-
绘图脚本说明
### 使用ggplot2绘制饼图
library(ggplot2)
## 读取数据：饼图
data <- read.table("pie.data",header=TRUE,sep="\t")
## 设置图例的标题：比如“microbe(n= reads)",sep=""分隔符为空白
title <- paste("microbe (n=",sum(data[,2])," reads)",sep="")
## 设置系列名称x：比如“Others (n=419, 0.89%)”,round()保留2位小数
x <- paste(data[,1]," (n=",data[,2],", ",round(data[,3]*100,2),"%",")",sep="")
## y值为第2列
y <- data[,2]
## 设置主题样式：theme_minimal()简约主题，去除xy轴字体和标题，去除边框线和网格线
blank_theme <- theme_minimal()+theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
panel.border = element_blank(),
panel.grid=element_blank())
## ggplot绘制饼图
p2 <- ggplot(data,aes(x=factor(1),fill=factor(x,order=TRUE,levels=x)))+
geom_bar(stat="identity",aes(y=y),size = 0.1, colour="black")+ #条形图转换成饼图，colour边框颜色
coord_polar(theta="y",start=1.5)+ #弯曲横纵坐标,start起始位置90度
labs(fill=title)+ #图例标题
blank_theme+
scale_fill_manual(values=c("grey70","grey60","grey50","grey40")) #根据检出病原的数量，更改颜色参数的个数
p2 ## 将结果保存为pdf格式，如下图所示：
pie.pdf
-
图片拼接：AI软件
使用AI软件打开上述两个结果图的PDF文件，调整为如下格式即可：
VM00986.pdf
网友评论