Minimap2是知名比对工具BWA的开发者Li Heng新开发的比对工具，它能够快速的将DNA或者mRNA序列比对到参考基因组上，使用场景有下面几种：
- 将PacBio或OXford Nanopore的read和已有参考基因组（如人类）进行比对
- 寻找高错误率read(15%)之间的overlap
- 将PacBio Iso-Seq 或Nanopore cDNA或RNA序列比对到参考基因组
- 将illumina 单端或者双端序列比对到参考基因组
- 组装之间的比对
- 临近物种的全基因组比对
最近用Canu的不同参数装了几个版本的基因组，希望比较不同组装之间是否连续，要是在组装结果A里面的是ABCD的排布，在结果B里面却出现了倒置变成了ACBD，那么我就要怀疑人生了。
我仔细看了一下PAF的输出格式，发现输出结果非常友好，但是不可能直接用read.table或者data.table::fread读取，所以我就自己写了一个解析函数。同时为了提高我对R语言中grid系统的理解，我又写了一个专门的画图代码
代码见最后，100多行的样子，保存为paf.R, 使用方法为:
source(".paf.R")
df <- read.paf("./asm_vs_ref.paf")
plot_synteny(df)
两个临近物种的结果图如下
效果图1局部
同一物种不同组装的效果图
效果图2局部
我原本以为画这种图需要累加坐标，后来发现原来grid系统里面的视图(viewport)可以继续分为很多个图层，每个图层可有不同的坐标系统。
read.paf <- function(file, sep = "\t",
header = FALSE,
MQ = 10,
...){
data <- readLines(file)
dataSize <- length(data)
# initialize
qName <- vector("character", dataSize)
qLength <- vector("integer", dataSize)
qStart <- vector("integer", dataSize)
qEnd <- vector("integer", dataSize)
strand <- vector("character", dataSize)
tName <- vector("character", dataSize)
tLength <- vector("integer", dataSize)
tStart <- vector("integer", dataSize)
tEnd <- vector("integer", dataSize)
reMatch <- vector("integer", dataSize)
bLength <- vector("integer", dataSize)
mQuality<- vector("integer", dataSize)
i <- 1
j <- 0
for (i in seq(dataSize)){
items <- strsplit(data[i], split = sep)[[1]]
quality <- as.numeric(items[12]) #Mapping Quality 0-255
if (quality < MQ)
next
j <- j + 1
qName[j] <- items[1] # Query sequence name
qLength[j] <- as.integer(items[2])
qStart[j] <- as.integer(items[3]) + 1L
qEnd[j] <- as.integer(items[4]) + 1L # convert 0-based to 1-based
strand[j] <- items[5] # Relative strand: "+" or "-"
tName[j] <- items[6] # Target sequence name
tLength[j] <- as.integer(items[7])
tStart[j] <- as.integer(items[8]) + 1L
tEnd[j] <- as.integer(items[9]) + 1L
reMatch[j] <- as.integer(items[10]) #Number of residue matches
bLength[j] <- as.integer(items[11]) #Alignment block length
mQuality[j]<- as.integer(items[12]) #Alignment block length
}
pafDataframe <- data.frame(qName = qName[1:j], qStart = qStart[1:j], qEnd = qEnd[1:j],
tName = tName[1:j], tStart = tStart[1:j], tEnd = tEnd[1:j],
qLength = qLength[1:j], tLength = tLength[1:j],
strand = strand[1:j],
reMatch = reMatch[1:j], bLength = bLength[1:j],
mQuality = mQuality[1:j],
stringsAsFactors = FALSE)
return(pafDataframe)
}
plot_synteny <- function(df, contigs = 20,
lineSize = 3,
borderCol = "#5496ff",... ){
# select the top N contig
x <- df[,c("qName","qLength")]
x <- x[!duplicated(x$qName),]
x <- x[order(x$qLength, decreasing = TRUE),][1:contigs,]
y <- df[,c("tName","tLength")]
y <- y[!duplicated(y$tName),]
y <- y[order(y$tLength, decreasing = TRUE),][1:contigs,]
# make new page for ploting
grid::grid.newpage()
# allot the ratio of each contig
x_frac <- x$qLength / sum(x$qLength)
y_frac <- y$tLength / sum(y$tLength)
# draw the contig name
grid::pushViewport(grid::viewport(height = 0.7,
width = 0.7,
gp = grid::gpar(cex = 0.75),
name = "contigName"))
x_pos <- c(0, cumsum(x_frac)[1:(contigs-1)])
y_pos <- 1- c(0, cumsum(y_frac)[1:(contigs-1)])
for (i in seq.int(1,contigs)){
grid::grid.text(label = x[i,1],
x = grid::unit(x_pos[i], "npc"),
y = grid::unit(1, "npc"),
just = c("left","bottom"),
rot = 35
)
grid::grid.text(label = y[i, 1],
x = grid::unit(1, "npc"),
y = grid::unit(y_pos[i], "npc"),
just = c("left")
)
}
vplay <- grid::grid.layout(contigs, contigs,
widths = x_frac,
heights = y_frac)
grid::pushViewport(grid::viewport(layout = vplay,
name = "vplay"))
# line represent the synteny
for ( i in seq(1, contigs)){
for (j in seq(1, contigs)){
# get the name and length for subsetting
xName <- x[i,1]
xLength <- x[i,2]
yName <- y[j,1]
yLength <- y[j,2]
# push view port for plot Collinearity
grid::pushViewport(grid::viewport(layout.pos.col = i,
layout.pos.row = j,
xscale = c(1, xLength),
yscale = c(1, yLength),
name = paste0("pos",i,j)))
grid::grid.rect(gp=grid::gpar(col=borderCol))
# select the data
plot_df <- df[df$qName == xName & df$tName == yName,]
blocks <- nrow(plot_df)
#cat(sprintf("block size is %d\n", blocks))
if (blocks == 0) {
grid::upViewport()
next
}
# plot the line
for (k in seq(1, blocks)){
if (plot_df$strand[k] == "+"){
grid::grid.lines(x = c(plot_df$qStart[k], plot_df$qEnd[k]),
y = c(plot_df$tStart[k], plot_df$tEnd[k]),
gp=grid::gpar(lwd = lineSize),
default.units = "native")
} else {
grid::grid.lines(x = c(plot_df$qStart[k], plot_df$qEnd[k]),
y = c(plot_df$tEnd[k], plot_df$tStart[k]),
gp=grid::gpar(lwd = lineSize),
default.units = "native")
}
}
grid::upViewport()
}
}
}
