Add support for extracting data out of Illumina's new RTA runfolder.
[htsworkflow.git] / htsworkflow / pipelines / runfolder.py
index fc2beeb425e3f1338874549e0f9b41502f392f87..14b7e07284bedcdb0507235cad9e7f4d021a6da6 100644 (file)
@@ -21,6 +21,7 @@ EUROPEAN_DATE_RE = "([0-9]{1,2}-[0-9]{1,2}-[0-9]{4,4})"
 VERSION_RE = "([0-9\.]+)"
 USER_RE = "([a-zA-Z0-9]+)"
 LANES_PER_FLOWCELL = 8
+LANE_LIST = range(1, LANES_PER_FLOWCELL+1)
 
 from htsworkflow.util.alphanum import alphanum
 from htsworkflow.util.ethelp import indent, flatten
@@ -136,6 +137,19 @@ class PipelineRun(object):
         tree = ElementTree.parse(filename).getroot()
         self.set_elements(tree)
 
+def load_pipeline_run_xml(pathname):
+    """
+    Load and instantiate a Pipeline run from a run xml file
+
+    :Parameters: 
+      - `pathname` : location of an run xml file
+
+    :Returns: initialized PipelineRun object
+    """
+    tree = ElementTree.parse(pathname).getroot()
+    run = PipelineRun(xml=tree)
+    return run
+
 def get_runs(runfolder):
     """
     Search through a run folder for all the various sub component runs
@@ -152,8 +166,10 @@ def get_runs(runfolder):
 
     def scan_post_image_analysis(runs, runfolder, image_analysis, pathname):
         logging.info("Looking for bustard directories in %s" % (pathname,))
-        bustard_glob = os.path.join(pathname, "Bustard*")
-        for bustard_pathname in glob(bustard_glob):
+        bustard_dirs = glob(os.path.join(pathname, "Bustard*"))
+        # RTA BaseCalls looks enough like Bustard.
+        bustard_dirs.extend(glob(os.path.join(pathname, "BaseCalls")))
+        for bustard_pathname in bustard_dirs:
             logging.info("Found bustard directory %s" % (bustard_pathname,))
             b = bustard.bustard(bustard_pathname)
             gerald_glob = os.path.join(bustard_pathname, 'GERALD*')
@@ -168,7 +184,7 @@ def get_runs(runfolder):
                     p.gerald = g
                     runs.append(p)
                 except IOError, e:
-                    print "Ignoring", str(e)
+                    logging.error("Ignoring " + str(e))
 
     datadir = os.path.join(runfolder, 'Data')
 
@@ -178,15 +194,99 @@ def get_runs(runfolder):
     for firecrest_pathname in glob(os.path.join(datadir,"*Firecrest*")):
         logging.info('Found firecrest in ' + datadir)
         image_analysis = firecrest.firecrest(firecrest_pathname)
-        scan_post_image_analysis(runs, runfolder, image_analysis, firecrest_pathname)
+        if image_analysis is None:
+            logging.warn(
+                "%s is an empty or invalid firecrest directory" % (firecrest_pathname,)
+            )
+        else:
+            scan_post_image_analysis(
+                runs, runfolder, image_analysis, firecrest_pathname
+            )
     # scan for IPAR directories
-    for ipar_pathname in glob(os.path.join(datadir,"IPAR_*")):
+    ipar_dirs = glob(os.path.join(datadir, "IPAR_*"))
+    # The Intensities directory from the RTA software looks a lot like IPAR
+    ipar_dirs.extend(glob(os.path.join(datadir, 'Intensities')))
+    for ipar_pathname in ipar_dirs:
         logging.info('Found ipar directories in ' + datadir)
         image_analysis = ipar.ipar(ipar_pathname)
-        scan_post_image_analysis(runs, runfolder, image_analysis, ipar_pathname)
+        if image_analysis is None:
+            logging.warn(
+                "%s is an empty or invalid IPAR directory" %(ipar_pathname,)
+            )
+        else:
+            scan_post_image_analysis(
+                runs, runfolder, image_analysis, ipar_pathname
+            )
 
     return runs
 
+def get_specific_run(gerald_dir):
+    """
+    Given a gerald directory, construct a PipelineRun out of its parents
+
+    Basically this allows specifying a particular run instead of the previous
+    get_runs which scans a runfolder for various combinations of
+    firecrest/ipar/bustard/gerald runs.
+    """
+    from htsworkflow.pipelines import firecrest
+    from htsworkflow.pipelines import ipar
+    from htsworkflow.pipelines import bustard
+    from htsworkflow.pipelines import gerald
+
+    gerald_dir = os.path.expanduser(gerald_dir)
+    bustard_dir = os.path.abspath(os.path.join(gerald_dir, '..'))
+    image_dir = os.path.abspath(os.path.join(gerald_dir, '..', '..'))
+
+    runfolder_dir = os.path.abspath(os.path.join(image_dir, '..','..'))
+   
+    logging.info('--- use-run detected options ---')
+    logging.info('runfolder: %s' % (runfolder_dir,))
+    logging.info('image_dir: %s' % (image_dir,))
+    logging.info('bustard_dir: %s' % (bustard_dir,))
+    logging.info('gerald_dir: %s' % (gerald_dir,))
+
+    # find our processed image dir
+    image_run = None
+    # split into parent, and leaf directory
+    # leaf directory should be an IPAR or firecrest directory
+    data_dir, short_image_dir = os.path.split(image_dir)
+    logging.info('data_dir: %s' % (data_dir,))
+    logging.info('short_iamge_dir: %s' %(short_image_dir,))
+
+    # guess which type of image processing directory we have by looking
+    # in the leaf directory name
+    if re.search('Firecrest', short_image_dir, re.IGNORECASE) is not None:
+        image_run = firecrest.firecrest(image_dir)
+    elif re.search('IPAR', short_image_dir, re.IGNORECASE) is not None:
+        image_run = ipar.ipar(image_dir)
+    elif re.search('Intensities', short_image_dir, re.IGNORECASE) is not None:
+        image_run = ipar.ipar(image_dir)
+
+    # if we din't find a run, report the error and return 
+    if image_run is None:
+        msg = '%s does not contain an image processing step' % (image_dir,)
+        logging.error(msg)
+        return None
+
+    # find our base calling
+    base_calling_run = bustard.bustard(bustard_dir)
+    if base_calling_run is None:
+        logging.error('%s does not contain a bustard run' % (bustard_dir,))
+        return None
+
+    # find alignments
+    gerald_run = gerald.gerald(gerald_dir)
+    if gerald_run is None:
+        logging.error('%s does not contain a gerald run' % (gerald_dir,))
+        return None
+
+    p = PipelineRun(runfolder_dir)
+    p.image_analysis = image_run
+    p.bustard = base_calling_run
+    p.gerald = gerald_run
+    
+    logging.info('Constructed PipelineRun from %s' % (gerald_dir,))
+    return p
 
 def extract_run_parameters(runs):
     """
@@ -195,7 +295,7 @@ def extract_run_parameters(runs):
     for run in runs:
       run.save()
 
-def summarize_mapped_reads(mapped_reads):
+def summarize_mapped_reads(genome_map, mapped_reads):
     """
     Summarize per chromosome reads into a genome count
     But handle spike-in/contamination symlinks seperately.
@@ -205,7 +305,7 @@ def summarize_mapped_reads(mapped_reads):
     genome = 'unknown'
     for k, v in mapped_reads.items():
         path, k = os.path.split(k)
-        if len(path) > 0:
+        if len(path) > 0 and not genome_map.has_key(path):
             genome = path
             genome_reads += v
         else:
@@ -216,28 +316,35 @@ def summarize_mapped_reads(mapped_reads):
 def summarize_lane(gerald, lane_id):
     report = []
     summary_results = gerald.summary.lane_results
-    eland_result = gerald.eland_results.results[lane_id]
-    report.append("Sample name %s" % (eland_result.sample_name))
-    report.append("Lane id %s" % (eland_result.lane_id,))
-    cluster = summary_results[eland_result.lane_id].cluster
-    report.append("Clusters %d +/- %d" % (cluster[0], cluster[1]))
-    report.append("Total Reads: %d" % (eland_result.reads))
-    mc = eland_result._match_codes
-    nm = mc['NM']
-    nm_percent = float(nm)/eland_result.reads  * 100
-    qc = mc['QC']
-    qc_percent = float(qc)/eland_result.reads * 100
-
-    report.append("No Match: %d (%2.2g %%)" % (nm, nm_percent))
-    report.append("QC Failed: %d (%2.2g %%)" % (qc, qc_percent))
-    report.append('Unique (0,1,2 mismatches) %d %d %d' % \
-                  (mc['U0'], mc['U1'], mc['U2']))
-    report.append('Repeat (0,1,2 mismatches) %d %d %d' % \
-                  (mc['R0'], mc['R1'], mc['R2']))
-    report.append("Mapped Reads")
-    mapped_reads = summarize_mapped_reads(eland_result.mapped_reads)
-    for name, counts in mapped_reads.items():
-      report.append("  %s: %d" % (name, counts))
+    for end in range(len(summary_results)):  
+      eland_result = gerald.eland_results.results[end][lane_id]
+      report.append("Sample name %s" % (eland_result.sample_name))
+      report.append("Lane id %s end %s" % (eland_result.lane_id, end))
+      cluster = summary_results[end][eland_result.lane_id].cluster
+      report.append("Clusters %d +/- %d" % (cluster[0], cluster[1]))
+      report.append("Total Reads: %d" % (eland_result.reads))
+
+      if hasattr(eland_result, 'match_codes'):
+          mc = eland_result.match_codes
+          nm = mc['NM']
+          nm_percent = float(nm)/eland_result.reads  * 100
+          qc = mc['QC']
+          qc_percent = float(qc)/eland_result.reads * 100
+
+          report.append("No Match: %d (%2.2g %%)" % (nm, nm_percent))
+          report.append("QC Failed: %d (%2.2g %%)" % (qc, qc_percent))
+          report.append('Unique (0,1,2 mismatches) %d %d %d' % \
+                        (mc['U0'], mc['U1'], mc['U2']))
+          report.append('Repeat (0,1,2 mismatches) %d %d %d' % \
+                        (mc['R0'], mc['R1'], mc['R2']))
+
+      if hasattr(eland_result, 'genome_map'):
+          report.append("Mapped Reads")
+          mapped_reads = summarize_mapped_reads(eland_result.genome_map, eland_result.mapped_reads)
+          for name, counts in mapped_reads.items():
+            report.append("  %s: %d" % (name, counts))
+
+      report.append('')
     return report
 
 def summary_report(runs):
@@ -249,7 +356,7 @@ def summary_report(runs):
         # print a run name?
         report.append('Summary for %s' % (run.name,))
        # sort the report
-       eland_keys = run.gerald.eland_results.results.keys()
+       eland_keys = run.gerald.eland_results.results[0].keys()
        eland_keys.sort(alphanum)
 
        for lane_id in eland_keys:
@@ -258,6 +365,14 @@ def summary_report(runs):
             report.append('')
         return os.linesep.join(report)
 
+def is_compressed(filename):
+    if os.path.splitext(filename)[1] == ".gz":
+        return True
+    elif os.path.splitext(filename)[1] == '.bz2':
+        return True
+    else:
+        return False
+
 def extract_results(runs, output_base_dir=None):
     if output_base_dir is None:
         output_base_dir = os.getcwd()
@@ -294,7 +409,15 @@ def extract_results(runs, output_base_dir=None):
 
       # tar score files
       score_files = []
-      for f in os.listdir(g.pathname):
+
+      # check for g.pathname/Temp a new feature of 1.1rc1
+      scores_path = g.pathname
+      scores_path_temp = os.path.join(scores_path, 'Temp')
+      if os.path.isdir(scores_path_temp):
+          scores_path = scores_path_temp
+
+      # hopefully we have a directory that contains s_*_score files
+      for f in os.listdir(scores_path):
           if re.match('.*_score.txt', f):
               score_files.append(f)
 
@@ -302,40 +425,85 @@ def extract_results(runs, output_base_dir=None):
       bzip_cmd = [ 'bzip2', '-9', '-c' ]
       tar_dest_name =os.path.join(cycle_dir, 'scores.tar.bz2')
       tar_dest = open(tar_dest_name, 'w')
-      logging.info("Compressing score files in %s" % (g.pathname,))
+      logging.info("Compressing score files from %s" % (scores_path,))
       logging.info("Running tar: " + " ".join(tar_cmd[:10]))
       logging.info("Running bzip2: " + " ".join(bzip_cmd))
       logging.info("Writing to %s" %(tar_dest_name))
 
-      tar = subprocess.Popen(tar_cmd, stdout=subprocess.PIPE, shell=False, cwd=g.pathname)
+      env = {'BZIP': '-9'}
+      tar = subprocess.Popen(tar_cmd, stdout=subprocess.PIPE, shell=False, env=env,
+                             cwd=scores_path)
       bzip = subprocess.Popen(bzip_cmd, stdin=tar.stdout, stdout=tar_dest)
       tar.wait()
 
       # copy & bzip eland files
-      for eland_lane in g.eland_results.values():
-          source_name = eland_lane.pathname
-          path, name = os.path.split(eland_lane.pathname)
-          dest_name = os.path.join(cycle_dir, name+'.bz2')
-
-          args = ['bzip2', '-9', '-c', source_name]
-          logging.info('Running: %s' % ( " ".join(args) ))
-          bzip_dest = open(dest_name, 'w')
-          bzip = subprocess.Popen(args, stdout=bzip_dest)
-          logging.info('Saving to %s' % (dest_name, ))
-          bzip.wait()
-
-def clean_runs(runs):
+      for lanes_dictionary in g.eland_results.results:
+          for eland_lane in lanes_dictionary.values():
+              source_name = eland_lane.pathname
+              path, name = os.path.split(eland_lane.pathname)
+              dest_name = os.path.join(cycle_dir, name)
+             logging.info("Saving eland file %s to %s" % \
+                          (source_name, dest_name))
+
+              if is_compressed(name):
+                logging.info('Already compressed, Saving to %s' % (dest_name, ))
+                shutil.copy(source_name, dest_name)
+              else:
+                # not compressed
+                dest_name += '.bz2'
+                args = ['bzip2', '-9', '-c', source_name]
+                logging.info('Running: %s' % ( " ".join(args) ))
+                bzip_dest = open(dest_name, 'w')
+                bzip = subprocess.Popen(args, stdout=bzip_dest)
+                logging.info('Saving to %s' % (dest_name, ))
+                bzip.wait()
+
+def rm_list(files, dry_run=True):
+    for f in files:
+        if os.path.exists(f):
+            logging.info('deleting %s' % (f,))
+            if not dry_run:
+                if os.path.isdir(f):
+                    shutil.rmtree(f)
+                else:
+                    os.unlink(f)
+        else:
+            logging.warn("%s doesn't exist."% (f,))
+
+def clean_runs(runs, dry_run=True):
     """
     Clean up run folders to optimize for compression.
     """
-    # TODO: implement this.
-    # rm RunLog*.xml
-    # rm pipeline_*.txt
-    # rm gclog.txt
-    # rm NetCopy.log
-    # rm nfn.log
-    # rm Images/L*
-    # cd Data/C1-*_Firecrest*
-    # make clean_intermediate
-
-    pass
+    if dry_run:
+        logging.info('In dry-run mode')
+
+    for run in runs:
+        logging.info('Cleaninging %s' % (run.pathname,))
+        # rm RunLog*.xml
+        runlogs = glob(os.path.join(run.pathname, 'RunLog*xml'))
+        rm_list(runlogs, dry_run)
+        # rm pipeline_*.txt
+        pipeline_logs = glob(os.path.join(run.pathname, 'pipeline*.txt'))
+        rm_list(pipeline_logs, dry_run)
+        # rm gclog.txt?
+        # rm NetCopy.log? Isn't this robocopy?
+        logs = glob(os.path.join(run.pathname, '*.log'))
+        rm_list(logs, dry_run)
+        # rm nfn.log?
+        # Calibration
+        calibration_dir = glob(os.path.join(run.pathname, 'Calibration_*'))
+        rm_list(calibration_dir, dry_run)
+        # rm Images/L*
+        logging.info("Cleaning images")
+        image_dirs = glob(os.path.join(run.pathname, 'Images', 'L*'))
+        rm_list(image_dirs, dry_run)
+        # cd Data/C1-*_Firecrest*
+        logging.info("Cleaning intermediate files")
+        # make clean_intermediate
+        if os.path.exists(os.path.join(run.image_analysis.pathname, 'Makefile')):
+            clean_process = subprocess.Popen(['make', 'clean_intermediate'], 
+                                             cwd=run.image_analysis.pathname,)
+            clean_process.wait()
+
+
+