SDbioinfo-数据可视化1-Circos图及R实现

library(tidyverse)library(IRanges)library(Biostrings) library(openxlsx)library(dplyr)library(data.table)library(rtracklayer)library(ComplexHeatmap)library(circlize)library(ggplot2)

####Step1: to get chromosome info(chr names + length)####第一步：获得染色体长度和名称信息作为circos plot圆盘极坐标Toxo.fasta_all readDNAStringSet("/Users/sidaye/Documents/R/AP2_all_analysis/Input/Toxo_ME49v62/ToxoDB-62_TgondiiME49_Genome.fasta") #names(Toxo.fasta_all)#rule out contigs/scaffolds and keep the chromosomeToxo.fasta c(1:14),]#head(Toxo.fasta)

Chr.len unlist(lapply(strsplit(names(Toxo.fasta), ' | '), '[[',7))Chr.len unlist(lapply(strsplit(Chr.len, '='), '[[',2))Chr.name unlist(lapply(strsplit(names(Toxo.fasta), ' | '), '[[',1))Chr.name unlist(lapply(strsplit(Chr.name, '_'), '[[',2))
Genome.len frame(Chr.name = Chr.name,                         Length = Chr.len)head(Genome.len) ###Datafrome contains info

####Step2. Make reference datafrome for the corresponding relationship between name of chromosome and length of chromosome####第二步：把Genome.len作为refref 

####3. Prepare the dataframe(bed file) for transcript track for genomicrectangular plot to show the transcript of whole genome####第三步：得到需要画的最外层圆轮的转录本/CDS区域的坐标储存在bed file##use the gtf file instead of gff file gtf_genome "/Users/sidaye/Documents/R/AP2_all_analysis/Input/Toxo_ME49v62/ToxoDB-62_TgondiiME49.gtf", sep = '\t')transcript_gtf % dplyr::filter(V3 == 'transcript')CDS_gtf% dplyr::filter(V3 == 'CDS')##transform gtf file into bed filefunction_gtf_to_bed   x % dplyr::select(V1, V4, V5, V9, V3, V7)  x % dplyr::rename(V1=V1, V2=V4, V3=V5, V4=V9, V5=V3, V6=V7)  x$V5 '0'    return(x)}transcript_bed transcript_names ' '), '[[',2))transcript_names ';'),'[[',1))transcript_bed$V4 #write.table(transcript_bed, '/Users/sidaye/Documents/R/AP2_all_analysis/Output/transcript_loci_Toxo_v62.bed', col.names = F, quote = F, row.names = F, sep = '\t')CDS_bed CDS_names ' '), '[[',4))CDS_names ';'), '[[',1))CDS_bed$V4 original_chr_names ' | '), '[[',1))original_chr_names 1:14]print(original_chr_names) ####这里是获取染色体的名称,需要注意染色体名称，因为后面到画图时染色体名称对应极坐标的x-axis

####Step4. To change the chromosomes names####第四步：确认和修改染色体名称，以免后面画图造成错误function_change_chr_names  function(x){  for (i in 1:length(original_chr_names)){    x$V1  gsub(original_chr_names[i], Chr.name[i], x$V1)  }  return(x)}head(transcript_bed)

transcript_bed_new transcript_bed_new %dplyr::filter(grepl('chr', V1)) #filter rows that contain a certain string 'chr'#####可以明显看到这里chromosome name 改过来了，例如TGME49_chrII 改成了chrII， chrII会作为后面画图的每一个染色体的x-axis名称head(transcript_bed_new ) 

CDS_bed_new CDS_bed_new %dplyr::filter(grepl('chr', V1))

####5. Prepare the dataframe for the peaking time of all the AP2s and genomic loci of each AP2####5.第五步：这里需要准备dot plots（显示每一个AP2基因的转录本表达peak time）的表格AP2s_transcript_loci  read.csv("/Users/sidaye/Documents/R/AP2_all_analysis/Input/ME49_AP2_genomicposition.csv")AP2s_list_source_id  AP2s_transcript_loci$source_idAP2s_transcript_loci  transcript_bed_new[(transcript_bed_new$V4 %in% AP2s_list_source_id),]max(AP2s_transcript_loci$V3-AP2s_transcript_loci$V2)

#21232head(AP2s_transcript_loci) ######显示需要画的所有AP2基因的坐标轴

####6. Prepare the dataframe for y-axis of dots####6. 第六步：这里根据每一一副图希望描述的特征和对象不同来进行修改，也就是储存基因特征的表格，这个表格是最终画图的表格####主要包含信息是每个基因的表达时间（这个可以是用数字表示的任何特征，比如相关性等等）和每个基因的坐标AP2s_transcript_loci$V4  unlist(lapply(strsplit(AP2s_transcript_loci$V4, "-"), "[[", 1))AP2_loci  data.frame(AP2 = AP2s_transcript_loci$V4,                       start = AP2s_transcript_loci$V2,                       end = AP2s_transcript_loci$V3,                       chr = AP2s_transcript_loci$V1)#####AP2_loci2  AP2_loci[c(53, 18:21, 63:67, 32:33, 42:44, 54:62, 22:24, 25:31, 45:52, 7:17, 1:6, 34:41),]AP2_Peak  data.frame(chr = AP2s_transcript_loci$V4,                       start = AP2s_transcript_loci$V2,                       end = AP2s_transcript_loci$V3,                       time = rep(0, 67))AP2_Peak2AP2_Peak[c(53, 18:21, 63:67, 32:33, 42:44, 54:62, 22:24, 25:31, 45:52, 7:17, 1:6, 34:41),]#AP2_Peak2$time head(AP2_Peak2)

AP2_microarray  read.xlsx("/Users/sidaye/Documents/R/AP2_all_analysis/Input/microarray_cyclic_timing2.xlsx")AP2_microarray$GeneID gsub('-','_', AP2_microarray$GeneID)AP2_microarray  AP2_microarray[AP2_microarray$GeneID %in% AP2_Peak2$chr, ]colnames(AP2_microarray)[1]  "chr"#####since the order of AP2 in AP2_Peak2 has already been adjusted to corresponding to the order of chromosom, so it should be fixed in order to make nested zooming figureAP2_Peak_all  left_join(AP2_Peak2, AP2_microarray, by="chr")AP2_Peak_constitutive  AP2_Peak_all %>% dplyr::filter(constitutive == 1)AP2_Peak_constitutive_list  AP2_Peak_constitutive$chrAP2_Peak_cyclic  AP2_Peak_all %>% dplyr::filter(cyclic == 1)AP2_Peak_noexpressed  AP2_Peak_all[is.na(AP2_Peak_all$expressed),]AP2_Peak_noexpressed_list AP2_Peak_noexpressed$chr ###########for those genes is constitutive, even though it has a peak time, just to make it as a line across the all cell cycleAP2_Peak2$time  AP2_Peak_all$peak.time##########label all the constitutive genes as peaking time "1.11111"AP2_Peak2  AP2_Peak2 %>% dplyr::mutate(time=ifelse(chr %in% AP2_Peak_constitutive_list, 1.11111, time))######Because the cell cycle is cyclic, so those AP2s have peaking time at 0 or 6 should indicate at both 0 and 6 (C1->G)AP2_Peak_cyclic  AP2_Peak2[AP2_Peak2$chr%in%AP2_Peak_cyclic$chr,]AP2_Peak_GtoC  AP2_Peak_cyclic %>% dplyr::filter(time == 0)AP2_Peak_GtoC$time  rep(6,nrow(AP2_Peak_GtoC))AP2_Peak_CtoG  AP2_Peak_cyclic %>% dplyr::filter(time == 6)AP2_Peak_CtoG$time  rep(0,nrow(AP2_Peak_CtoG))AP2_Peak3  rbind(AP2_Peak2,AP2_Peak_GtoC,AP2_Peak_CtoG)head(AP2_Peak3) ####准备好的最终画图表格格式

####7.prepare the corresponding dataframe for nested zooming connection in the plot####第七步.因为我们需要花nested zoom in circos plot，需要做一个另外的table提供nested zoom in的坐标信息####前三列是最外层一轮的坐标，后三列是里面nested circos plot的坐标，因为nested circos plot外层和里层极坐标不一样####需要设置关系表格correspondance，来表明两套极坐标的关系correspondance data.frame(chr = AP2s_transcript_loci$V1,                            start = AP2s_transcript_loci$V2,                            end = AP2s_transcript_loci$V3,                            AP2 = AP2s_transcript_loci$V4,                            start.1 = AP2s_transcript_loci$V2,                            end.1 = AP2s_transcript_loci$V3)###调整顺序correspondance2  correspondance[c(53, 18:21, 63:67, 32:33, 42:44, 54:62, 22:24, 25:31, 45:52, 7:17, 1:6, 34:41),]head(correspondance2)

########################################################Starting to plot circos########################################################Starting to plot circos########################################################Starting to plot circos#clear all the previous circos plotcircos.clear()#set color for text in the circos plotcol_text  "black"#set up 12 random colors # chrIa       chrIb       chrII      chrIII       chrIV        chrV       chrVI     chrVIIa     chrVIIb     chrVIII       chrIX #"#AD2D0533" "#FFC78833" "#1018E833" "#C0670E33" "#6EFFC833" "#95387433" "#E4BF8D33" "#0992C633" "#7E30CB33" "#A1884333" "#3FC7F433" # chrX       chrXI      chrXII #"#34E28D33" "#363ED033" "#78ECDE33" #chr_bg_color = rand_color(14, transparency = 0.8)chr_bg_color =c("#AD2D0533" , "#FFC78833", "#1018E833", "#C0670E33", "#6EFFC833", "#95387433", "#E4BF8D33", "#0992C633", "#7E30CB33", "#A1884333","#3FC7F433", "#34E28D33", "#363ED033", "#78ECDE33")names(chr_bg_color) = ref$Chr.name

#Sectors in `f1()` and `f2()` should be in the same order, or else the connection lines may overlapf1 = function() {  circos.par("track.height"=0.8,gap.degree=2,start.degree = 90,clock.wise = T,             cell.padding=c(0,0,0,0))  #create circos object, 14=number of chromosome  circos.initialize(factors=ref$Chr.name,                              xlim=matrix(c(rep(0,14),ref$Length),ncol=2)  )  #cex，Specify the size of the title text with a numeric value of length 1    #######Track1  #set the chromosome name as track1, set the bg.col as white  #change font by cex  circos.track(ylim=c(0,1),panel.fun=function(x,y) {    Genome=CELL_META$sector.index    xlim=CELL_META$xlim    ylim=CELL_META$ylim    circos.text(mean(xlim),mean(ylim),Genome,cex=3,col=col_text,                facing="outside",niceFacing=TRUE,font =2)  },bg.col="white",bg.border=F,track.height=0.05)    #set gap between tracks  set_track_gap(mm_h(1)) # 0.5cm  #set_track_gap(mm_h(2)) # 2mm    #set track2 and add scale to the track  circos.track(ylim=c(0,1),panel.fun=function(x,y) {    Genome=CELL_META$sector.index    xlim=CELL_META$xlim    ylim=CELL_META$ylim    circos.text(mean(xlim),mean(ylim),Genome,cex=0.03,col="#808080",                facing="outside",niceFacing=TRUE)  },bg.col=chr_bg_color,bg.border=F,track.height=0.02)    #ref check the largest length of the chromosome  #max(ref$Length)  #7353086  #brk 0.0, 1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0)*10^6  #line width=lwd  #circos.track(track.index = get.current.track.index(), panel.fun = function(x, y) {  #circos.axis(h="top",major.at=brk,labels=round(brk/10^6,1),labels.cex=1.2,  #col=col_text,labels.col=col_text,lwd=0.8, labels.facing="clockwise")  #},bg.border=F)    set_track_gap(mm_h(0.5))  ##track3 shows the positions for all genes transcript  circos.genomicTrackPlotRegion(    transcript_bed_new, track.height = 0.05, stack = F, bg.border = NA, ylim = c(0,1),     panel.fun = function(region, value, ...) {      circos.genomicRect(region, value, col = 'black', border = NA, ...)    } )}circos.yaxis.scale as.numeric(c(0, 3, 4.7, 5, 6))#sector.index for constitutive genessector.index.constitutive in%AP2_Peak_constitutive_list) sector.index.noexpressed in%AP2_Peak_noexpressed_list) ###################To change all the GeneID in input dataframe into AP2 namesAP2_geneID2name "/Users/sidaye/Documents/R/AP2_all_analysis/Input/AP2_geneID2name.xlsx")colnames(AP2_geneID2name)[1] "chr"AP2_Peak_all2 by="chr")AP2_Peak_all2$Ap2Name2 'P2'), '[[',2))#gap.after is a vector means the gap between each neighboring sectors#Total track.height =1#Input dataframe has standardized format:chr, start, end, timef2 = function() {  #this track shows the position of each AP2s on the genome and the nested zooming track for the peaking time during CDC  circos.par(cell.padding = c(0, 0, 0, 0), gap.after = c(rep(1, nrow(AP2_loci)-1), 10))  #To create new sectors for the second part  circos.genomicInitialize(AP2_loci2, plotType = NULL)  circos.track(ylim=c(0,1),bg.col="white",bg.border=F,track.height=0.18)  AP2_Peak3_new   AP2_Peak3_new % dplyr::mutate(x.center=round((start+end)/2))  sector_indices get.all.sector.index()  for (i in 1:67){    circos.text(AP2_Peak3_new$x.center[i], 0.5, AP2_Peak_all2$Ap2Name2[i],                sector.index = sector_indices[i],                track.index = get.current.track.index(),                cex=1.8,                col=ifelse(i %in% sector.index.constitutive,"#d42729" ,ifelse(i %in% sector.index.noexpressed ,"grey16" ,"#0080ff")),                font = 2,                facing="reverse.clockwise",niceFacing=TRUE)}   set_track_gap(cm_h(0.8))  #set_track_gap(mm_h(1.2))  #track.height here change the height of dot block  circos.genomicTrack(AP2_Peak3, ylim = c(0, 6),track.height=0.4,                       #the panel.fun argument is used to specify a function that will be called for each sector of a circular plot                      panel.fun = function(region, value, ...) {                        #here, the sector.index is the name of each AP2 genes                        sector_indices get.all.sector.index()                        rec_colors "#f8c09b", "#e7db9f", "#bfc9a8", "#dea9f9")                        circos.rect(CELL_META$cell.xlim[1],0,CELL_META$cell.xlim[2],3, col=rec_colors[1], border = NA)                        circos.rect(CELL_META$cell.xlim[1],3,CELL_META$cell.xlim[2],4.7, col=rec_colors[2], border = NA)                        circos.rect(CELL_META$cell.xlim[1],4.7,CELL_META$cell.xlim[2],5, col=rec_colors[3], border = NA)                        circos.rect(CELL_META$cell.xlim[1],5,CELL_META$cell.xlim[2],6, col=rec_colors[4], border = NA)                        for(h in circos.yaxis.scale) {                          #CELL_META is an easy way to get meta data in the current cell                          #The cell.xlim element specifies the range of x-axis values that are visible within the cell.                          circos.lines(CELL_META$cell.xlim, c(h, h), lty = 4, col = "grey16")                        }                        circos.genomicPoints(region, value,                                              col = ifelse(value[[1]] == 0|value[[1]] == 6, "#0080ff",ifelse(value[[1]] ==1.11111, "white", "#0080ff")),                                             pch = 16,                                             cex = ifelse(value[[1]] == 0|value[[1]] == 6, 1.8,ifelse(value[[1]] ==1.11111, 0.01, 1.8))                        )                                              }, bg.col = "white", track.margin = c(0.02, 0))      #for (i in sector.index.constitutive){  #circos.lines(c(AP2_Peak3_new$x.center[i], AP2_Peak3_new$x.center[i]), CELL_META$cell.ylim,   #sector.index = sector_indices[i],  #track.index = get.current.track.index(),  #lty = 1,   #col ="#d42729")  #}  circos.yaxis(side = "left", at = circos.yaxis.scale,                sector.index = get.all.sector.index()[1], labels = FALSE)}#circos.nested shows the zooming connecting parts, the order of sector should be the same as the df "AP2_loci2" and "AP2_Peak3"circos.nested(f1, f2, correspondance2, connection_col = chr_bg_color[correspondance2[[1]]])

#15inches X 15inches#legend("right", inset=.02, title="Category",#     c("Constitutive","Not expressed","Cyclic"), fill=c('#d42729', 'grey16', '#377EB8'), horiz=F, cex=0.8)circos.info()

#rand_color : generate random colors in circlize#circos.par: Parameters for the circular layoutcircos.clear()