如何使用R語言包c(diǎn)irclize可視化展示blast雙序列比對結(jié)果

circlize這個包挺強(qiáng)大的，R語言里用來畫圈圖非常方便。 

今天這篇文章記錄用circlize這個包畫圈圖展示blast雙序列比對結(jié)果的代碼

植物線粒體基因組類的文章通常會分析細(xì)胞器基因組間基因轉(zhuǎn)移情況，基本的分析方法就是blast比對?？梢暬故究梢赃x擇用這個圈圖來做

首先是使用blast建庫比對

makeblastdb -in mt.fasta -dbtype nucl -out mt
blastn -query cp.fasta -db mt -outfmt 6 > output.txt

然后是準(zhǔn)備數(shù)據(jù)用來畫最外圈用來展示兩條序列的部分

df<-data.frame(chr=c(rep("chloroplast",2),rep("mitochondrial",2)),
               x=c(1,131478,1,444567),
               y=c(0,1,0,1))

> df
            chr      x y
1   chloroplast      1 0
2   chloroplast 131478 1
3 mitochondrial      1 0
4 mitochondrial 444567 1

然后是讀入blast的輸出結(jié)果

df1<-read.csv("output6.txt",stringsAsFactors = F,header=F,sep="\t")

作圖用到的代碼

library(circlize)
library(RColorBrewer)
library(ComplexHeatmap)
col<-RColorBrewer::brewer.pal(6,"Paired")
circos.par("start.degree" = 130)
circos.initialize(factors = df$chr,x=df$x)
circos.trackPlotRegion(factors = df$chr,y=df$y,
                       panel.fun = function(x,y){
                         circos.axis()
                       },track.height = 0.1)
highlight.sector(sector.index = "chloroplast",col=col[1])
highlight.sector(sector.index = "mitochondrial",col=col[2])
circos.text(x=70000,y=0.5,
            labels = "chloroplast",
            sector.index = "chloroplast")
circos.text(x=220000,y=0.5,
            labels = "mitochondrial",
            sector.index = "mitochondrial",
            facing = "outside")
col_fun = colorRamp2(c(70,90,100),
                     c("green", "yellow", "red"))
for (i in 1:13){
  x<-sort(c(df1[i,8],df1[i,7]))
  y<-sort(c(df1[i,10],df1[i,9]))
  z<-df1[i,3]
  circos.link("chloroplast",x,"mitochondrial",y,
              col=add_transparency(col_fun(z)))
}
circos.clear()
lgd_links = Legend(at = c(70, 80, 90, 100), 
                   col_fun = col_fun, 
                   title_position = "topleft",
                   title = "identity(%)")
lgd_list_vertical = packLegend(lgd_links)

draw(lgd_list_vertical, x = unit(10, "mm"), 
     y = unit(10, "mm"), just = c("left", "bottom"))

調(diào)整整體的角度

circos.par("start.degree" = 130)

調(diào)整用來表示染色體的外圈粗細(xì)

circos.trackPlotRegion(factors = df$chr,y=df$y,
                       panel.fun = function(x,y){
                         circos.axis()
                       },track.height = 0.1)

畫圖的時候可以加一個track.height參數(shù)

調(diào)整外圈的刻度，現(xiàn)在展示的有點(diǎn)多，我想增大間隔，減少展示的數(shù)字，暫時不知道如何實(shí)現(xiàn)。