『代码记录』circlize篇

之前用TBtools做了一次圈图，颜值确实不低，circos早就折腾过，奈何没学过perl，用起来特别不顺手，顺便吐槽下circos的网站，跟着教程来，确实..一言难尽...
毕竟在接触过的这些语言中R还是最顺手的，早就听说circlize这个包挺不错，还没用过，今天花了一天时间，拿着cookbook现学现卖，基本100%复刻了之前那个TBtools的circos图。
不过折腾过circlize，circos最后得到一个结论，对于简单的圈图来说，TBtools真香！！
这里只放代码，有空了详细写下解读，不过我做的感觉已经够傻瓜的了..
还是国际惯例，先上图

20200911181244.jpg

###########################################
#       Prj: Heterosis ==> circRNA
#       Assignment: circos plot
#       Author: Shawn Wang
#       Date: Sep 11,2020
###########################################
library(circlize)
## 1st ==> need chr length to plot the chr track.
## import chrlen
chrlen <- read.table("~/02.Project/03.circRNA/CircRNAVersion2/02.data/04.other/Gh.Chr.Len.xls",
                     header = F,
                     sep = "\t",
                     stringsAsFactors = F)
head(chrlen)
## start from 0
chr = data.frame(Chr = chrlen$V1,
                 start = 0,
                 end = chrlen$V2)
## remove contig
chr = filter(chr,str_detect(chr$Chr,"^[A-D][0-9][0-9]"))
## 2nd ==> gene density plot. need chr bin with gene number
## import gene bin file
demo <- read.table("~/02.Project/03.circRNA/CircRNAVersion2/02.data/04.other/Gh.gene_density.txt",
                   header = F,
                   sep = "\t",
                   stringsAsFactors = F)
demo = data.frame(chr = demo$V1,
                  start = demo$V2 - demo[1,2],
                  end = demo$V2,
                  value = demo$V3)
## 3rd ==> circRNA and miRNA postion
## in this cycle, I'll plot miRNA as blue dot and circRNA as green dot in same layer.
# import miRNA postion file 
mi <- read.delim("~/02.Project/03.circRNA/CircRNAVersion2/02.data/05.circos/tmp/04.mi-pos.xls",
                 header = F,
                 sep = "\t",
                 stringsAsFactors = F)
# add values
mi.pos = data.frame(Chr = mi$V2,
                    start = mi$V3,
                    end = mi$V4,
                    value = 10)
## adjust start and end position
for (i in 1:nrow(mi.pos)) {
  start = mi.pos[i,2]
  end = mi.pos[i,3]
  if (start < end) {
    mi.pos[i,2] = start
    mi.pos[i,3] = end
  } else{
    mi.pos[i,2] = end
    mi.pos[i,3] = start
  }
}
## double check the start-end position
mi.pos[,2] < mi.pos[,3]
## import circ.pos
circ = read.delim("~/02.Project/03.circRNA/CircRNAVersion2/02.data/05.circos/02.workingDir/highlight/gh.circRNA.pos",
                  header = F,
                  sep = "\t",
                  stringsAsFactors = F)
circ.pos = data.frame(Chr = circ$V1,
                      start = circ$V2,
                      end = circ$V3,
                      value = 20)
head(circ.pos)
## merge two position file 
bed1 = circ.pos
bed2 = mi.pos
bed_list = list(bed1,bed2)
## links
## At2At
At2At = read.table("~/02.Project/03.circRNA/CircRNAVersion2/02.data/05.circos/02.workingDir/links/AtvsAt.link",
                   header = F,sep = "\t",stringsAsFactors = F)
fromAt1= At2At[,c(1:3)]
toAt1 = At2At[,c(4:6)]
## At2Dt
At2Dt = read.table("~/02.Project/03.circRNA/CircRNAVersion2/02.data/05.circos/02.workingDir/links/AtvsDt.link",
                   header = F,sep = "\t",stringsAsFactors = F)
fromAt2= At2Dt[,c(1:3)]
toDt2 = At2Dt[,c(4:6)]
## Dt2Dt
Dt2Dt = read.table("~/02.Project/03.circRNA/CircRNAVersion2/02.data/05.circos/02.workingDir/links/DtvsDt.link",
                   header = F,sep = "\t",stringsAsFactors = F)
fromDt3= Dt2Dt[,c(1:3)]
toDt3 = Dt2Dt[,c(4:6)]
## Dt2At
Dt2At = read.table("~/02.Project/03.circRNA/CircRNAVersion2/02.data/05.circos/02.workingDir/links/DtvsAt.link",
                   header = F,sep = "\t",stringsAsFactors = F)
fromDt4= Dt2At[,c(1:3)]
toAt4 = Dt2At[,c(4:6)]
## plot
circos.clear() ## cleaning last work
## track
circos.genomicInitialize(chr,major.by = 30000000)
## colored At and Dt subgenome with different color
chrcolor = rep(c("turquoise","#FF6347"),each = 13)
circos.genomicTrackPlotRegion(chr,ylim=c(0,0.1),track.height=0.03,bg.col=chrcolor,cell.padding=c(0.01, 0.5, 0.01, 0.5),track.margin=c(0, 0))
## gene density heatmap
min1 = quantile(demo$value)[1]
quant25 = quantile(demo$value)[2]
mid = quantile(demo$value)[3]
quant75 = quantile(demo$value)[4]
max = quantile(demo$value)[5]
col_fun = colorRamp2(c(min1,quant75, max), c("green","black","red"))
circos.genomicTrack(demo, numeric.column = 4,track.height = 0.06,ylim = c(0,1.2),
                    panel.fun = function(region,value,...){
                      circos.genomicRect(region, value, col = col_fun(value[[1]]), border = NA)
                    })
circos.genomicTrack(bed_list, track.height = 0.06,
                    panel.fun = function(region, value, ...) {
                      i = getI(...)
                      circos.genomicPoints(region, value, pch = 16, cex = 0.3, col = i, ...)
                    })
## links
# At2Dt
circos.genomicLink(fromAt2, toDt2, col = "#1E90FF", border = "#1E90FF", lwd = 0.15)
# Dt2At
circos.genomicLink(fromDt4, toAt4, col = "#FFD700", border = "#FFD700", lwd = 0.15)
# At2At
circos.genomicLink(fromAt1, toAt1, col = "#7FFFAA", border = "#7FFFAA",lwd = 0.12)
# Dt2Dt
circos.genomicLink(fromDt3, toDt3, col = "#FF6347", border = "#FF6347", lwd = 0.12)