美文网首页
批量运行弦图,封装函数

批量运行弦图,封装函数

作者: 千容安 | 来源:发表于2023-02-26 19:36 被阅读0次

    一共30个样本,要出弦图,先在jupyter里批量运行了PlotFancyVJUsage:

    def circos():
        import os
        for i in (1,7,35,189,217):
            for j in (1003,1501,1503,2001,2501,2503):
                x = "day" + str(i) + "-" + str(j)
                print(x)
                cmd_string = "java -jar ./vdjtools-1.2.1.jar PlotFancyVJUsage "+x+".txt "+ x
                print('x:{}'.format(cmd_string))
                print(os.popen(cmd_string).read())
                    
    circos()
    

    之前是在产生的错误文件(vj_pairing_plot.r)之里修改内容后一个个样本文件进行运行,想到这样太麻烦,就想写成一个函数,变化的是输入文件的名字
    输入的文件名称的规律:

    代码里输入文件的地方在:

    args <- c("day1-2501.fancyvj.wt.txt", "day1-2501.fancyvj.wt.pdf")
    temp <- read.table("day1-2501.fancyvj.wt.txt", sep="\t", comment="")
    

    本来这个循环我是写成这样的:

    但是要么报这个错

    要么报这个错

    因为手动输入文件名是可以运行的,所以觉得200G这个不存在……一直在找文件名哪里输的不对
    然后改成只运行一个循环,就成功了
    可能真的是需要那么多运行内存吧...

    circle_plot <- function(){
      for(i in c(1,7,35,189,217)){
          x = paste('day',i,'-2503',sep = "")
          y = '.fancyvj.wt.txt'
          z = '.fancyvj.wt.pdf'
          a = paste(x,y,sep = "")
          b = paste(x,z,sep = "")
          args <- c(a,b)
          #args<-commandArgs(TRUE)
        
          file_in  <- args[1]
          file_out <- args[2]
      
          require(circlize); require(RColorBrewer)
    
          # load data and preproc to fit formats
        
          temp <- read.table(a, sep="\t", comment="")
          n <- nrow(temp)
          m <- ncol(temp)
          rn = as.character(temp[2:n,1])
          cn = apply(temp[1,2:m], 2 , as.character)
          mat <- matrix(apply(temp[2:n, 2:m], 1:2, as.numeric), n - 1, m-1) * 100
          
          n <- nrow(temp)
          m <- ncol(temp)
          
          # Here columns and rows correspond to V and J segments respectively
          # Also replace possible duplicates (undef, '.', ...)
          
          duplicates <- intersect(rn, cn)
          
          rownames(mat) <- replace(rn, rn==duplicates, paste("V", duplicates, sep=""))
          colnames(mat) <- replace(cn, cn==duplicates, paste("J", duplicates, sep=""))
          
          # sort
          
          col_sum = apply(mat, 2, sum)
          row_sum = apply(mat, 1, sum)
          
          mat <- mat[order(row_sum), order(col_sum)]
          
          # equal number of characters for visualizaiton
          
          rn <- rownames(mat)
          cn <- colnames(mat)
          
          maxrn <- max(nchar(rn))
          maxcn <- max(nchar(cn))
          
          for(i in seq_len(length(rn))) {
                rn[i] <- paste(rn[i], paste(rep(" ", maxrn - nchar(rn[i])), collapse = ''))
          }
          
          for(i in seq_len(length(cn))) {
                cn[i] <- paste(cn[i], paste(rep(" ", maxcn - nchar(cn[i])), collapse = ''))
          }
          
          rownames(mat) <- rn
          colnames(mat) <- cn
          
          # viz using circlize
          
          if (grepl("\\.pdf$",file_out)){
             pdf(file_out)
          } else if (grepl("\\.png$",file_out)) {
             png(file_out, width     = 3.25,
                           height    = 3.25,
                           units     = "in",
                           res       = 1200,
                           pointsize = 4)
          } else {
             stop('Unknown plotting format')
          }
          
          circos.par(gap.degree = c(rep(1, nrow(mat)-1), 10, rep(1, ncol(mat)-1), 15), start.degree = 5)
          
          rcols <- rep(brewer.pal(12, "Paired"), nrow(mat)/12 + 1)[1:nrow(mat)]
          ccols <- rep(brewer.pal(12, "Paired"), ncol(mat)/12 + 1)[1:ncol(mat)]
          
          names(rcols) <- sort(rownames(mat))
          names(ccols) <- sort(colnames(mat))
          
          chordDiagram(mat, annotationTrack = "grid",
                       reduce = 0,
                       grid.col = c(rcols, ccols),
                       preAllocateTracks = list(track.height = 0.2), transparency = 0.5)
          
          circos.trackPlotRegion(track.index = 1, bg.border = NA,
                 panel.fun = function(x, y) {
                             sector.name = get.cell.meta.data("sector.index")
                             xlim = get.cell.meta.data("xlim")
                             ylim = get.cell.meta.data("ylim")
                             circos.text(mean(xlim), ylim[1], cex = 0.5, sector.name, facing = "clockwise", adj = c(0, 0.5))
                             }
                 )
          
          circos.clear()
          
          dev.off()
        }
    }
      
    circle_plot()
    
    

    相关文章

      网友评论

          本文标题:批量运行弦图,封装函数

          本文链接:https://www.haomeiwen.com/subject/qmycldtx.html