scripts/rerun_eland.py

   1 #!/usr/bin/env python
   2
   3 import logging
   4 from optparse import OptionParser
   5 import os
   6 import subprocess
   7 import sys
   8
   9 from gaworkflow.pipeline import gerald
  10 from gaworkflow.pipeline import runfolder
  11
  12 def make_query_filename(eland_obj, output_dir):
  13     query_name = '%s_%s_eland_query.txt'
  14     query_name %= (eland_obj.sample_name, eland_obj.lane_id)
  15
  16     query_pathname = os.path.join(output_dir, query_name)
  17
  18     if os.path.exists(query_pathname):
  19         logging.warn("overwriting %s" % (query_pathname,))
  20
  21     return query_pathname
  22
  23 def make_result_filename(eland_obj, output_dir):
  24     result_name = '%s_%s_eland_result.txt'
  25     result_name %= (eland_obj.sample_name, eland_obj.lane_id)
  26
  27     result_pathname = os.path.join(output_dir, result_name)
  28
  29     if os.path.exists(result_pathname):
  30         logging.warn("overwriting %s" % (result_pathname,))
  31
  32     return result_pathname
  33
  34 def extract_sequence(inpathname, query_pathname, length, dry_run=False):
  35     logging.info('extracting %d bases' %(length,))
  36     logging.info('extracting from %s' %(inpathname,))
  37     logging.info('extracting to %s' %(query_pathname,))
  38
  39     if not dry_run:
  40         try:
  41             instream = open(inpathname, 'r')
  42             outstream = open(query_pathname, 'w')
  43             gerald.extract_eland_sequence(instream, outstream, 0, length)
  44         finally:
  45             outstream.close()
  46             instream.close()
  47
  48 def run_eland(length, query_name, genome, result_name, multi=False, dry_run=False):
  49     cmdline = ['eland_%d' % (length,), query_name, genome, result_name]
  50     if multi:
  51         cmdline += ['--multi']
  52
  53     logging.info('running eland: ' + " ".join(cmdline))
  54     if not dry_run:
  55         return subprocess.Popen(cmdline)
  56     else:
  57         return None
  58
  59
  60 def rerun(gerald_dir, output_dir, length=25, dry_run=False):
  61     """
  62     look for eland files in gerald_dir and write a subset to output_dir
  63     """
  64     logging.info("Extracting %d bp from files in %s" % (length, gerald_dir))
  65     g = gerald.gerald(gerald_dir)
  66
  67     # this will only work if we're only missing the last dir in output_dir
  68     if not os.path.exists(output_dir):
  69         logging.info("Making %s" %(output_dir,))
  70         if not dry_run: os.mkdir(output_dir)
  71
  72     processes = []
  73     for lane_id, lane_param in g.lanes.items():
  74         eland = g.eland_results[lane_id]
  75
  76         inpathname = eland.pathname
  77         query_pathname = make_query_filename(eland, output_dir)
  78         result_pathname = make_result_filename(eland, output_dir)
  79
  80         extract_sequence(inpathname, query_pathname, length, dry_run=dry_run)
  81
  82         p = run_eland(length,
  83                       query_pathname,
  84                       lane_param.eland_genome,
  85                       result_pathname,
  86                       dry_run=dry_run)
  87         if p is not None:
  88             processes.append(p)
  89
  90     for p in processes:
  91         p.wait()
  92
  93 def make_parser():
  94     usage = '%prog: [options] runfolder'
  95
  96     parser = OptionParser(usage)
  97
  98     parser.add_option('--gerald',
  99                       help='specify location of GERALD directory',
 100                       default=None)
 101     parser.add_option('-o', '--output',
 102                       help='specify output location of files',
 103                       default=None)
 104     parser.add_option('-l', '--read-length', type='int',
 105                       help='specify new eland length',
 106                       dest='length',
 107                       default=25)
 108     parser.add_option('--dry-run', action='store_true',
 109                       help='only pretend to run',
 110                       default=False)
 111     parser.add_option('-v', '--verbose', action='store_true',
 112                       help='increase verbosity',
 113                       default=False)
 114
 115     return parser
 116
 117
 118 def main(cmdline=None):
 119     logging.basicConfig(level=logging.WARNING)
 120
 121     parser = make_parser()
 122     opts, args = parser.parse_args(cmdline)
 123
 124     if opts.length < 16 or opts.length > 32:
 125         parser.error("eland can only process reads in the range 16-32")
 126
 127     if len(args) > 1:
 128         parser.error("Can only process one runfolder directory")
 129     elif len(args) == 1:
 130         runs = runfolder.get_runs(args[0])
 131         if len(runs) != 1:
 132             parser.error("Not a runfolder")
 133         opts.gerald = runs[0].gerald.pathname
 134         if opts.output is None:
 135             opts.output = os.path.join(
 136                 runs[0].pathname,
 137                 'Data',
 138                 # pythons 0..n ==> elands 1..n+1
 139                 'C1-%d' % (opts.length+1,)
 140             )
 141
 142     elif opts.gerald is None:
 143         parser.error("need gerald directory")
 144
 145     if opts.output is None:
 146         parser.error("specify location for the new eland files")
 147
 148     if opts.verbose:
 149         root_logger = logging.getLogger()
 150         root_logger.setLevel(logging.INFO)
 151
 152     rerun(opts.gerald, opts.output, opts.length, dry_run=opts.dry_run)
 153
 154     return 0
 155
 156 if __name__ == "__main__":
 157     sys.exit(main(sys.argv[1:]))