R script for Circos plot !

Abhi — Tue, 11 Jul 2023 01:41:03 -0500
#!/usr/bin/env Rscript
library(RCircos)

# usage: Rscript make_circos.r     

# parse args
args = commandArgs(trailingOnly=TRUE)
sv.file <- args[1]
sample.name <- args[2]
gene.label.file <- args[3]
cnv.file <- args[4]
out.file <- args[5]
# TMP <- Sys.getenv("TMP_DIR") 
# tmp.bed = paste0(TMP ,"/" , sample.name, "_bkpts.bed")
tmp.bed = paste0(sample.name, "_bkpts.bed")

# load prereq data
data(UCSC.HG19.Human.CytoBandIdeogram)

# set core parameters
chr.exclude <- NULL;
cyto.info <- UCSC.HG19.Human.CytoBandIdeogram;
tracks.inside <- 10;
tracks.outside <- 5;
RCircos.Set.Core.Components(cyto.info, chr.exclude, tracks.inside, tracks.outside);
rcircos.params <- RCircos.Get.Plot.Parameters();
rcircos.params$text.size <- 1
RCircos.Reset.Plot.Parameters(rcircos.params)
rcircos.cyto <- RCircos.Get.Plot.Ideogram();
rcircos.position <- RCircos.Get.Plot.Positions();
RCircos.List.Plot.Parameters()

link.data <- tryCatch(read.table(sv.file, sep = ',', stringsAsFactors = F, header = T), error=function(e) data.frame())
                    
if (nrow(link.data) != 0) {
  
  link.data <- transform(link.data,
                         chromStart = as.numeric(chromStart),
                         chromEnd = as.numeric(chromEnd),
                         chromStart.1 = as.numeric(chromStart.1),
                         chromEnd.1 = as.numeric(chromEnd.1))
  
  # write a bed file of all breakpoints to intersect with gene label table
  bkpts.1 <- link.data[c("Chromosome", "chromStart", "chromEnd")]
  bkpts.2 <- link.data[c("Chromosome.1", "chromStart.1", "chromEnd.1")]
  colnames(bkpts.2) <- colnames(bkpts.1)
  write.table(rbind(bkpts.1, bkpts.2), tmp.bed, sep = '\t', quote = F, col.names = F, row.names = F)
  
  # only keep labels that fall within an event
  print(paste0("bedtools intersect -wb -a ", tmp.bed, " -b ", gene.label.file))
  gene.labels <- system(paste0("bedtools intersect -wb -a ", tmp.bed, " -b ", gene.label.file), intern = T)
  gene.labels <- data.frame(do.call('rbind', strsplit(gene.labels, '\t', fixed=TRUE)), stringsAsFactors = F)
  if (nrow(gene.labels) > 0) {
    gene.labels <- gene.labels[,4:ncol(gene.labels)]
    
    # deduplicate labels
    gene.labels <- gene.labels[!duplicated(gene.labels),]
    colnames(gene.labels) <- c("Chromosome", "chromStart", "chromEnd", "gene")
    gene.labels <- transform(gene.labels,
                             chromStart = as.numeric(chromStart),
                             chromEnd = as.numeric(chromEnd))
  }
  
  # make the plot
  png(file=out.file, height=3000, width=3000, res = 500)
  RCircos.Set.Plot.Area()
  RCircos.Chromosome.Ideogram.Plot()
  track.num <- 2
  RCircos.Link.Plot(link.data, track.num, TRUE)
  title(sample.name, line=-1)
  
  # label the genes
  if (nrow(gene.labels) > 0) {
    name.col <- 4
    side <- "out"
    track.num <- 1
    RCircos.Gene.Connector.Plot(gene.labels, track.num, side);
    track.num <- 2
    RCircos.Gene.Name.Plot(gene.labels, name.col, track.num, side);
  }
    
  # remove intermediate file
  system(paste0("rm -f ", tmp.bed))
  
} else {
  # make empty plot
  png(file=out.file, height=3000, width=3000, res = 500)
  RCircos.Set.Plot.Area()
  RCircos.Chromosome.Ideogram.Plot()
  title(sample.name, line=-1)
}
                    
# parse cnv data
cnv = tryCatch(read.csv(cnv.file, stringsAsFactors = F), error=function(e) data.frame())
               
if (nrow(cnv) != 0) {
    colnames(cnv) <- c("Chromosome", "chromStart", "chromEnd", "cnv")
    cnv$Chromosome <- paste0('chr', cnv$Chromosome)
    cnv$GeneName <- "gene"
    cnv <- cnv[, c("Chromosome", "chromStart", "chromEnd", "GeneName", "cnv")]
}
                    
# add CNV heatmap track
if (nrow(cnv) != 0) {
  RCircos.Heatmap.Plot(cnv, data.col = 5, track.num = 1, side = "in")
}
                   
dev.off()

#-------- DATA FORMAT ------
chr1	11869	14412	DDX11L1
chr1	14363	29806	WASH7P
chr1	29554	31109	MIR1302-10
chr1	34554	36081	FAM138A
chr1	52473	54936	OR4G4P
chr1	62948	63887	OR4G11P
chr1	69091	70008	OR4F5
chr1	131025	134836	CICP27
chr1	134901	139379	AL627309.1
chr1	157784	157887	RNU6-1100P
chr1	227615	267253	AP006222.2
chr1	228292	228775	AP006222.1
chr1	317720	453948	RP4-669L17.10
chr1	326096	328112	RP4-669L17.8
chr1	329431	332236	CICP7
chr1	334126	334305	RP4-669L17.4
chr1	367640	368634	OR4F29
chr1	379105	379467	WBP1LP7
BOL: R script for Circos plot !

R script for Circos plot !
