在2年之前我写过一篇教程介绍MUMmer软件的使用方法,可以通过如何使用MUMmer比对大片段序列阅读。
MUMmer作为一个比对工具,它的主要功能就是寻找两个序列的相似之处,至于如何展示结果,并不是它的主要目标。这篇文章将会介绍如何基于MUMMER的输出结果进行可视化。
首先是下载数据,我们用细菌的基因组作为案例
wget http://mummer.sourceforge.net/examples/data/H_pylori26695_Eslice.fasta
wget http://mummer.sourceforge.net/examples/data/H_pyloriJ99_Eslice.fasta
然后使用nucmer
进行比对
nucmer H_pylori26695_Eslice.fasta H_pyloriJ99_Eslice.fasta
只保留1对1的最佳联配
delta-filter -1 out.delta > filter.delta
输出容易展示的信息
show-coords -T -q -H filter.delta > coord.txt
上面这几个程序的参数需要根据具体需求进行修改,并不是固定用法
后续就是在R语言里面继续绘图。 先加载数据并根据数据格式命名列
df <- read.table("~/coord.txt", sep = "\t")
colnames(df) <- c("ref_start", "ref_end", "qry_start", "qry_end", "ref_len", "qry_len",
"identiy", "ref_tag","qry_tag" )
获取x和y轴的范围
x_range <- range(c(df$ref_start, df$ref_end))
y_range <- range(c(df$qry_start, df$qry_end))
新建一个画图设备,并设置好画图区域
plot.new()
plot.window(xlim = x_range,
ylim = y_range)
之后逐行绘制每个联配结果,用不同的颜色来展示倒置的情况
for( i in 1:nrow(df)){
if (df[i,3] < df[i,4]){
lines(x = df[i,1:2], y = df[i,3:4], col = "red")
} else{
lines(x = df[i,1:2], y = df[i,3:4], col = "blue")
}
}
最后加上x轴和y轴的标签
box()
axis(1, at = seq(0, x_range[2], 10000), labels = seq(0, x_range[2], 10000) / 10000)
axis(2, at = seq(0, y_range[2], 10000), labels = seq(0, y_range[2], 10000) / 10000)
结果如下
上面这种作图方式就是R语言的基础作图模式,从绘图思想上叫做画笔模式,这区别于ggplot2所代表的图形语法。
我还写过一篇「R绘图」minimap2的PAF文件如何进行可视化?,介绍如何用grid系统进行画图。