MACS3/IO/OutputWriter.py

   1 # Time-stamp: <2020-12-01 12:57:23 Tao Liu>
   2
   3 """Module Description: Functions to write file
   4
   5 This code is free software; you can redistribute it and/or modify it
   6 under the terms of the BSD License (see the file LICENSE included with
   7 the distribution).
   8 """
   9
  10 # ------------------------------------
  11 # python modules
  12 # ------------------------------------
  13 import os
  14 import sys
  15 from array import array as pyarray
  16
  17 from MACS3.Utilities.Constants import MACS_VERSION
  18
  19 # ------------------------------------
  20 # constants
  21 # ------------------------------------
  22
  23 # ------------------------------------
  24 # Misc functions
  25 # ------------------------------------
  26 def zwig_write (trackI, subdir, fileprefix, d, log=None,space=10, single=False):
  27     """Write shifted tags information in wiggle file in a given
  28     step. Then compress it using 'gzip' program.
  29
  30     trackI: shifted tags from PeakDetect object
  31     subdir: directory where to put the wiggle file
  32     fileprefix: wiggle file prefix
  33     d     : d length
  34     log   : logging function, default is sys.stderr.write
  35     space : space to write tag number on spots, default 10
  36     """
  37     if not log:
  38         log = lambda x: sys.stderr.write(x+"\n")
  39     chrs = trackI.get_chr_names()
  40     os.makedirs (subdir)
  41     step = 10000000 + 2*d
  42
  43     if single:
  44         log("write to a wiggle file")
  45         f = os.path.join(subdir,fileprefix+"_all"+".wig")
  46         wigfhd = open(f,"w")
  47         wigfhd.write("track type=wiggle_0 name=\"%s_all\" description=\"Extended tag pileup from MACS version %s for every %d bp\"\n" % (fileprefix.replace('_afterfiting',''), MACS_VERSION, space)) # data type line
  48
  49     for chrom in chrs:
  50         if not single:
  51             f = os.path.join(subdir,fileprefix+"_"+chrom+".wig")
  52             log("write to "+f+" for chromosome "+chrom)
  53             wigfhd = open(f,"w")
  54             # suggested by dawe
  55             wigfhd.write("track type=wiggle_0 name=\"%s_%s\" description=\"Extended tag pileup from MACS version %s for every %d bp\"\n" % ( fileprefix.replace('_afterfiting',''), chrom, MACS_VERSION, space)) # data type line
  56         else:
  57             log("write data for chromosome "+chrom)
  58
  59         wigfhd.write("variableStep chrom=%s span=%d\n" % (chrom,space))
  60         tags = trackI.get_locations_by_chr(chrom)[0]
  61         l = len(tags)
  62         window_counts = pyarray('L',[0]*step)
  63         startp = -1*d
  64         endp   = startp+step
  65         index_tag = 0
  66
  67         while index_tag<l:
  68             s = tags[index_tag]-d/2     # start of tag
  69             e = s+d                     # end of tag
  70
  71             if e < endp:
  72                 # project tag to window_counts line
  73                 ps = s-startp # projection start
  74                 pe = ps+d     # projection end
  75                 for i in range(ps,pe):
  76                     window_counts[i] += 1
  77                 index_tag += 1
  78             else:
  79                 # write it to zwig file then reset parameters
  80                 # keep this tag for next window
  81                 for i in range(d,step-d,space):
  82                     if window_counts[i] == 0:
  83                         pass
  84                     else:
  85                         wigfhd.write("%d\t%d\n" % (i+startp+1,window_counts[i]))
  86                 # reset
  87                 window_counts_next = pyarray('L',[0]*step)
  88                 # copy d values from the tail of previous window to next window
  89                 for n,i in enumerate(range(step-2*d,step)): # debug
  90                     window_counts_next[n] = window_counts[i]
  91                 window_counts = window_counts_next
  92                 startp = endp - 2*d
  93                 endp = startp+step
  94         # last window
  95         for i in range(d,step-d,space):
  96             if window_counts[i] == 0:
  97                 pass
  98             else:
  99                 wigfhd.write("%d\t%d\n" % (i+startp+1,window_counts[i]))
 100         if not single:
 101             wigfhd.close()
 102             log("compress the wiggle file using gzip...")
 103             os.system("gzip "+f)
 104     if single:
 105         wigfhd.close()
 106         log("compress the wiggle file using gzip...")
 107         os.system("gzip "+f)
 108
 109
 110 def zbdg_write (trackI, subdir, fileprefix, d, log=None, single=False):
 111     """Write shifted tags information in wiggle file in a given
 112     step. Then compress it using 'gzip' program.
 113
 114     trackI: shifted tags from PeakDetect object
 115     subdir: directory where to put the wiggle file
 116     fileprefix: wiggle file prefix
 117     d     : d length
 118     log   : logging function, default is sys.stderr.write
 119     space : space to write tag number on spots, default 10
 120     """
 121     if not log:
 122         log = lambda x: sys.stderr.write(x+"\n")
 123     chrs = trackI.get_chr_names()
 124     os.makedirs (subdir)
 125     step = 10000000 + 2*d
 126
 127     if single:
 128         log("write to a bedGraph file")
 129         f = os.path.join(subdir,fileprefix+"_all"+".bdg")
 130         bdgfhd = open(f,"w")
 131         bdgfhd.write("track type=bedGraph name=\"%s_all\" description=\"Extended tag pileup from MACS version %s\"\n" % (fileprefix.replace('_afterfiting',''), MACS_VERSION)) # data type line
 132
 133     for chrom in chrs:
 134         if not single:
 135             f = os.path.join(subdir,fileprefix+"_"+chrom+".bdg")
 136             log("write to "+f+" for chromosome "+chrom)
 137             bdgfhd = open(f,"w")
 138             bdgfhd.write("track type=bedGraph name=\"%s_%s\" description=\"Extended tag pileup from MACS version %s\"\n" % (fileprefix.replace('_afterfiting',''), chrom, MACS_VERSION)) # data type line
 139         else:
 140             log("write data for chromosome "+chrom)
 141
 142         tags = trackI.get_locations_by_chr(chrom)[0]
 143         l = len(tags)
 144         window_counts = pyarray('L',[0]*step)
 145         startp = -1*d
 146         endp   = startp+step
 147         index_tag = 0
 148
 149         while index_tag<l:
 150             s = tags[index_tag]-d/2     # start of tag
 151             e = s+d                     # end of tag
 152
 153             if e < endp:
 154                 # project tag to window_counts line
 155                 ps = s-startp # projection start
 156                 pe = ps+d     # projection end
 157                 for i in range(ps,pe):
 158                     window_counts[i] += 1
 159                 index_tag += 1
 160             else:
 161                 # write it to zbdg file then reset parameters
 162                 # keep this tag for next window
 163                 prev = window_counts[d]
 164                 left = startp+d
 165                 right = left+1
 166                 for i in range(d+1,step-d):
 167                     if window_counts[i] == prev:
 168                         # same value, extend
 169                         right += 1
 170                     else:
 171                         # diff value, close
 172                         if prev != 0:
 173                             bdgfhd.write("%s\t%d\t%d\t%d\n" % (chrom,left,right,prev))
 174                         prev = window_counts[i]
 175                         left = right
 176                         right = left + 1
 177                 # last bin
 178                 if prev != 0:
 179                     bdgfhd.write("%s\t%d\t%d\t%d\n" % (chrom,left,right,prev))
 180
 181                 # reset
 182                 window_counts_next = pyarray('L',[0]*step)
 183                 # copy d values from the tail of previous window to next window
 184                 for n,i in enumerate(range(step-2*d,step)): # debug
 185                     window_counts_next[n] = window_counts[i]
 186                 window_counts = window_counts_next
 187                 startp = endp - 2*d
 188                 endp = startp+step
 189         # last window
 190         prev = window_counts[d]
 191         left = startp+d
 192         right = left+1
 193         for i in range(d+1,step-d):
 194             if window_counts[i] == prev:
 195                 # same value, exrend
 196                 right += 1
 197             else:
 198                 # diff value, close
 199                 if prev != 0:
 200                     bdgfhd.write("%s\t%d\t%d\t%d\n" % (chrom,left,right,prev))
 201                 prev = window_counts[i]
 202                 left = right
 203                 right = left + 1
 204         # last bin
 205         if prev != 0:
 206             bdgfhd.write("%s\t%d\t%d\t%d\n" % (chrom,left,right,prev))
 207
 208         if not single:
 209             bdgfhd.close()
 210             log("compress the bedGraph file using gzip...")
 211             os.system("gzip "+f)
 212     if single:
 213         bdgfhd.close()
 214         log("compress the bedGraph file using gzip...")
 215         os.system("gzip "+f)
 216
 217
 218 def model2r_script(model,filename,name):
 219     rfhd = open(filename,"w")
 220     p = model.plus_line
 221     m = model.minus_line
 222     ycorr = model.ycorr
 223     xcorr = model.xcorr
 224     alt_d = model.alternative_d
 225     #s = model.shifted_line
 226     d = model.d
 227     w = len(p)
 228     norm_p = [0]*w
 229     norm_m = [0]*w
 230     #norm_s = [0]*w
 231     sum_p = sum(p)
 232     sum_m = sum(m)
 233     #sum_s = sum(s)
 234     for i in range(w):
 235         norm_p[i] = float(p[i])*100/sum_p
 236         norm_m[i] = float(m[i])*100/sum_m
 237         #norm_s[i] = float(s[i])*100/sum_s
 238     rfhd.write("# R script for Peak Model\n")
 239     rfhd.write("#  -- generated by MACS\n")
 240
 241     rfhd.write("""p <- c(%s)
 242 m <- c(%s)
 243 ycorr <- c(%s)
 244 xcorr <- c(%s)
 245 altd  <- c(%s)
 246 x <- seq.int((length(p)-1)/2*-1,(length(p)-1)/2)
 247 pdf('%s_model.pdf',height=6,width=6)
 248 plot(x,p,type='l',col=c('red'),main='Peak Model',xlab='Distance to the middle',ylab='Percentage')
 249 lines(x,m,col=c('blue'))
 250 legend('topleft',c('forward tags','reverse tags'),lty=c(1,1,1),col=c('red','blue'))
 251 plot(xcorr,ycorr,type='l',col=c('black'),main='Cross-Correlation',xlab='Lag between + and - tags',ylab='Correlation')
 252 abline(v=altd,lty=2,col=c('red'))
 253 legend('topleft','alternative lag(s)',lty=2,col='red')
 254 legend('right','alt lag(s) : %s',bty='n')
 255 dev.off()
 256 """ % (','.join(map(str,norm_p)),
 257        ','.join(map(str,norm_m)),
 258        ','.join(map(str,ycorr)),
 259        ','.join(map(str,xcorr)),
 260        ', '.join(map(str,alt_d)),
 261        name,
 262        ','.join(map(str,alt_d))
 263        ))
 264     rfhd.close()
 265
 266 def diag_write (filename, diag_result):
 267     ofhd_diag = open(filename,"w")
 268     a = diag_result[0]
 269     l = len(a)-2
 270     s = [90-x*10 for x in range(l)]
 271     ofhd_diag.write("FC range\t# of Peaks\tcover by sampling %s\n" % ("%\t".join (map(str,s))+"%"))
 272     format = "%s\t%d"+"\t%.2f"*l+"\n"
 273     ofhd_diag.write( "".join( [format % tuple(x) for x in diag_result])  )
 274     ofhd_diag.close()
 275
 276 # ------------------------------------
 277 # Classes
 278 # ------------------------------------