test/jaccard.py

   1 #!/usr/bin/env python
   2 # Time-stamp: <2024-02-12 16:46:59 Tao Liu>
   3
   4 import sys
   5
   6 from MACS3.Signal.Region import *
   7
   8
   9 # ------------------------------------
  10 # Main function
  11 # ------------------------------------
  12 def main():
  13     if len(sys.argv) < 3:
  14         sys.stderr.write("Calculate Jaccard Index of two peak files and return the index.\n\nneed 2 paras: %s <peak1> <peak2>\n" % sys.argv[0])
  15         sys.exit(1)
  16
  17     #options.info("#  Read peak file 1...")
  18     peak1 = open( sys.argv[1] )
  19     regions1 = Regions()
  20     for l in peak1:
  21         if l.startswith("track"):
  22             continue
  23         fs = l.rstrip().split()
  24         regions1.add_loc( fs[0].encode(), int(fs[1]), int(fs[2]) )
  25     regions1.sort()
  26
  27     #options.info("#  Read peak file 2...")
  28     peak2 = open( sys.argv[2] )
  29     regions2 = Regions()
  30     for l in peak2:
  31         if l.startswith("track"):
  32             continue
  33         fs = l.rstrip().split()
  34         regions2.add_loc( fs[0].encode(), int(fs[1]), int(fs[2]) )
  35     regions2.sort()
  36
  37     tl1 = regions1.total_length()
  38     tl2 = regions2.total_length()
  39     inter = regions1.intersect(regions2)
  40     tl_inter = inter.total_length()
  41
  42     if (tl1 + tl2 - tl_inter) == 0:
  43         # this means both of the files are empty
  44         print ( 1.0 )
  45     else:
  46         jaccard_index = tl_inter / (tl1 + tl2 - tl_inter)
  47         print ( jaccard_index )
  48
  49 if __name__ == '__main__':
  50     main()