import time
import re
-from htswdataprod.runfolder import \
+from htswdataprod.illumina.runfolder import \
ElementTree, \
VERSION_RE, \
EUROPEAN_STRPTIME
import re
import time
-from htswdataprod.runfolder import \
+from htswdataprod.illumina.runfolder import \
ElementTree, \
VERSION_RE, \
EUROPEAN_STRPTIME
import time
import types
-from htswdataprod.runfolder import \
+from htswdataprod.illumina.runfolder import \
ElementTree, \
EUROPEAN_STRPTIME, \
LANES_PER_FLOWCELL, \
--- /dev/null
+"""
+Utility functions to make bedfiles.
+"""
+import os
+import re
+
+# map eland_result.txt sense
+sense_map = { 'F': '+', 'R': '-'}
+sense_color = { 'F': '0,0,255', 'R': '255,255,0' }
+
+def write_bed_header(outstream, name, description):
+ """
+ Produce the headerline for a bedfile
+ """
+ # provide default track names
+ if name is None: name = "track"
+ if description is None: description = "eland result file"
+ bed_header = 'track name="%s" description="%s" visibility=4 itemRgb="ON"'
+ bed_header += os.linesep
+ outstream.write(bed_header % (name, description))
+
+def make_bed_from_eland_stream(instream, outstream, name, description, chromosome_prefix='chr'):
+ """
+ read an eland result file from instream and write a bedfile to outstream
+ """
+ # indexes into fields in eland_result.txt file
+ SEQ = 1
+ CHR = 6
+ START = 7
+ SENSE = 8
+
+ write_bed_header(outstream, name, description)
+
+ for line in instream:
+ fields = line.split()
+ # we need more than the CHR field, and it needs to match a chromosome
+ if len(fields) <= CHR or \
+ (chromosome_prefix is not None and \
+ fields[CHR][:3] != chromosome_prefix):
+ continue
+ start = fields[START]
+ stop = int(start) + len(fields[SEQ])
+ chromosome, extension = fields[CHR].split('.')
+ assert extension == "fa"
+ outstream.write('%s %s %d read 0 %s - - %s%s' % (
+ chromosome,
+ start,
+ stop,
+ sense_map[fields[SENSE]],
+ sense_color[fields[SENSE]],
+ os.linesep
+ ))
+
+
+def make_bed_from_multi_eland_stream(
+ instream,
+ outstream,
+ name,
+ description,
+ chr_prefix='chr',
+ max_reads=255
+ ):
+ """
+ read a multi eland stream and write a bedfile
+ """
+ write_bed_header(outstream, name, description)
+ parse_multi_eland(instream, outstream, chr_prefix, max_reads)
+
+def parse_multi_eland(instream, outstream, chr_prefix, max_reads=255):
+
+ loc_pattern = '(?P<fullloc>(?P<start>[0-9]+)(?P<dir>[FR])(?P<count>[0-9]+))'
+ other_pattern = '(?P<chr>[^:,]+)'
+ split_re = re.compile('(%s|%s)' % (loc_pattern, other_pattern))
+
+ for line in instream:
+ rec = line.split()
+ if len(rec) > 3:
+ # colony_id = rec[0]
+ seq = rec[1]
+ # number of matches for 0, 1, and 2 mismatches
+ # m0, m1, m2 = [int(x) for x in rec[2].split(':')]
+ compressed_reads = rec[3]
+ cur_chr = ""
+ reads = {0: [], 1: [], 2:[]}
+
+ for token in split_re.finditer(compressed_reads):
+ if token.group('chr') is not None:
+ cur_chr = token.group('chr')[:-3] # strip off .fa
+ elif token.group('fullloc') is not None:
+ matches = int(token.group('count'))
+ # only emit a bed line if
+ # our current chromosome starts with chromosome pattern
+ if chr_prefix is None or cur_chr.startswith(chr_prefix):
+ start = int(token.group('start'))
+ stop = start + len(seq)
+ orientation = token.group('dir')
+ strand = sense_map[orientation]
+ color = sense_color[orientation]
+ # build up list of reads for this record
+ reads[matches].append((cur_chr, start, stop, strand, color))
+
+ # report up to our max_read threshold reporting the fewer-mismatch
+ # matches first
+ reported_reads = 0
+ keys = [0,1,2]
+ for mismatch, read_list in ((k, reads[k]) for k in keys):
+ reported_reads += len(read_list)
+ if reported_reads <= max_reads:
+ for cur_chr, start, stop, strand, color in read_list:
+ reported_reads += 1
+ outstream.write('%s %d %d read 0 %s - - %s%s' % (
+ cur_chr,
+ start,
+ stop,
+ sense_map[orientation],
+ sense_color[orientation],
+ os.linesep
+ ))
+
+def make_description(database, flowcell_id, lane):
+ """
+ compute a bedfile name and description from the fctracker database
+ """
+ from gaworkflow.util.fctracker import fctracker
+
+ fc = fctracker(database)
+ cells = fc._get_flowcells("where flowcell_id='%s'" % (flowcell_id))
+ if len(cells) != 1:
+ raise RuntimeError("couldn't find flowcell id %s" % (flowcell_id))
+ lane = int(lane)
+ if lane < 1 or lane > 8:
+ raise RuntimeError("flowcells only have lanes 1-8")
+
+ name = "%s-%s" % (flowcell_id, lane)
+
+ cell_id, cell = cells.items()[0]
+ assert cell_id == flowcell_id
+
+ cell_library_id = cell['lane_%d_library_id' %(lane,)]
+ cell_library = cell['lane_%d_library' %(lane,)]
+ description = "%s-%s" % (cell_library['library_name'], cell_library_id)
+ return name, description
--- /dev/null
+"""
+Core information needed to inspect a runfolder.
+"""
+from glob import glob
+import logging
+import os
+import re
+import shutil
+import stat
+import subprocess
+import sys
+import time
+
+try:
+ from xml.etree import ElementTree
+except ImportError, e:
+ from elementtree import ElementTree
+
+EUROPEAN_STRPTIME = "%d-%m-%Y"
+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
+
+from htswcommon.util.alphanum import alphanum
+from htswcommon.util.ethelp import indent, flatten
+
+
+class PipelineRun(object):
+ """
+ Capture "interesting" information about a pipeline run
+ """
+ XML_VERSION = 1
+ PIPELINE_RUN = 'PipelineRun'
+ FLOWCELL_ID = 'FlowcellID'
+
+ def __init__(self, pathname=None, firecrest=None, bustard=None, gerald=None, xml=None):
+ if pathname is not None:
+ self.pathname = os.path.normpath(pathname)
+ else:
+ self.pathname = None
+ self._name = None
+ self._flowcell_id = None
+ self.firecrest = firecrest
+ self.bustard = bustard
+ self.gerald = gerald
+
+ if xml is not None:
+ self.set_elements(xml)
+
+ def _get_flowcell_id(self):
+ # extract flowcell ID
+ if self._flowcell_id is None:
+ config_dir = os.path.join(self.pathname, 'Config')
+ flowcell_id_path = os.path.join(config_dir, 'FlowcellId.xml')
+ if os.path.exists(flowcell_id_path):
+ flowcell_id_tree = ElementTree.parse(flowcell_id_path)
+ self._flowcell_id = flowcell_id_tree.findtext('Text')
+ else:
+ path_fields = self.pathname.split('_')
+ if len(path_fields) > 0:
+ # guessing last element of filename
+ flowcell_id = path_fields[-1]
+ else:
+ flowcell_id = 'unknown'
+
+ logging.warning(
+ "Flowcell id was not found, guessing %s" % (
+ flowcell_id))
+ self._flowcell_id = flowcell_id
+ return self._flowcell_id
+ flowcell_id = property(_get_flowcell_id)
+
+ def get_elements(self):
+ """
+ make one master xml file from all of our sub-components.
+ """
+ root = ElementTree.Element(PipelineRun.PIPELINE_RUN)
+ flowcell = ElementTree.SubElement(root, PipelineRun.FLOWCELL_ID)
+ flowcell.text = self.flowcell_id
+ root.append(self.firecrest.get_elements())
+ root.append(self.bustard.get_elements())
+ root.append(self.gerald.get_elements())
+ return root
+
+ def set_elements(self, tree):
+ # this file gets imported by all the others,
+ # so we need to hide the imports to avoid a cyclic imports
+ from htswdataprod.illumina import firecrest
+ from htswdataprod.illumina import bustard
+ from htswdataprod.illumina import gerald
+
+ tag = tree.tag.lower()
+ if tag != PipelineRun.PIPELINE_RUN.lower():
+ raise ValueError('Pipeline Run Expecting %s got %s' % (
+ PipelineRun.PIPELINE_RUN, tag))
+ for element in tree:
+ tag = element.tag.lower()
+ if tag == PipelineRun.FLOWCELL_ID.lower():
+ self._flowcell_id = element.text
+ #ok the xword.Xword.XWORD pattern for module.class.constant is lame
+ elif tag == firecrest.Firecrest.FIRECREST.lower():
+ self.firecrest = firecrest.Firecrest(xml=element)
+ elif tag == bustard.Bustard.BUSTARD.lower():
+ self.bustard = bustard.Bustard(xml=element)
+ elif tag == gerald.Gerald.GERALD.lower():
+ self.gerald = gerald.Gerald(xml=element)
+ else:
+ logging.warn('PipelineRun unrecognized tag %s' % (tag,))
+
+ def _get_run_name(self):
+ """
+ Given a run tuple, find the latest date and use that as our name
+ """
+ if self._name is None:
+ tmax = max(self.firecrest.time, self.bustard.time, self.gerald.time)
+ timestamp = time.strftime('%Y-%m-%d', time.localtime(tmax))
+ self._name = 'run_'+self.flowcell_id+"_"+timestamp+'.xml'
+ return self._name
+ name = property(_get_run_name)
+
+ def save(self, destdir=None):
+ if destdir is None:
+ destdir = ''
+ logging.info("Saving run report "+ self.name)
+ xml = self.get_elements()
+ indent(xml)
+ dest_pathname = os.path.join(destdir, self.name)
+ ElementTree.ElementTree(xml).write(dest_pathname)
+
+ def load(self, filename):
+ logging.info("Loading run report from " + filename)
+ tree = ElementTree.parse(filename).getroot()
+ self.set_elements(tree)
+
+def get_runs(runfolder):
+ """
+ Search through a run folder for all the various sub component runs
+ and then return a PipelineRun for each different combination.
+
+ For example if there are two different GERALD runs, this will
+ generate two different PipelineRun objects, that differ
+ in there gerald component.
+ """
+ from htswdataprod.illumina import firecrest
+ from htswdataprod.illumina import bustard
+ from htswdataprod.illumina import gerald
+
+ datadir = os.path.join(runfolder, 'Data')
+
+ logging.info('Searching for runs in ' + datadir)
+ runs = []
+ for firecrest_pathname in glob(os.path.join(datadir,"*Firecrest*")):
+ f = firecrest.firecrest(firecrest_pathname)
+ bustard_glob = os.path.join(firecrest_pathname, "Bustard*")
+ for bustard_pathname in glob(bustard_glob):
+ b = bustard.bustard(bustard_pathname)
+ gerald_glob = os.path.join(bustard_pathname, 'GERALD*')
+ for gerald_pathname in glob(gerald_glob):
+ try:
+ g = gerald.gerald(gerald_pathname)
+ runs.append(PipelineRun(runfolder, f, b, g))
+ except IOError, e:
+ print "Ignoring", str(e)
+ return runs
+
+
+def extract_run_parameters(runs):
+ """
+ Search through runfolder_path for various runs and grab their parameters
+ """
+ for run in runs:
+ run.save()
+
+def summarize_mapped_reads(mapped_reads):
+ """
+ Summarize per chromosome reads into a genome count
+ But handle spike-in/contamination symlinks seperately.
+ """
+ summarized_reads = {}
+ genome_reads = 0
+ genome = 'unknown'
+ for k, v in mapped_reads.items():
+ path, k = os.path.split(k)
+ if len(path) > 0:
+ genome = path
+ genome_reads += v
+ else:
+ summarized_reads[k] = summarized_reads.setdefault(k, 0) + v
+ summarized_reads[genome] = genome_reads
+ return summarized_reads
+
+def summary_report(runs):
+ """
+ Summarize cluster numbers and mapped read counts for a runfolder
+ """
+ report = []
+ for run in 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.sort(alphanum)
+
+ lane_results = run.gerald.summary.lane_results
+ for lane_id in eland_keys:
+ result = run.gerald.eland_results.results[lane_id]
+ report.append("Sample name %s" % (result.sample_name))
+ report.append("Lane id %s" % (result.lane_id,))
+ cluster = lane_results[result.lane_id].cluster
+ report.append("Clusters %d +/- %d" % (cluster[0], cluster[1]))
+ report.append("Total Reads: %d" % (result.reads))
+ mc = result._match_codes
+ nm = mc['NM']
+ nm_percent = float(nm)/result.reads * 100
+ qc = mc['QC']
+ qc_percent = float(qc)/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(result.mapped_reads)
+ for name, counts in mapped_reads.items():
+ report.append(" %s: %d" % (name, counts))
+ report.append('---')
+ report.append('')
+ return os.linesep.join(report)
+
+def extract_results(runs, output_base_dir=None):
+ if output_base_dir is None:
+ output_base_dir = os.getcwd()
+
+ for r in runs:
+ result_dir = os.path.join(output_base_dir, r.flowcell_id)
+ logging.info("Using %s as result directory" % (result_dir,))
+ if not os.path.exists(result_dir):
+ os.mkdir(result_dir)
+
+ # create cycle_dir
+ cycle = "C%d-%d" % (r.firecrest.start, r.firecrest.stop)
+ logging.info("Filling in %s" % (cycle,))
+ cycle_dir = os.path.join(result_dir, cycle)
+ if os.path.exists(cycle_dir):
+ logging.error("%s already exists, not overwriting" % (cycle_dir,))
+ continue
+ else:
+ os.mkdir(cycle_dir)
+
+ # copy stuff out of the main run
+ g = r.gerald
+
+ # save run file
+ r.save(cycle_dir)
+
+ # Copy Summary.htm
+ summary_path = os.path.join(r.gerald.pathname, 'Summary.htm')
+ if os.path.exists(summary_path):
+ logging.info('Copying %s to %s' % (summary_path, cycle_dir))
+ shutil.copy(summary_path, cycle_dir)
+ else:
+ logging.info('Summary file %s was not found' % (summary_path,))
+
+ # tar score files
+ score_files = []
+ for f in os.listdir(g.pathname):
+ if re.match('.*_score.txt', f):
+ score_files.append(f)
+
+ tar_cmd = ['/bin/tar', 'c'] + score_files
+ 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("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)
+ 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):
+ """
+ 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
--- /dev/null
+import unittest
+
+from StringIO import StringIO
+from htswdataprod import genome_mapper
+
+class testGenomeMapper(unittest.TestCase):
+ def test_construct_mapper(self):
+ genomes = {
+ 'Arabidopsis thaliana': {'v01212004': '/arabidopsis'},
+ 'Homo sapiens': {'hg18': '/hg18'},
+ 'Mus musculus': {'mm8': '/mm8',
+ 'mm9': '/mm9',
+ 'mm10': '/mm10'},
+ 'Phage': {'174': '/phi'},
+ }
+ genome_map = genome_mapper.constructMapperDict(genomes)
+
+ self.failUnlessEqual("%(Mus musculus|mm8)s" % (genome_map), "/mm8")
+ self.failUnlessEqual("%(Phage|174)s" % (genome_map), "/phi")
+ self.failUnlessEqual("%(Mus musculus)s" % (genome_map), "/mm10")
+ self.failUnlessEqual("%(Mus musculus|mm8)s" % (genome_map), "/mm8")
+ self.failUnlessEqual("%(Mus musculus|mm10)s" % (genome_map), "/mm10")
+
+ self.failUnlessEqual(len(genome_map.keys()), 6)
+ self.failUnlessEqual(len(genome_map.values()), 6)
+ self.failUnlessEqual(len(genome_map.items()), 6)
+
+
+def suite():
+ return unittest.makeSuite(testGenomeMapper,'test')
+
+if __name__ == "__main__":
+ unittest.main(defaultTest="suite")
--- /dev/null
+import os
+from StringIO import StringIO
+import unittest
+
+from htswdataprod.illumina import makebed
+
+class testMakeBed(unittest.TestCase):
+ def test_multi_1_0_0_limit_1(self):
+ instream = StringIO('>HWI-EAS229_26_209LVAAXX:7:3:112:383 TCAAATCTTATGCTANGAATCNCAAATTTTCT 1:0:0 mm9_chr13_random.fa:1240R0')
+ out = StringIO()
+
+ makebed.parse_multi_eland(instream, out, 'mm9_chr', 1)
+ self.failUnlessEqual(out.getvalue(), 'mm9_chr13_random 1240 1272 read 0 - - - 255,255,0\n')
+
+ def test_multi_1_0_0_limit_255(self):
+ instream = StringIO('>HWI-EAS229_26_209LVAAXX:7:3:112:383 TCAAATCTTATGCTANGAATCNCAAATTTTCT 1:0:0 mm9_chr13_random.fa:1240R0')
+ out = StringIO()
+
+ makebed.parse_multi_eland(instream, out, 'mm9_chr', 255)
+ self.failUnlessEqual(out.getvalue(), 'mm9_chr13_random 1240 1272 read 0 - - - 255,255,0\n')
+
+
+ def test_multi_2_0_0_limit_1(self):
+ instream = StringIO('>HWI-EAS229_26_209LVAAXX:7:3:104:586 GTTCTCGCATAAACTNACTCTNAATAGATTCA 2:0:0 mm9_chr4.fa:42995432F0,mm9_chrX.fa:101541458F0')
+ out = StringIO()
+
+ makebed.parse_multi_eland(instream, out, 'mm9_chr', 1)
+ self.failUnlessEqual(out.len, 0)
+
+ def test_multi_2_0_0_limit_255(self):
+ instream = StringIO('>HWI-EAS229_26_209LVAAXX:7:3:104:586 GTTCTCGCATAAACTNACTCTNAATAGATTCA 2:0:0 mm9_chr4.fa:42995432F0,mm9_chrX.fa:101541458F0')
+ out = StringIO()
+
+ makebed.parse_multi_eland(instream, out, 'mm9_chr', 255)
+ self.failUnlessEqual(out.len, 98)
+
+ def test_multi_0_2_0_limit_1(self):
+ instream = StringIO('>HWI-EAS229_26_209LVAAXX:7:3:115:495 TCTCCCTGAAAAATANAAGTGNTGTTGGTGAG 0:2:1 mm9_chr14.fa:104434729F2,mm9_chr16.fa:63263818R1,mm9_chr2.fa:52265438R1')
+ out = StringIO()
+
+ makebed.parse_multi_eland(instream, out, 'mm9_chr', 1)
+ print out.getvalue()
+ self.failUnlessEqual(out.len, 0)
+
+def suite():
+ return unittest.makeSuite(testMakeBed, 'test')
+
+if __name__ == "__main__":
+ unittest.main(defaultTest='suite')
+
+
--- /dev/null
+#!/usr/bin/env python
+
+from datetime import datetime, date
+import os
+import tempfile
+import shutil
+import unittest
+
+from htswdataprod.illumina import firecrest
+from htswdataprod.illumina import bustard
+from htswdataprod.illumina import gerald
+from htswdataprod.illumina import runfolder
+from htswdataprod.illumina.runfolder import ElementTree
+
+
+def make_flowcell_id(runfolder_dir, flowcell_id=None):
+ if flowcell_id is None:
+ flowcell_id = '207BTAAXY'
+
+ config = """<?xml version="1.0"?>
+<FlowcellId>
+ <Text>%s</Text>
+</FlowcellId>""" % (flowcell_id,)
+ config_dir = os.path.join(runfolder_dir, 'Config')
+
+ if not os.path.exists(config_dir):
+ os.mkdir(config_dir)
+ pathname = os.path.join(config_dir, 'FlowcellId.xml')
+ f = open(pathname,'w')
+ f.write(config)
+ f.close()
+
+def make_matrix(matrix_dir):
+ contents = """# Auto-generated frequency response matrix
+> A
+> C
+> G
+> T
+0.77 0.15 -0.04 -0.04
+0.76 1.02 -0.05 -0.06
+-0.10 -0.10 1.17 -0.03
+-0.13 -0.12 0.80 1.27
+"""
+ s_matrix = os.path.join(matrix_dir, 's_matrix.txt')
+ f = open(s_matrix, 'w')
+ f.write(contents)
+ f.close()
+
+def make_phasing_params(bustard_dir):
+ for lane in range(1,9):
+ pathname = os.path.join(bustard_dir, 'params%d.xml' % (lane))
+ f = open(pathname, 'w')
+ f.write("""<Parameters>
+ <Phasing>0.009900</Phasing>
+ <Prephasing>0.003500</Prephasing>
+</Parameters>
+""")
+ f.close()
+
+def make_gerald_config(gerald_dir):
+ config_xml = """<RunParameters>
+<ChipWideRunParameters>
+ <ANALYSIS>default</ANALYSIS>
+ <BAD_LANES></BAD_LANES>
+ <BAD_TILES></BAD_TILES>
+ <CONTAM_DIR></CONTAM_DIR>
+ <CONTAM_FILE></CONTAM_FILE>
+ <ELAND_GENOME>Need_to_specify_ELAND_genome_directory</ELAND_GENOME>
+ <ELAND_MULTIPLE_INSTANCES>8</ELAND_MULTIPLE_INSTANCES>
+ <ELAND_REPEAT></ELAND_REPEAT>
+ <EMAIL_DOMAIN>domain.com</EMAIL_DOMAIN>
+ <EMAIL_LIST>diane</EMAIL_LIST>
+ <EMAIL_SERVER>localhost:25</EMAIL_SERVER>
+ <EXPT_DIR>/home/diane/gec/080416_HWI-EAS229_0024_207BTAAXX/Data/C1-33_Firecrest1.8.28_19-04-2008_diane/Bustard1.8.28_19-04-2008_diane</EXPT_DIR>
+ <EXPT_DIR_ROOT>/home/diane/gec</EXPT_DIR_ROOT>
+ <FORCE>1</FORCE>
+ <GENOME_DIR>/home/diane/proj/SolexaPipeline-0.2.2.6/Goat/../Gerald/../../Genomes</GENOME_DIR>
+ <GENOME_FILE>Need_to_specify_genome_file_name</GENOME_FILE>
+ <HAMSTER_FLAG>genome</HAMSTER_FLAG>
+ <OUT_DIR>/home/diane/gec/080416_HWI-EAS229_0024_207BTAAXX/Data/C1-33_Firecrest1.8.28_19-04-2008_diane/Bustard1.8.28_19-04-2008_diane/GERALD_19-04-2008_diane</OUT_DIR>
+ <POST_RUN_COMMAND></POST_RUN_COMMAND>
+ <PRB_FILE_SUFFIX>_prb.txt</PRB_FILE_SUFFIX>
+ <PURE_BASES>12</PURE_BASES>
+ <QF_PARAMS>'((CHASTITY>=0.6))'</QF_PARAMS>
+ <QHG_FILE_SUFFIX>_qhg.txt</QHG_FILE_SUFFIX>
+ <QUALITY_FORMAT>--symbolic</QUALITY_FORMAT>
+ <READ_LENGTH>32</READ_LENGTH>
+ <SEQUENCE_FORMAT>--scarf</SEQUENCE_FORMAT>
+ <SEQ_FILE_SUFFIX>_seq.txt</SEQ_FILE_SUFFIX>
+ <SIG_FILE_SUFFIX_DEPHASED>_sig2.txt</SIG_FILE_SUFFIX_DEPHASED>
+ <SIG_FILE_SUFFIX_NOT_DEPHASED>_sig.txt</SIG_FILE_SUFFIX_NOT_DEPHASED>
+ <SOFTWARE_VERSION>@(#) Id: GERALD.pl,v 1.68.2.2 2007/06/13 11:08:49 km Exp</SOFTWARE_VERSION>
+ <TILE_REGEX>s_[1-8]_[0-9][0-9][0-9][0-9]</TILE_REGEX>
+ <TILE_ROOT>s</TILE_ROOT>
+ <TIME_STAMP>Sat Apr 19 19:08:30 2008</TIME_STAMP>
+ <TOOLS_DIR>/home/diane/proj/SolexaPipeline-0.2.2.6/Goat/../Gerald</TOOLS_DIR>
+ <USE_BASES>all</USE_BASES>
+ <WEB_DIR_ROOT>http://host.domain.com/yourshare/</WEB_DIR_ROOT>
+</ChipWideRunParameters>
+<LaneSpecificRunParameters>
+ <ANALYSIS>
+ <s_1>eland</s_1>
+ <s_2>eland</s_2>
+ <s_3>eland</s_3>
+ <s_4>eland</s_4>
+ <s_5>eland</s_5>
+ <s_6>eland</s_6>
+ <s_7>eland</s_7>
+ <s_8>eland</s_8>
+ </ANALYSIS>
+ <ELAND_GENOME>
+ <s_1>/g/dm3</s_1>
+ <s_2>/g/equcab1</s_2>
+ <s_3>/g/equcab1</s_3>
+ <s_4>/g/canfam2</s_4>
+ <s_5>/g/hg18</s_5>
+ <s_6>/g/hg18</s_6>
+ <s_7>/g/hg18</s_7>
+ <s_8>/g/hg18</s_8>
+ </ELAND_GENOME>
+ <READ_LENGTH>
+ <s_1>32</s_1>
+ <s_2>32</s_2>
+ <s_3>32</s_3>
+ <s_4>32</s_4>
+ <s_5>32</s_5>
+ <s_6>32</s_6>
+ <s_7>32</s_7>
+ <s_8>32</s_8>
+ </READ_LENGTH>
+ <USE_BASES>
+ <s_1>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_1>
+ <s_2>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_2>
+ <s_3>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_3>
+ <s_4>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_4>
+ <s_5>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_5>
+ <s_6>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_6>
+ <s_7>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_7>
+ <s_8>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_8>
+ </USE_BASES>
+</LaneSpecificRunParameters>
+</RunParameters>
+"""
+ pathname = os.path.join(gerald_dir, 'config.xml')
+ f = open(pathname,'w')
+ f.write(config_xml)
+ f.close()
+
+
+def make_summary_htm(gerald_dir):
+ summary_htm = """<!--RUN_TIME Mon Apr 21 11:52:25 2008 -->
+<!--SOFTWARE_VERSION @(#) $Id: jerboa.pl,v 1.31 2007/03/05 17:52:15 km Exp $-->
+<html>
+<body>
+
+<a name="Top"><h2><title>080416_HWI-EAS229_0024_207BTAAXX Summary</title></h2></a>
+<h1>Summary Information For Experiment 080416_HWI-EAS229_0024_207BTAAXX on Machine HWI-EAS229</h1>
+<h2><br></br>Chip Summary<br></br></h2>
+<table border="1" cellpadding="5">
+<tr><td>Machine</td><td>HWI-EAS229</td></tr>
+<tr><td>Run Folder</td><td>080416_HWI-EAS229_0024_207BTAAXX</td></tr>
+<tr><td>Chip ID</td><td>unknown</td></tr>
+</table>
+<h2><br></br>Lane Parameter Summary<br></br></h2>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane</td>
+<td>Sample ID</td>
+<td>Sample Target</td>
+<td>Sample Type</td>
+<td>Length</td>
+<td>Filter</td>
+<td>Tiles</td>
+</tr>
+<tr>
+<td>1</td>
+<td>unknown</td>
+<td>dm3</td>
+<td>ELAND</td>
+<td>32</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td><a href="#Lane1">Lane 1</a></td>
+</tr>
+<tr>
+<td>2</td>
+<td>unknown</td>
+<td>equcab1</td>
+<td>ELAND</td>
+<td>32</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td><a href="#Lane2">Lane 2</a></td>
+</tr>
+<tr>
+<td>3</td>
+<td>unknown</td>
+<td>equcab1</td>
+<td>ELAND</td>
+<td>32</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td><a href="#Lane3">Lane 3</a></td>
+</tr>
+<tr>
+<td>4</td>
+<td>unknown</td>
+<td>canfam2</td>
+<td>ELAND</td>
+<td>32</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td><a href="#Lane4">Lane 4</a></td>
+</tr>
+<tr>
+<td>5</td>
+<td>unknown</td>
+<td>hg18</td>
+<td>ELAND</td>
+<td>32</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td><a href="#Lane5">Lane 5</a></td>
+</tr>
+<tr>
+<td>6</td>
+<td>unknown</td>
+<td>hg18</td>
+<td>ELAND</td>
+<td>32</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td><a href="#Lane6">Lane 6</a></td>
+</tr>
+<tr>
+<td>7</td>
+<td>unknown</td>
+<td>hg18</td>
+<td>ELAND</td>
+<td>32</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td><a href="#Lane7">Lane 7</a></td>
+</tr>
+<tr>
+<td>8</td>
+<td>unknown</td>
+<td>hg18</td>
+<td>ELAND</td>
+<td>32</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td><a href="#Lane8">Lane 8</a></td>
+</tr>
+</table>
+<h2><br></br>Lane Results Summary<br></br></h2>
+<table border="1" cellpadding="5">
+<tr>
+
+<td>Lane </td>
+<td>Clusters </td>
+<td>Av 1st Cycle Int </td>
+<td>% intensity after 20 cycles </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Av Alignment Score (PF) </td>
+<td> % Error Rate (PF) </td>
+</tr>
+<tr>
+<td>1</td>
+<td>17421 +/- 2139</td>
+<td>7230 +/- 801</td>
+<td>23.73 +/- 10.79</td>
+<td>13.00 +/- 22.91</td>
+<td>32.03 +/- 18.45</td>
+<td>6703.57 +/- 3753.85</td>
+<td>4.55 +/- 4.81</td>
+</tr>
+<tr>
+<td>2</td>
+<td>20311 +/- 2402</td>
+<td>7660 +/- 678</td>
+<td>17.03 +/- 4.40</td>
+<td>40.74 +/- 30.33</td>
+<td>29.54 +/- 9.03</td>
+<td>5184.02 +/- 1631.54</td>
+<td>3.27 +/- 3.94</td>
+</tr>
+<tr>
+<td>3</td>
+<td>20193 +/- 2399</td>
+<td>7700 +/- 797</td>
+<td>15.75 +/- 3.30</td>
+<td>56.56 +/- 17.16</td>
+<td>27.33 +/- 7.48</td>
+<td>4803.49 +/- 1313.31</td>
+<td>3.07 +/- 2.86</td>
+</tr>
+<tr>
+<td>4</td>
+<td>15537 +/- 2531</td>
+<td>7620 +/- 1392</td>
+<td>15.37 +/- 3.79</td>
+<td>63.05 +/- 18.30</td>
+<td>15.88 +/- 4.99</td>
+<td>3162.13 +/- 962.59</td>
+<td>3.11 +/- 2.22</td>
+</tr>
+<tr>
+<td>5</td>
+<td>32047 +/- 3356</td>
+<td>8093 +/- 831</td>
+<td>23.79 +/- 6.18</td>
+<td>53.36 +/- 18.06</td>
+<td>48.04 +/- 13.77</td>
+<td>9866.23 +/- 2877.30</td>
+<td>2.26 +/- 1.16</td>
+</tr>
+<tr>
+<td>6</td>
+<td>32946 +/- 4753</td>
+<td>8227 +/- 736</td>
+<td>24.07 +/- 4.69</td>
+<td>54.65 +/- 12.57</td>
+<td>50.98 +/- 10.54</td>
+<td>10468.86 +/- 2228.53</td>
+<td>2.21 +/- 2.33</td>
+</tr>
+<tr>
+<td>7</td>
+<td>39504 +/- 4171</td>
+<td>8401 +/- 785</td>
+<td>22.55 +/- 4.56</td>
+<td>45.22 +/- 10.34</td>
+<td>48.41 +/- 9.67</td>
+<td>9829.40 +/- 1993.20</td>
+<td>2.26 +/- 1.11</td>
+</tr>
+<tr>
+<td>8</td>
+<td>37998 +/- 3792</td>
+<td>8443 +/- 1211</td>
+<td>39.03 +/- 7.52</td>
+<td>42.16 +/- 12.35</td>
+<td>40.98 +/- 14.89</td>
+<td>8128.87 +/- 3055.34</td>
+<td>3.57 +/- 2.77</td>
+</tr>
+</table>
+</body>
+</html>
+"""
+ pathname = os.path.join(gerald_dir, 'Summary.htm')
+ f = open(pathname, 'w')
+ f.write(summary_htm)
+ f.close()
+
+def make_eland_results(gerald_dir):
+ eland_result = """>HWI-EAS229_24_207BTAAXX:1:7:599:759 ACATAGNCACAGACATAAACATAGACATAGAC U0 1 1 3 chrUextra.fa 28189829 R D.
+>HWI-EAS229_24_207BTAAXX:1:7:205:842 AAACAANNCTCCCAAACACGTAAACTGGAAAA U1 0 1 0 chr2L.fa 8796855 R DD 24T
+>HWI-EAS229_24_207BTAAXX:1:7:776:582 AGCTCANCCGATCGAAAACCTCNCCAAGCAAT NM 0 0 0
+>HWI-EAS229_24_207BTAAXX:1:7:205:842 AAACAANNCTCCCAAACACGTAAACTGGAAAA U1 0 1 0 Lambda.fa 8796855 R DD 24T
+"""
+ for i in range(1,9):
+ pathname = os.path.join(gerald_dir,
+ 's_%d_eland_result.txt' % (i,))
+ f = open(pathname, 'w')
+ f.write(eland_result)
+ f.close()
+
+class RunfolderTests(unittest.TestCase):
+ """
+ Test components of the runfolder processing code
+ which includes firecrest, bustard, and gerald
+ """
+ def setUp(self):
+ # make a fake runfolder directory
+ self.temp_dir = tempfile.mkdtemp(prefix='tmp_runfolder_')
+
+ self.runfolder_dir = os.path.join(self.temp_dir,
+ '080102_HWI-EAS229_0010_207BTAAXX')
+ os.mkdir(self.runfolder_dir)
+
+ self.data_dir = os.path.join(self.runfolder_dir, 'Data')
+ os.mkdir(self.data_dir)
+
+ self.firecrest_dir = os.path.join(self.data_dir,
+ 'C1-33_Firecrest1.8.28_12-04-2008_diane'
+ )
+ os.mkdir(self.firecrest_dir)
+ self.matrix_dir = os.path.join(self.firecrest_dir, 'Matrix')
+ os.mkdir(self.matrix_dir)
+ make_matrix(self.matrix_dir)
+
+ self.bustard_dir = os.path.join(self.firecrest_dir,
+ 'Bustard1.8.28_12-04-2008_diane')
+ os.mkdir(self.bustard_dir)
+ make_phasing_params(self.bustard_dir)
+
+ self.gerald_dir = os.path.join(self.bustard_dir,
+ 'GERALD_12-04-2008_diane')
+ os.mkdir(self.gerald_dir)
+ make_gerald_config(self.gerald_dir)
+ make_summary_htm(self.gerald_dir)
+ make_eland_results(self.gerald_dir)
+
+ def tearDown(self):
+ shutil.rmtree(self.temp_dir)
+
+ def test_firecrest(self):
+ """
+ Construct a firecrest object
+ """
+ f = firecrest.firecrest(self.firecrest_dir)
+ self.failUnlessEqual(f.version, '1.8.28')
+ self.failUnlessEqual(f.start, 1)
+ self.failUnlessEqual(f.stop, 33)
+ self.failUnlessEqual(f.user, 'diane')
+ self.failUnlessEqual(f.date, date(2008,4,12))
+
+ xml = f.get_elements()
+ # just make sure that element tree can serialize the tree
+ xml_str = ElementTree.tostring(xml)
+
+ f2 = firecrest.Firecrest(xml=xml)
+ self.failUnlessEqual(f.version, f2.version)
+ self.failUnlessEqual(f.start, f2.start)
+ self.failUnlessEqual(f.stop, f2.stop)
+ self.failUnlessEqual(f.user, f2.user)
+ self.failUnlessEqual(f.date, f2.date)
+
+ def test_bustard(self):
+ """
+ construct a bustard object
+ """
+ b = bustard.bustard(self.bustard_dir)
+ self.failUnlessEqual(b.version, '1.8.28')
+ self.failUnlessEqual(b.date, date(2008,4,12))
+ self.failUnlessEqual(b.user, 'diane')
+ self.failUnlessEqual(len(b.phasing), 8)
+ self.failUnlessAlmostEqual(b.phasing[8].phasing, 0.0099)
+
+ xml = b.get_elements()
+ b2 = bustard.Bustard(xml=xml)
+ self.failUnlessEqual(b.version, b2.version)
+ self.failUnlessEqual(b.date, b2.date )
+ self.failUnlessEqual(b.user, b2.user)
+ self.failUnlessEqual(len(b.phasing), len(b2.phasing))
+ for key in b.phasing.keys():
+ self.failUnlessEqual(b.phasing[key].lane,
+ b2.phasing[key].lane)
+ self.failUnlessEqual(b.phasing[key].phasing,
+ b2.phasing[key].phasing)
+ self.failUnlessEqual(b.phasing[key].prephasing,
+ b2.phasing[key].prephasing)
+
+ def test_gerald(self):
+ # need to update gerald and make tests for it
+ g = gerald.gerald(self.gerald_dir)
+
+ self.failUnlessEqual(g.version,
+ '@(#) Id: GERALD.pl,v 1.68.2.2 2007/06/13 11:08:49 km Exp')
+ self.failUnlessEqual(g.date, datetime(2008,4,19,19,8,30))
+ self.failUnlessEqual(len(g.lanes), len(g.lanes.keys()))
+ self.failUnlessEqual(len(g.lanes), len(g.lanes.items()))
+
+
+ # list of genomes, matches what was defined up in
+ # make_gerald_config.
+ # the first None is to offset the genomes list to be 1..9
+ # instead of pythons default 0..8
+ genomes = [None, '/g/dm3', '/g/equcab1', '/g/equcab1', '/g/canfam2',
+ '/g/hg18', '/g/hg18', '/g/hg18', '/g/hg18', ]
+
+ # test lane specific parameters from gerald config file
+ for i in range(1,9):
+ cur_lane = g.lanes[str(i)]
+ self.failUnlessEqual(cur_lane.analysis, 'eland')
+ self.failUnlessEqual(cur_lane.eland_genome, genomes[i])
+ self.failUnlessEqual(cur_lane.read_length, '32')
+ self.failUnlessEqual(cur_lane.use_bases, 'Y'*32)
+
+ # test data extracted from summary file
+ clusters = [None,
+ (17421, 2139), (20311, 2402), (20193, 2399), (15537, 2531),
+ (32047, 3356), (32946, 4753), (39504, 4171), (37998, 3792)]
+
+ for i in range(1,9):
+ summary_lane = g.summary[str(i)]
+ self.failUnlessEqual(summary_lane.cluster, clusters[i])
+ self.failUnlessEqual(summary_lane.lane, str(i))
+
+ xml = g.get_elements()
+ # just make sure that element tree can serialize the tree
+ xml_str = ElementTree.tostring(xml)
+ g2 = gerald.Gerald(xml=xml)
+
+ # do it all again after extracting from the xml file
+ self.failUnlessEqual(g.version, g2.version)
+ self.failUnlessEqual(g.date, g2.date)
+ self.failUnlessEqual(len(g.lanes.keys()), len(g2.lanes.keys()))
+ self.failUnlessEqual(len(g.lanes.items()), len(g2.lanes.items()))
+
+ # test lane specific parameters from gerald config file
+ for i in range(1,9):
+ g_lane = g.lanes[str(i)]
+ g2_lane = g2.lanes[str(i)]
+ self.failUnlessEqual(g_lane.analysis, g2_lane.analysis)
+ self.failUnlessEqual(g_lane.eland_genome, g2_lane.eland_genome)
+ self.failUnlessEqual(g_lane.read_length, g2_lane.read_length)
+ self.failUnlessEqual(g_lane.use_bases, g2_lane.use_bases)
+
+ # test (some) summary elements
+ for i in range(1,9):
+ g_summary = g.summary[str(i)]
+ g2_summary = g2.summary[str(i)]
+ self.failUnlessEqual(g_summary.cluster, g2_summary.cluster)
+ self.failUnlessEqual(g_summary.lane, g2_summary.lane)
+
+ g_eland = g.eland_results
+ g2_eland = g2.eland_results
+ for lane in g_eland.keys():
+ self.failUnlessEqual(g_eland[lane].reads,
+ g2_eland[lane].reads)
+ self.failUnlessEqual(len(g_eland[lane].mapped_reads),
+ len(g2_eland[lane].mapped_reads))
+ for k in g_eland[lane].mapped_reads.keys():
+ self.failUnlessEqual(g_eland[lane].mapped_reads[k],
+ g2_eland[lane].mapped_reads[k])
+
+ self.failUnlessEqual(len(g_eland[lane].match_codes),
+ len(g2_eland[lane].match_codes))
+ for k in g_eland[lane].match_codes.keys():
+ self.failUnlessEqual(g_eland[lane].match_codes[k],
+ g2_eland[lane].match_codes[k])
+
+
+ def test_eland(self):
+ dm3_map = { 'chrUextra.fa' : 'dm3/chrUextra.fa',
+ 'chr2L.fa': 'dm3/chr2L.fa',
+ 'Lambda.fa': 'Lambda.fa'}
+ genome_maps = { '1':dm3_map, '2':dm3_map, '3':dm3_map, '4':dm3_map,
+ '5':dm3_map, '6':dm3_map, '7':dm3_map, '8':dm3_map }
+ eland = gerald.eland(self.gerald_dir, genome_maps=genome_maps)
+
+ for i in range(1,9):
+ lane = eland[str(i)]
+ self.failUnlessEqual(lane.reads, 4)
+ self.failUnlessEqual(lane.sample_name, "s")
+ self.failUnlessEqual(lane.lane_id, unicode(i))
+ self.failUnlessEqual(len(lane.mapped_reads), 3)
+ self.failUnlessEqual(lane.mapped_reads['Lambda.fa'], 1)
+ self.failUnlessEqual(lane.mapped_reads['dm3/chr2L.fa'], 1)
+ self.failUnlessEqual(lane.match_codes['U1'], 2)
+ self.failUnlessEqual(lane.match_codes['NM'], 1)
+
+ xml = eland.get_elements()
+ # just make sure that element tree can serialize the tree
+ xml_str = ElementTree.tostring(xml)
+ e2 = gerald.ELAND(xml=xml)
+
+ for i in range(1,9):
+ l1 = eland[str(i)]
+ l2 = e2[str(i)]
+ self.failUnlessEqual(l1.reads, l2.reads)
+ self.failUnlessEqual(l1.sample_name, l2.sample_name)
+ self.failUnlessEqual(l1.lane_id, l2.lane_id)
+ self.failUnlessEqual(len(l1.mapped_reads), len(l2.mapped_reads))
+ self.failUnlessEqual(len(l1.mapped_reads), 3)
+ for k in l1.mapped_reads.keys():
+ self.failUnlessEqual(l1.mapped_reads[k],
+ l2.mapped_reads[k])
+
+ self.failUnlessEqual(len(l1.match_codes), 9)
+ self.failUnlessEqual(len(l1.match_codes), len(l2.match_codes))
+ for k in l1.match_codes.keys():
+ self.failUnlessEqual(l1.match_codes[k],
+ l2.match_codes[k])
+
+ def test_runfolder(self):
+ runs = runfolder.get_runs(self.runfolder_dir)
+
+ # do we get the flowcell id from the filename?
+ self.failUnlessEqual(len(runs), 1)
+ self.failUnlessEqual(runs[0].name, 'run_207BTAAXX_2008-04-19.xml')
+
+ # do we get the flowcell id from the FlowcellId.xml file
+ make_flowcell_id(self.runfolder_dir, '207BTAAXY')
+ runs = runfolder.get_runs(self.runfolder_dir)
+ self.failUnlessEqual(len(runs), 1)
+ self.failUnlessEqual(runs[0].name, 'run_207BTAAXY_2008-04-19.xml')
+
+ r1 = runs[0]
+ xml = r1.get_elements()
+ xml_str = ElementTree.tostring(xml)
+
+ r2 = runfolder.PipelineRun(xml=xml)
+ self.failUnlessEqual(r1.name, r2.name)
+ self.failIfEqual(r2.firecrest, None)
+ self.failIfEqual(r2.bustard, None)
+ self.failIfEqual(r2.gerald, None)
+
+
+def suite():
+ return unittest.makeSuite(RunfolderTests,'test')
+
+if __name__ == "__main__":
+ unittest.main(defaultTest="suite")
+
--- /dev/null
+#!/usr/bin/env python
+
+from datetime import datetime, date
+import os
+import tempfile
+import shutil
+import unittest
+
+from htswdataprod.illumina import firecrest
+from htswdataprod.illumina import bustard
+from htswdataprod.illumina import gerald
+from htswdataprod.illumina import runfolder
+from htswdataprod.illumina.runfolder import ElementTree
+
+
+def make_flowcell_id(runfolder_dir, flowcell_id=None):
+ if flowcell_id is None:
+ flowcell_id = '207BTAAXY'
+
+ config = """<?xml version="1.0"?>
+<FlowcellId>
+ <Text>%s</Text>
+</FlowcellId>""" % (flowcell_id,)
+ config_dir = os.path.join(runfolder_dir, 'Config')
+
+ if not os.path.exists(config_dir):
+ os.mkdir(config_dir)
+ pathname = os.path.join(config_dir, 'FlowcellId.xml')
+ f = open(pathname,'w')
+ f.write(config)
+ f.close()
+
+def make_matrix(matrix_dir):
+ contents = """# Auto-generated frequency response matrix
+> A
+> C
+> G
+> T
+0.77 0.15 -0.04 -0.04
+0.76 1.02 -0.05 -0.06
+-0.10 -0.10 1.17 -0.03
+-0.13 -0.12 0.80 1.27
+"""
+ s_matrix = os.path.join(matrix_dir, 's_matrix.txt')
+ f = open(s_matrix, 'w')
+ f.write(contents)
+ f.close()
+
+def make_phasing_params(bustard_dir):
+ for lane in range(1,9):
+ pathname = os.path.join(bustard_dir, 'params%d.xml' % (lane))
+ f = open(pathname, 'w')
+ f.write("""<Parameters>
+ <Phasing>0.009900</Phasing>
+ <Prephasing>0.003500</Prephasing>
+</Parameters>
+""")
+ f.close()
+
+def make_gerald_config(gerald_dir):
+ config_xml = """<RunParameters>
+<ChipWideRunParameters>
+ <ANALYSIS>default</ANALYSIS>
+ <BAD_LANES></BAD_LANES>
+ <BAD_TILES></BAD_TILES>
+ <CONTAM_DIR></CONTAM_DIR>
+ <CONTAM_FILE></CONTAM_FILE>
+ <ELAND_GENOME>Need_to_specify_ELAND_genome_directory</ELAND_GENOME>
+ <ELAND_MULTIPLE_INSTANCES>8</ELAND_MULTIPLE_INSTANCES>
+ <ELAND_REPEAT></ELAND_REPEAT>
+ <EMAIL_DOMAIN>domain.com</EMAIL_DOMAIN>
+ <EMAIL_LIST>diane</EMAIL_LIST>
+ <EMAIL_SERVER>localhost:25</EMAIL_SERVER>
+ <EXPT_DIR>/home/diane/gec/080416_HWI-EAS229_0024_207BTAAXX/Data/C1-33_Firecrest1.8.28_19-04-2008_diane/Bustard1.8.28_19-04-2008_diane</EXPT_DIR>
+ <EXPT_DIR_ROOT>/home/diane/gec</EXPT_DIR_ROOT>
+ <FORCE>1</FORCE>
+ <GENOME_DIR>/home/diane/proj/SolexaPipeline-0.2.2.6/Goat/../Gerald/../../Genomes</GENOME_DIR>
+ <GENOME_FILE>Need_to_specify_genome_file_name</GENOME_FILE>
+ <HAMSTER_FLAG>genome</HAMSTER_FLAG>
+ <OUT_DIR>/home/diane/gec/080416_HWI-EAS229_0024_207BTAAXX/Data/C1-33_Firecrest1.8.28_19-04-2008_diane/Bustard1.8.28_19-04-2008_diane/GERALD_19-04-2008_diane</OUT_DIR>
+ <POST_RUN_COMMAND></POST_RUN_COMMAND>
+ <PRB_FILE_SUFFIX>_prb.txt</PRB_FILE_SUFFIX>
+ <PURE_BASES>12</PURE_BASES>
+ <QF_PARAMS>'((CHASTITY>=0.6))'</QF_PARAMS>
+ <QHG_FILE_SUFFIX>_qhg.txt</QHG_FILE_SUFFIX>
+ <QUALITY_FORMAT>--symbolic</QUALITY_FORMAT>
+ <READ_LENGTH>32</READ_LENGTH>
+ <SEQUENCE_FORMAT>--scarf</SEQUENCE_FORMAT>
+ <SEQ_FILE_SUFFIX>_seq.txt</SEQ_FILE_SUFFIX>
+ <SIG_FILE_SUFFIX_DEPHASED>_sig2.txt</SIG_FILE_SUFFIX_DEPHASED>
+ <SIG_FILE_SUFFIX_NOT_DEPHASED>_sig.txt</SIG_FILE_SUFFIX_NOT_DEPHASED>
+ <SOFTWARE_VERSION>@(#) Id: GERALD.pl,v 1.68.2.2 2007/06/13 11:08:49 km Exp</SOFTWARE_VERSION>
+ <TILE_REGEX>s_[1-8]_[0-9][0-9][0-9][0-9]</TILE_REGEX>
+ <TILE_ROOT>s</TILE_ROOT>
+ <TIME_STAMP>Sat Apr 19 19:08:30 2008</TIME_STAMP>
+ <TOOLS_DIR>/home/diane/proj/SolexaPipeline-0.2.2.6/Goat/../Gerald</TOOLS_DIR>
+ <USE_BASES>all</USE_BASES>
+ <WEB_DIR_ROOT>http://host.domain.com/yourshare/</WEB_DIR_ROOT>
+</ChipWideRunParameters>
+<LaneSpecificRunParameters>
+ <ANALYSIS>
+ <s_1>eland</s_1>
+ <s_2>eland</s_2>
+ <s_3>eland</s_3>
+ <s_4>eland</s_4>
+ <s_5>eland</s_5>
+ <s_6>eland</s_6>
+ <s_7>eland</s_7>
+ <s_8>eland</s_8>
+ </ANALYSIS>
+ <ELAND_GENOME>
+ <s_1>/g/dm3</s_1>
+ <s_2>/g/equcab1</s_2>
+ <s_3>/g/equcab1</s_3>
+ <s_4>/g/canfam2</s_4>
+ <s_5>/g/hg18</s_5>
+ <s_6>/g/hg18</s_6>
+ <s_7>/g/hg18</s_7>
+ <s_8>/g/hg18</s_8>
+ </ELAND_GENOME>
+ <READ_LENGTH>
+ <s_1>32</s_1>
+ <s_2>32</s_2>
+ <s_3>32</s_3>
+ <s_4>32</s_4>
+ <s_5>32</s_5>
+ <s_6>32</s_6>
+ <s_7>32</s_7>
+ <s_8>32</s_8>
+ </READ_LENGTH>
+ <USE_BASES>
+ <s_1>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_1>
+ <s_2>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_2>
+ <s_3>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_3>
+ <s_4>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_4>
+ <s_5>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_5>
+ <s_6>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_6>
+ <s_7>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_7>
+ <s_8>YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY</s_8>
+ </USE_BASES>
+</LaneSpecificRunParameters>
+</RunParameters>
+"""
+ pathname = os.path.join(gerald_dir, 'config.xml')
+ f = open(pathname,'w')
+ f.write(config_xml)
+ f.close()
+
+def make_summary_htm(gerald_dir):
+ summary_htm="""<!--RUN_TIME Wed Jul 2 06:47:44 2008 -->
+<!--SOFTWARE_VERSION @(#) $Id: jerboa.pl,v 1.94 2007/12/04 09:59:07 rshaw Exp $-->
+<html>
+<body>
+
+<a name="Top"><h2><title>080627_HWI-EAS229_0036_3055HAXX Summary</title></h2></a>
+<h1>Summary Information For Experiment 080627_HWI-EAS229_0036_3055HAXX on Machine HWI-EAS229</h1>
+<h2><br></br>Chip Summary<br></br></h2>
+<table border="1" cellpadding="5">
+<tr><td>Machine</td><td>HWI-EAS229</td></tr>
+<tr><td>Run Folder</td><td>080627_HWI-EAS229_0036_3055HAXX</td></tr>
+<tr><td>Chip ID</td><td>unknown</td></tr>
+</table>
+<h2><br></br>Chip Results Summary<br></br></h2>
+<table border="1" cellpadding="5">
+<tr>
+<td>Clusters</td>
+<td>Clusters (PF)</td>
+<td>Yield (kbases)</td>
+</tr>
+<tr><td>80933224</td>
+<td>43577803</td>
+<td>1133022</td>
+</tr>
+</table>
+<h2><br></br>Lane Parameter Summary<br></br></h2>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane</td>
+<td>Sample ID</td>
+<td>Sample Target</td>
+<td>Sample Type</td>
+<td>Length</td>
+<td>Filter</td>
+<td>Num Tiles</td>
+<td>Tiles</td>
+</tr>
+<tr>
+<td>1</td>
+<td>unknown</td>
+<td>mm9</td>
+<td>ELAND</td>
+<td>26</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td>100</td>
+<td><a href="#Lane1">Lane 1</a></td>
+</tr>
+<tr>
+<td>2</td>
+<td>unknown</td>
+<td>mm9</td>
+<td>ELAND</td>
+<td>26</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td>100</td>
+<td><a href="#Lane2">Lane 2</a></td>
+</tr>
+<tr>
+<td>3</td>
+<td>unknown</td>
+<td>mm9</td>
+<td>ELAND</td>
+<td>26</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td>100</td>
+<td><a href="#Lane3">Lane 3</a></td>
+</tr>
+<tr>
+<td>4</td>
+<td>unknown</td>
+<td>elegans170</td>
+<td>ELAND</td>
+<td>26</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td>100</td>
+<td><a href="#Lane4">Lane 4</a></td>
+</tr>
+<tr>
+<td>5</td>
+<td>unknown</td>
+<td>elegans170</td>
+<td>ELAND</td>
+<td>26</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td>100</td>
+<td><a href="#Lane5">Lane 5</a></td>
+</tr>
+<tr>
+<td>6</td>
+<td>unknown</td>
+<td>elegans170</td>
+<td>ELAND</td>
+<td>26</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td>100</td>
+<td><a href="#Lane6">Lane 6</a></td>
+</tr>
+<tr>
+<td>7</td>
+<td>unknown</td>
+<td>elegans170</td>
+<td>ELAND</td>
+<td>26</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td>100</td>
+<td><a href="#Lane7">Lane 7</a></td>
+</tr>
+<tr>
+<td>8</td>
+<td>unknown</td>
+<td>elegans170</td>
+<td>ELAND</td>
+<td>26</td>
+<td>'((CHASTITY>=0.6))'</td>
+<td>100</td>
+<td><a href="#Lane8">Lane 8</a></td>
+</tr>
+</table>
+<h2><br></br>Lane Results Summary<br></br></h2>
+<table border="1" cellpadding="5">
+<tr>
+<td colspan="2">Lane Info</td>
+<td colspan="8">Tile Mean +/- SD for Lane</td>
+</tr>
+<tr>
+<td>Lane </td>
+<td>Lane Yield (kbases) </td>
+<td>Clusters (raw)</td>
+<td>Clusters (PF) </td>
+<td>1st Cycle Int (PF) </td>
+<td>% intensity after 20 cycles (PF) </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Alignment Score (PF) </td>
+<td> % Error Rate (PF) </td>
+</tr>
+<tr>
+<td>1</td>
+<td>158046</td>
+<td>96483 +/- 9074</td>
+<td>60787 +/- 4240</td>
+<td>329 +/- 35</td>
+<td>101.88 +/- 6.03</td>
+<td>63.21 +/- 3.29</td>
+<td>70.33 +/- 0.24</td>
+<td>9054.08 +/- 59.16</td>
+<td>0.46 +/- 0.18</td>
+</tr>
+<tr>
+<td>2</td>
+<td>156564</td>
+<td>133738 +/- 7938</td>
+<td>60217 +/- 1926</td>
+<td>444 +/- 39</td>
+<td>92.62 +/- 7.58</td>
+<td>45.20 +/- 3.31</td>
+<td>51.98 +/- 0.74</td>
+<td>6692.04 +/- 92.49</td>
+<td>0.46 +/- 0.09</td>
+</tr>
+<tr>
+<td>3</td>
+<td>185818</td>
+<td>152142 +/- 10002</td>
+<td>71468 +/- 2827</td>
+<td>366 +/- 36</td>
+<td>91.53 +/- 8.66</td>
+<td>47.19 +/- 3.80</td>
+<td>82.24 +/- 0.44</td>
+<td>10598.68 +/- 64.13</td>
+<td>0.41 +/- 0.04</td>
+</tr>
+<tr>
+<td>4</td>
+<td>34953</td>
+<td>15784 +/- 2162</td>
+<td>13443 +/- 1728</td>
+<td>328 +/- 40</td>
+<td>97.53 +/- 9.87</td>
+<td>85.29 +/- 1.91</td>
+<td>80.02 +/- 0.53</td>
+<td>10368.82 +/- 71.08</td>
+<td>0.15 +/- 0.05</td>
+</tr>
+<tr>
+<td>5</td>
+<td>167936</td>
+<td>119735 +/- 8465</td>
+<td>64590 +/- 2529</td>
+<td>417 +/- 37</td>
+<td>88.69 +/- 14.79</td>
+<td>54.10 +/- 2.59</td>
+<td>76.95 +/- 0.32</td>
+<td>9936.47 +/- 65.75</td>
+<td>0.28 +/- 0.02</td>
+</tr>
+<tr>
+<td>6</td>
+<td>173463</td>
+<td>152177 +/- 8146</td>
+<td>66716 +/- 2493</td>
+<td>372 +/- 39</td>
+<td>87.06 +/- 9.86</td>
+<td>43.98 +/- 3.12</td>
+<td>78.80 +/- 0.43</td>
+<td>10162.28 +/- 49.65</td>
+<td>0.38 +/- 0.03</td>
+</tr>
+<tr>
+<td>7</td>
+<td>149287</td>
+<td>84649 +/- 7325</td>
+<td>57418 +/- 3617</td>
+<td>295 +/- 28</td>
+<td>89.40 +/- 8.23</td>
+<td>67.97 +/- 1.82</td>
+<td>33.38 +/- 0.25</td>
+<td>4247.92 +/- 32.37</td>
+<td>1.00 +/- 0.03</td>
+</tr>
+<tr>
+<td>8</td>
+<td>106953</td>
+<td>54622 +/- 4812</td>
+<td>41136 +/- 3309</td>
+<td>284 +/- 37</td>
+<td>90.21 +/- 9.10</td>
+<td>75.39 +/- 2.27</td>
+<td>48.33 +/- 0.29</td>
+<td>6169.21 +/- 169.50</td>
+<td>0.86 +/- 1.22</td>
+</tr>
+<tr><td colspan="13">Tile mean across chip</td></tr>
+<tr>
+<td>Av.</td>
+<td></td>
+<td>101166</td>
+<td>54472</td>
+<td>354</td>
+<td>92.36</td>
+<td>60.29</td>
+<td>65.25</td>
+<td>8403.69</td>
+<td>0.50</td>
+</tr>
+</table>
+<h2><br></br>Expanded Lane Summary<br></br></h2>
+<table border="1" cellpadding="5">
+<tr>
+
+<tr><td colspan="2">Lane Info</td>
+<td colspan="2">Phasing Info</td>
+<td colspan="2">Raw Data (tile mean)</td>
+<td colspan="7">Filtered Data (tile mean)</td></tr>
+<td>Lane </td>
+<td>Clusters (tile mean) (raw)</td>
+<td>% Phasing </td>
+<td>% Prephasing </td>
+<td>% Error Rate (raw) </td>
+<td> Equiv Perfect Clusters (raw) </td>
+<td>% retained </td>
+<td>Cycle 2-4 Av Int (PF) </td>
+<td>Cycle 2-10 Av % Loss (PF) </td>
+<td>Cycle 10-20 Av % Loss (PF) </td>
+<td>% Align (PF) </td>
+<td>% Error Rate (PF) </td>
+<td> Equiv Perfect Clusters (PF) </td>
+</tr>
+<tr>
+<td>1</td>
+<td>96483</td>
+<td>0.7700</td>
+<td>0.3100</td>
+<td>1.00</td>
+<td>49676</td>
+<td>63.21</td>
+<td>317 +/- 32</td>
+<td>0.13 +/- 0.44</td>
+<td>-1.14 +/- 0.34</td>
+<td>70.33</td>
+<td>0.46</td>
+<td>41758</td>
+</tr>
+<tr>
+<td>2</td>
+<td>133738</td>
+<td>0.7700</td>
+<td>0.3100</td>
+<td>1.22</td>
+<td>40467</td>
+<td>45.20</td>
+<td>415 +/- 33</td>
+<td>0.29 +/- 0.40</td>
+<td>-0.79 +/- 0.35</td>
+<td>51.98</td>
+<td>0.46</td>
+<td>30615</td>
+</tr>
+<tr>
+<td>3</td>
+<td>152142</td>
+<td>0.7700</td>
+<td>0.3100</td>
+<td>1.30</td>
+<td>78588</td>
+<td>47.19</td>
+<td>344 +/- 26</td>
+<td>0.68 +/- 0.51</td>
+<td>-0.77 +/- 0.42</td>
+<td>82.24</td>
+<td>0.41</td>
+<td>57552</td>
+</tr>
+<tr>
+<td>4</td>
+<td>15784</td>
+<td>0.7700</td>
+<td>0.3100</td>
+<td>0.29</td>
+<td>11095</td>
+<td>85.29</td>
+<td>306 +/- 34</td>
+<td>0.20 +/- 0.69</td>
+<td>-1.28 +/- 0.66</td>
+<td>80.02</td>
+<td>0.15</td>
+<td>10671</td>
+</tr>
+<tr>
+<td>5</td>
+<td>119735</td>
+<td>0.7700</td>
+<td>0.3100</td>
+<td>0.85</td>
+<td>60335</td>
+<td>54.10</td>
+<td>380 +/- 32</td>
+<td>0.34 +/- 0.49</td>
+<td>-1.55 +/- 4.69</td>
+<td>76.95</td>
+<td>0.28</td>
+<td>49015</td>
+</tr>
+<tr>
+<td>6</td>
+<td>152177</td>
+<td>0.7700</td>
+<td>0.3100</td>
+<td>1.21</td>
+<td>70905</td>
+<td>43.98</td>
+<td>333 +/- 27</td>
+<td>0.57 +/- 0.50</td>
+<td>-0.91 +/- 0.39</td>
+<td>78.80</td>
+<td>0.38</td>
+<td>51663</td>
+</tr>
+<tr>
+<td>7</td>
+<td>84649</td>
+<td>0.7700</td>
+<td>0.3100</td>
+<td>1.38</td>
+<td>21069</td>
+<td>67.97</td>
+<td>272 +/- 20</td>
+<td>1.15 +/- 0.52</td>
+<td>-0.84 +/- 0.58</td>
+<td>33.38</td>
+<td>1.00</td>
+<td>18265</td>
+</tr>
+<tr>
+<td>8</td>
+<td>54622</td>
+<td>0.7700</td>
+<td>0.3100</td>
+<td>1.17</td>
+<td>21335</td>
+<td>75.39</td>
+<td>262 +/- 31</td>
+<td>1.10 +/- 0.59</td>
+<td>-1.01 +/- 0.47</td>
+<td>48.33</td>
+<td>0.86</td>
+<td>19104</td>
+</tr>
+</table>
+<b><br></br>IVC Plots</b>
+<p> <a href='IVC.htm' target="_blank"> IVC.htm
+ </a></p>
+<b><br></br>All Intensity Plots</b>
+<p> <a href='All.htm' target="_blank"> All.htm
+ </a></p>
+<b><br></br>Error graphs: </b>
+<p> <a href='Error.htm' target="_blank"> Error.htm
+ </a></p>
+<td><a href="#Top">Back to top</a></td>
+<a name="Lane1"><h2><br></br>Lane 1<br></br></h2></a>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane </td>
+<td>Tile </td>
+<td>Clusters (raw)</td>
+<td>Av 1st Cycle Int (PF) </td>
+<td>Av % intensity after 20 cycles (PF) </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Av Alignment Score (PF) </td>
+<td>% Error Rate (PF) </td>
+</tr>
+<tr>
+<td>1</td>
+<td>0001</td>
+<td>114972</td>
+<td>326.48</td>
+<td>94.39</td>
+<td>57.44</td>
+<td>70.2</td>
+<td>9038.6</td>
+<td>0.44</td>
+</tr>
+</table>
+<td><a href="#Top">Back to top</a></td>
+<a name="Lane2"><h2><br></br>Lane 2<br></br></h2></a>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane </td>
+<td>Tile </td>
+<td>Clusters (raw)</td>
+<td>Av 1st Cycle Int (PF) </td>
+<td>Av % intensity after 20 cycles (PF) </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Av Alignment Score (PF) </td>
+<td>% Error Rate (PF) </td>
+</tr>
+<tr>
+<td>2</td>
+<td>0001</td>
+<td>147793</td>
+<td>448.12</td>
+<td>83.68</td>
+<td>38.57</td>
+<td>53.7</td>
+<td>6905.4</td>
+<td>0.54</td>
+</tr>
+</table>
+<td><a href="#Top">Back to top</a></td>
+<a name="Lane3"><h2><br></br>Lane 3<br></br></h2></a>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane </td>
+<td>Tile </td>
+<td>Clusters (raw)</td>
+<td>Av 1st Cycle Int (PF) </td>
+<td>Av % intensity after 20 cycles (PF) </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Av Alignment Score (PF) </td>
+<td>% Error Rate (PF) </td>
+</tr>
+<tr>
+<td>3</td>
+<td>0001</td>
+<td>167904</td>
+<td>374.05</td>
+<td>86.91</td>
+<td>40.36</td>
+<td>81.3</td>
+<td>10465.0</td>
+<td>0.47</td>
+</tr>
+</table>
+<td><a href="#Top">Back to top</a></td>
+<a name="Lane4"><h2><br></br>Lane 4<br></br></h2></a>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane </td>
+<td>Tile </td>
+<td>Clusters (raw)</td>
+<td>Av 1st Cycle Int (PF) </td>
+<td>Av % intensity after 20 cycles (PF) </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Av Alignment Score (PF) </td>
+<td>% Error Rate (PF) </td>
+</tr>
+<tr>
+<td>4</td>
+<td>0001</td>
+<td>20308</td>
+<td>276.85</td>
+<td>92.87</td>
+<td>84.26</td>
+<td>80.4</td>
+<td>10413.8</td>
+<td>0.16</td>
+</tr>
+</table>
+<td><a href="#Top">Back to top</a></td>
+<a name="Lane5"><h2><br></br>Lane 5<br></br></h2></a>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane </td>
+<td>Tile </td>
+<td>Clusters (raw)</td>
+<td>Av 1st Cycle Int (PF) </td>
+<td>Av % intensity after 20 cycles (PF) </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Av Alignment Score (PF) </td>
+<td>% Error Rate (PF) </td>
+</tr>
+</table>
+<td><a href="#Top">Back to top</a></td>
+<a name="Lane6"><h2><br></br>Lane 6<br></br></h2></a>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane </td>
+<td>Tile </td>
+<td>Clusters (raw)</td>
+<td>Av 1st Cycle Int (PF) </td>
+<td>Av % intensity after 20 cycles (PF) </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Av Alignment Score (PF) </td>
+<td>% Error Rate (PF) </td>
+</tr>
+<tr>
+<td>6</td>
+<td>0001</td>
+<td>166844</td>
+<td>348.12</td>
+<td>77.59</td>
+<td>38.13</td>
+<td>79.7</td>
+<td>10264.4</td>
+<td>0.44</td>
+</tr>
+</table>
+<td><a href="#Top">Back to top</a></td>
+<a name="Lane7"><h2><br></br>Lane 7<br></br></h2></a>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane </td>
+<td>Tile </td>
+<td>Clusters (raw)</td>
+<td>Av 1st Cycle Int (PF) </td>
+<td>Av % intensity after 20 cycles (PF) </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Av Alignment Score (PF) </td>
+<td>% Error Rate (PF) </td>
+</tr>
+<tr>
+<td>7</td>
+<td>0001</td>
+<td>98913</td>
+<td>269.90</td>
+<td>86.66</td>
+<td>64.55</td>
+<td>33.2</td>
+<td>4217.5</td>
+<td>1.02</td>
+</tr>
+</table>
+<td><a href="#Top">Back to top</a></td>
+<a name="Lane8"><h2><br></br>Lane 8<br></br></h2></a>
+<table border="1" cellpadding="5">
+<tr>
+<td>Lane </td>
+<td>Tile </td>
+<td>Clusters (raw)</td>
+<td>Av 1st Cycle Int (PF) </td>
+<td>Av % intensity after 20 cycles (PF) </td>
+<td>% PF Clusters </td>
+<td>% Align (PF) </td>
+<td>Av Alignment Score (PF) </td>
+<td>% Error Rate (PF) </td>
+</tr>
+<tr>
+<td>8</td>
+<td>0001</td>
+<td>64972</td>
+<td>243.60</td>
+<td>89.40</td>
+<td>73.17</td>
+<td>48.3</td>
+<td>6182.8</td>
+<td>0.71</td>
+</tr>
+</table>
+<td><a href="#Top">Back to top</a></td>
+</body>
+</html>
+"""
+ pathname = os.path.join(gerald_dir, 'Summary.htm')
+ f = open(pathname, 'w')
+ f.write(summary_htm)
+ f.close()
+
+def make_eland_results(gerald_dir):
+ eland_result = """>HWI-EAS229_24_207BTAAXX:1:7:599:759 ACATAGNCACAGACATAAACATAGACATAGAC U0 1 1 3 chrUextra.fa 28189829 R D.
+>HWI-EAS229_24_207BTAAXX:1:7:205:842 AAACAANNCTCCCAAACACGTAAACTGGAAAA U1 0 1 0 chr2L.fa 8796855 R DD 24T
+>HWI-EAS229_24_207BTAAXX:1:7:776:582 AGCTCANCCGATCGAAAACCTCNCCAAGCAAT NM 0 0 0
+>HWI-EAS229_24_207BTAAXX:1:7:205:842 AAACAANNCTCCCAAACACGTAAACTGGAAAA U1 0 1 0 Lambda.fa 8796855 R DD 24T
+"""
+ for i in range(1,9):
+ pathname = os.path.join(gerald_dir,
+ 's_%d_eland_result.txt' % (i,))
+ f = open(pathname, 'w')
+ f.write(eland_result)
+ f.close()
+
+def make_runfolder(obj=None):
+ """
+ Make a fake runfolder, attach all the directories to obj if defined
+ """
+ # make a fake runfolder directory
+ temp_dir = tempfile.mkdtemp(prefix='tmp_runfolder_')
+
+ runfolder_dir = os.path.join(temp_dir,
+ '080102_HWI-EAS229_0010_207BTAAXX')
+ os.mkdir(runfolder_dir)
+
+ data_dir = os.path.join(runfolder_dir, 'Data')
+ os.mkdir(data_dir)
+
+ firecrest_dir = os.path.join(data_dir,
+ 'C1-33_Firecrest1.8.28_12-04-2008_diane'
+ )
+ os.mkdir(firecrest_dir)
+ matrix_dir = os.path.join(firecrest_dir, 'Matrix')
+ os.mkdir(matrix_dir)
+ make_matrix(matrix_dir)
+
+ bustard_dir = os.path.join(firecrest_dir,
+ 'Bustard1.8.28_12-04-2008_diane')
+ os.mkdir(bustard_dir)
+ make_phasing_params(bustard_dir)
+
+ gerald_dir = os.path.join(bustard_dir,
+ 'GERALD_12-04-2008_diane')
+ os.mkdir(gerald_dir)
+ make_gerald_config(gerald_dir)
+ make_summary_htm(gerald_dir)
+ make_eland_results(gerald_dir)
+
+ if obj is not None:
+ obj.temp_dir = temp_dir
+ obj.runfolder_dir = runfolder_dir
+ obj.data_dir = data_dir
+ obj.firecrest_dir = firecrest_dir
+ obj.matrix_dir = matrix_dir
+ obj.bustard_dir = bustard_dir
+ obj.gerald_dir = gerald_dir
+
+
+class RunfolderTests(unittest.TestCase):
+ """
+ Test components of the runfolder processing code
+ which includes firecrest, bustard, and gerald
+ """
+ def setUp(self):
+ # attaches all the directories to the object passed in
+ make_runfolder(self)
+
+ def tearDown(self):
+ shutil.rmtree(self.temp_dir)
+
+ def test_firecrest(self):
+ """
+ Construct a firecrest object
+ """
+ f = firecrest.firecrest(self.firecrest_dir)
+ self.failUnlessEqual(f.version, '1.8.28')
+ self.failUnlessEqual(f.start, 1)
+ self.failUnlessEqual(f.stop, 33)
+ self.failUnlessEqual(f.user, 'diane')
+ self.failUnlessEqual(f.date, date(2008,4,12))
+
+ xml = f.get_elements()
+ # just make sure that element tree can serialize the tree
+ xml_str = ElementTree.tostring(xml)
+
+ f2 = firecrest.Firecrest(xml=xml)
+ self.failUnlessEqual(f.version, f2.version)
+ self.failUnlessEqual(f.start, f2.start)
+ self.failUnlessEqual(f.stop, f2.stop)
+ self.failUnlessEqual(f.user, f2.user)
+ self.failUnlessEqual(f.date, f2.date)
+
+ def test_bustard(self):
+ """
+ construct a bustard object
+ """
+ b = bustard.bustard(self.bustard_dir)
+ self.failUnlessEqual(b.version, '1.8.28')
+ self.failUnlessEqual(b.date, date(2008,4,12))
+ self.failUnlessEqual(b.user, 'diane')
+ self.failUnlessEqual(len(b.phasing), 8)
+ self.failUnlessAlmostEqual(b.phasing[8].phasing, 0.0099)
+
+ xml = b.get_elements()
+ b2 = bustard.Bustard(xml=xml)
+ self.failUnlessEqual(b.version, b2.version)
+ self.failUnlessEqual(b.date, b2.date )
+ self.failUnlessEqual(b.user, b2.user)
+ self.failUnlessEqual(len(b.phasing), len(b2.phasing))
+ for key in b.phasing.keys():
+ self.failUnlessEqual(b.phasing[key].lane,
+ b2.phasing[key].lane)
+ self.failUnlessEqual(b.phasing[key].phasing,
+ b2.phasing[key].phasing)
+ self.failUnlessEqual(b.phasing[key].prephasing,
+ b2.phasing[key].prephasing)
+
+ def test_gerald(self):
+ # need to update gerald and make tests for it
+ g = gerald.gerald(self.gerald_dir)
+
+ self.failUnlessEqual(g.version,
+ '@(#) Id: GERALD.pl,v 1.68.2.2 2007/06/13 11:08:49 km Exp')
+ self.failUnlessEqual(g.date, datetime(2008,4,19,19,8,30))
+ self.failUnlessEqual(len(g.lanes), len(g.lanes.keys()))
+ self.failUnlessEqual(len(g.lanes), len(g.lanes.items()))
+
+
+ # list of genomes, matches what was defined up in
+ # make_gerald_config.
+ # the first None is to offset the genomes list to be 1..9
+ # instead of pythons default 0..8
+ genomes = [None, '/g/dm3', '/g/equcab1', '/g/equcab1', '/g/canfam2',
+ '/g/hg18', '/g/hg18', '/g/hg18', '/g/hg18', ]
+
+ # test lane specific parameters from gerald config file
+ for i in range(1,9):
+ cur_lane = g.lanes[str(i)]
+ self.failUnlessEqual(cur_lane.analysis, 'eland')
+ self.failUnlessEqual(cur_lane.eland_genome, genomes[i])
+ self.failUnlessEqual(cur_lane.read_length, '32')
+ self.failUnlessEqual(cur_lane.use_bases, 'Y'*32)
+
+ # test data extracted from summary file
+ clusters = [None,
+ (96483, 9074), (133738, 7938),
+ (152142, 10002), (15784, 2162),
+ (119735, 8465), (152177, 8146),
+ (84649, 7325), (54622, 4812),]
+
+ for i in range(1,9):
+ summary_lane = g.summary[str(i)]
+ self.failUnlessEqual(summary_lane.cluster, clusters[i])
+ self.failUnlessEqual(summary_lane.lane, str(i))
+
+ xml = g.get_elements()
+ # just make sure that element tree can serialize the tree
+ xml_str = ElementTree.tostring(xml)
+ g2 = gerald.Gerald(xml=xml)
+
+ # do it all again after extracting from the xml file
+ self.failUnlessEqual(g.version, g2.version)
+ self.failUnlessEqual(g.date, g2.date)
+ self.failUnlessEqual(len(g.lanes.keys()), len(g2.lanes.keys()))
+ self.failUnlessEqual(len(g.lanes.items()), len(g2.lanes.items()))
+
+ # test lane specific parameters from gerald config file
+ for i in range(1,9):
+ g_lane = g.lanes[str(i)]
+ g2_lane = g2.lanes[str(i)]
+ self.failUnlessEqual(g_lane.analysis, g2_lane.analysis)
+ self.failUnlessEqual(g_lane.eland_genome, g2_lane.eland_genome)
+ self.failUnlessEqual(g_lane.read_length, g2_lane.read_length)
+ self.failUnlessEqual(g_lane.use_bases, g2_lane.use_bases)
+
+ # test (some) summary elements
+ for i in range(1,9):
+ g_summary = g.summary[str(i)]
+ g2_summary = g2.summary[str(i)]
+ self.failUnlessEqual(g_summary.cluster, g2_summary.cluster)
+ self.failUnlessEqual(g_summary.lane, g2_summary.lane)
+
+ g_eland = g.eland_results
+ g2_eland = g2.eland_results
+ for lane in g_eland.keys():
+ self.failUnlessEqual(g_eland[lane].reads,
+ g2_eland[lane].reads)
+ self.failUnlessEqual(len(g_eland[lane].mapped_reads),
+ len(g2_eland[lane].mapped_reads))
+ for k in g_eland[lane].mapped_reads.keys():
+ self.failUnlessEqual(g_eland[lane].mapped_reads[k],
+ g2_eland[lane].mapped_reads[k])
+
+ self.failUnlessEqual(len(g_eland[lane].match_codes),
+ len(g2_eland[lane].match_codes))
+ for k in g_eland[lane].match_codes.keys():
+ self.failUnlessEqual(g_eland[lane].match_codes[k],
+ g2_eland[lane].match_codes[k])
+
+
+ def test_eland(self):
+ dm3_map = { 'chrUextra.fa' : 'dm3/chrUextra.fa',
+ 'chr2L.fa': 'dm3/chr2L.fa',
+ 'Lambda.fa': 'Lambda.fa'}
+ genome_maps = { '1':dm3_map, '2':dm3_map, '3':dm3_map, '4':dm3_map,
+ '5':dm3_map, '6':dm3_map, '7':dm3_map, '8':dm3_map }
+ eland = gerald.eland(self.gerald_dir, genome_maps=genome_maps)
+
+ for i in range(1,9):
+ lane = eland[str(i)]
+ self.failUnlessEqual(lane.reads, 4)
+ self.failUnlessEqual(lane.sample_name, "s")
+ self.failUnlessEqual(lane.lane_id, unicode(i))
+ self.failUnlessEqual(len(lane.mapped_reads), 3)
+ self.failUnlessEqual(lane.mapped_reads['Lambda.fa'], 1)
+ self.failUnlessEqual(lane.mapped_reads['dm3/chr2L.fa'], 1)
+ self.failUnlessEqual(lane.match_codes['U1'], 2)
+ self.failUnlessEqual(lane.match_codes['NM'], 1)
+
+ xml = eland.get_elements()
+ # just make sure that element tree can serialize the tree
+ xml_str = ElementTree.tostring(xml)
+ e2 = gerald.ELAND(xml=xml)
+
+ for i in range(1,9):
+ l1 = eland[str(i)]
+ l2 = e2[str(i)]
+ self.failUnlessEqual(l1.reads, l2.reads)
+ self.failUnlessEqual(l1.sample_name, l2.sample_name)
+ self.failUnlessEqual(l1.lane_id, l2.lane_id)
+ self.failUnlessEqual(len(l1.mapped_reads), len(l2.mapped_reads))
+ self.failUnlessEqual(len(l1.mapped_reads), 3)
+ for k in l1.mapped_reads.keys():
+ self.failUnlessEqual(l1.mapped_reads[k],
+ l2.mapped_reads[k])
+
+ self.failUnlessEqual(len(l1.match_codes), 9)
+ self.failUnlessEqual(len(l1.match_codes), len(l2.match_codes))
+ for k in l1.match_codes.keys():
+ self.failUnlessEqual(l1.match_codes[k],
+ l2.match_codes[k])
+
+ def test_runfolder(self):
+ runs = runfolder.get_runs(self.runfolder_dir)
+
+ # do we get the flowcell id from the filename?
+ self.failUnlessEqual(len(runs), 1)
+ self.failUnlessEqual(runs[0].name, 'run_207BTAAXX_2008-04-19.xml')
+
+ # do we get the flowcell id from the FlowcellId.xml file
+ make_flowcell_id(self.runfolder_dir, '207BTAAXY')
+ runs = runfolder.get_runs(self.runfolder_dir)
+ self.failUnlessEqual(len(runs), 1)
+ self.failUnlessEqual(runs[0].name, 'run_207BTAAXY_2008-04-19.xml')
+
+ r1 = runs[0]
+ xml = r1.get_elements()
+ xml_str = ElementTree.tostring(xml)
+
+ r2 = runfolder.PipelineRun(xml=xml)
+ self.failUnlessEqual(r1.name, r2.name)
+ self.failIfEqual(r2.firecrest, None)
+ self.failIfEqual(r2.bustard, None)
+ self.failIfEqual(r2.gerald, None)
+
+
+def suite():
+ return unittest.makeSuite(RunfolderTests,'test')
+
+if __name__ == "__main__":
+ unittest.main(defaultTest="suite")
+
+++ /dev/null
-"""
-Core information needed to inspect a runfolder.
-"""
-from glob import glob
-import logging
-import os
-import re
-import shutil
-import stat
-import subprocess
-import sys
-import time
-
-try:
- from xml.etree import ElementTree
-except ImportError, e:
- from elementtree import ElementTree
-
-EUROPEAN_STRPTIME = "%d-%m-%Y"
-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
-
-from gaworkflow.util.alphanum import alphanum
-from gaworkflow.util.ethelp import indent, flatten
-
-
-class PipelineRun(object):
- """
- Capture "interesting" information about a pipeline run
- """
- XML_VERSION = 1
- PIPELINE_RUN = 'PipelineRun'
- FLOWCELL_ID = 'FlowcellID'
-
- def __init__(self, pathname=None, firecrest=None, bustard=None, gerald=None, xml=None):
- if pathname is not None:
- self.pathname = os.path.normpath(pathname)
- else:
- self.pathname = None
- self._name = None
- self._flowcell_id = None
- self.firecrest = firecrest
- self.bustard = bustard
- self.gerald = gerald
-
- if xml is not None:
- self.set_elements(xml)
-
- def _get_flowcell_id(self):
- # extract flowcell ID
- if self._flowcell_id is None:
- config_dir = os.path.join(self.pathname, 'Config')
- flowcell_id_path = os.path.join(config_dir, 'FlowcellId.xml')
- if os.path.exists(flowcell_id_path):
- flowcell_id_tree = ElementTree.parse(flowcell_id_path)
- self._flowcell_id = flowcell_id_tree.findtext('Text')
- else:
- path_fields = self.pathname.split('_')
- if len(path_fields) > 0:
- # guessing last element of filename
- flowcell_id = path_fields[-1]
- else:
- flowcell_id = 'unknown'
-
- logging.warning(
- "Flowcell id was not found, guessing %s" % (
- flowcell_id))
- self._flowcell_id = flowcell_id
- return self._flowcell_id
- flowcell_id = property(_get_flowcell_id)
-
- def get_elements(self):
- """
- make one master xml file from all of our sub-components.
- """
- root = ElementTree.Element(PipelineRun.PIPELINE_RUN)
- flowcell = ElementTree.SubElement(root, PipelineRun.FLOWCELL_ID)
- flowcell.text = self.flowcell_id
- root.append(self.firecrest.get_elements())
- root.append(self.bustard.get_elements())
- root.append(self.gerald.get_elements())
- return root
-
- def set_elements(self, tree):
- # this file gets imported by all the others,
- # so we need to hide the imports to avoid a cyclic imports
- from gaworkflow.pipeline import firecrest
- from gaworkflow.pipeline import bustard
- from gaworkflow.pipeline import gerald
-
- tag = tree.tag.lower()
- if tag != PipelineRun.PIPELINE_RUN.lower():
- raise ValueError('Pipeline Run Expecting %s got %s' % (
- PipelineRun.PIPELINE_RUN, tag))
- for element in tree:
- tag = element.tag.lower()
- if tag == PipelineRun.FLOWCELL_ID.lower():
- self._flowcell_id = element.text
- #ok the xword.Xword.XWORD pattern for module.class.constant is lame
- elif tag == firecrest.Firecrest.FIRECREST.lower():
- self.firecrest = firecrest.Firecrest(xml=element)
- elif tag == bustard.Bustard.BUSTARD.lower():
- self.bustard = bustard.Bustard(xml=element)
- elif tag == gerald.Gerald.GERALD.lower():
- self.gerald = gerald.Gerald(xml=element)
- else:
- logging.warn('PipelineRun unrecognized tag %s' % (tag,))
-
- def _get_run_name(self):
- """
- Given a run tuple, find the latest date and use that as our name
- """
- if self._name is None:
- tmax = max(self.firecrest.time, self.bustard.time, self.gerald.time)
- timestamp = time.strftime('%Y-%m-%d', time.localtime(tmax))
- self._name = 'run_'+self.flowcell_id+"_"+timestamp+'.xml'
- return self._name
- name = property(_get_run_name)
-
- def save(self, destdir=None):
- if destdir is None:
- destdir = ''
- logging.info("Saving run report "+ self.name)
- xml = self.get_elements()
- indent(xml)
- dest_pathname = os.path.join(destdir, self.name)
- ElementTree.ElementTree(xml).write(dest_pathname)
-
- def load(self, filename):
- logging.info("Loading run report from " + filename)
- tree = ElementTree.parse(filename).getroot()
- self.set_elements(tree)
-
-def get_runs(runfolder):
- """
- Search through a run folder for all the various sub component runs
- and then return a PipelineRun for each different combination.
-
- For example if there are two different GERALD runs, this will
- generate two different PipelineRun objects, that differ
- in there gerald component.
- """
- from gaworkflow.pipeline import firecrest
- from gaworkflow.pipeline import bustard
- from gaworkflow.pipeline import gerald
-
- datadir = os.path.join(runfolder, 'Data')
-
- logging.info('Searching for runs in ' + datadir)
- runs = []
- for firecrest_pathname in glob(os.path.join(datadir,"*Firecrest*")):
- f = firecrest.firecrest(firecrest_pathname)
- bustard_glob = os.path.join(firecrest_pathname, "Bustard*")
- for bustard_pathname in glob(bustard_glob):
- b = bustard.bustard(bustard_pathname)
- gerald_glob = os.path.join(bustard_pathname, 'GERALD*')
- for gerald_pathname in glob(gerald_glob):
- try:
- g = gerald.gerald(gerald_pathname)
- runs.append(PipelineRun(runfolder, f, b, g))
- except IOError, e:
- print "Ignoring", str(e)
- return runs
-
-
-def extract_run_parameters(runs):
- """
- Search through runfolder_path for various runs and grab their parameters
- """
- for run in runs:
- run.save()
-
-def summarize_mapped_reads(mapped_reads):
- """
- Summarize per chromosome reads into a genome count
- But handle spike-in/contamination symlinks seperately.
- """
- summarized_reads = {}
- genome_reads = 0
- genome = 'unknown'
- for k, v in mapped_reads.items():
- path, k = os.path.split(k)
- if len(path) > 0:
- genome = path
- genome_reads += v
- else:
- summarized_reads[k] = summarized_reads.setdefault(k, 0) + v
- summarized_reads[genome] = genome_reads
- return summarized_reads
-
-def summary_report(runs):
- """
- Summarize cluster numbers and mapped read counts for a runfolder
- """
- report = []
- for run in 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.sort(alphanum)
-
- lane_results = run.gerald.summary.lane_results
- for lane_id in eland_keys:
- result = run.gerald.eland_results.results[lane_id]
- report.append("Sample name %s" % (result.sample_name))
- report.append("Lane id %s" % (result.lane_id,))
- cluster = lane_results[result.lane_id].cluster
- report.append("Clusters %d +/- %d" % (cluster[0], cluster[1]))
- report.append("Total Reads: %d" % (result.reads))
- mc = result._match_codes
- nm = mc['NM']
- nm_percent = float(nm)/result.reads * 100
- qc = mc['QC']
- qc_percent = float(qc)/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(result.mapped_reads)
- for name, counts in mapped_reads.items():
- report.append(" %s: %d" % (name, counts))
- report.append('---')
- report.append('')
- return os.linesep.join(report)
-
-def extract_results(runs, output_base_dir=None):
- if output_base_dir is None:
- output_base_dir = os.getcwd()
-
- for r in runs:
- result_dir = os.path.join(output_base_dir, r.flowcell_id)
- logging.info("Using %s as result directory" % (result_dir,))
- if not os.path.exists(result_dir):
- os.mkdir(result_dir)
-
- # create cycle_dir
- cycle = "C%d-%d" % (r.firecrest.start, r.firecrest.stop)
- logging.info("Filling in %s" % (cycle,))
- cycle_dir = os.path.join(result_dir, cycle)
- if os.path.exists(cycle_dir):
- logging.error("%s already exists, not overwriting" % (cycle_dir,))
- continue
- else:
- os.mkdir(cycle_dir)
-
- # copy stuff out of the main run
- g = r.gerald
-
- # save run file
- r.save(cycle_dir)
-
- # Copy Summary.htm
- summary_path = os.path.join(r.gerald.pathname, 'Summary.htm')
- if os.path.exists(summary_path):
- logging.info('Copying %s to %s' % (summary_path, cycle_dir))
- shutil.copy(summary_path, cycle_dir)
- else:
- logging.info('Summary file %s was not found' % (summary_path,))
-
- # tar score files
- score_files = []
- for f in os.listdir(g.pathname):
- if re.match('.*_score.txt', f):
- score_files.append(f)
-
- tar_cmd = ['/bin/tar', 'c'] + score_files
- 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("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)
- 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):
- """
- 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
--- /dev/null
+#!/usr/bin/python
+import optparse
+import sys
+import os
+
+from gaworkflow.util.makebed import make_bed_from_eland_stream, make_bed_from_multi_eland_stream, make_description
+
+def make_parser():
+ parser = optparse.OptionParser()
+ parser.add_option('-e', '--eland', dest='inname',
+ help='specify input eland filename')
+ parser.add_option('-b', '--bed', dest='outname',
+ help='specify output befilename')
+ parser.add_option('-n', '--name', dest='name',
+ help='specify the track (short) name.',
+ default=None)
+ parser.add_option('-d', '--description', dest='description',
+ help='specify the track description',
+ default=None)
+ parser.add_option('--chromosome', dest='prefix',
+ help='Set the chromosome prefix name. defaults to "chr"',
+ default='chr')
+ parser.add_option("--database", dest='database',
+ help="specify location of fctracker database",
+ default=None)
+ parser.add_option("--flowcell", dest='flowcell',
+ help="compute name and description from database using flowcell id",
+ default=None)
+ parser.add_option("--lane", dest='lane',
+ help='specify which lane to use when retrieving description from database',
+ default=None)
+
+ multi = optparse.OptionGroup(parser, 'Multi-read ELAND support')
+
+ multi.add_option('-m', '--multi', action='store_true',
+ help='Enable parsing multi-read eland files',
+ default=False)
+ multi.add_option('--reads', type='int',
+ help='limit reporting multi reads to this many reads'
+ '(most usefully --reads=1 will turn a multi-read '
+ 'file into a single read file)',
+ default=255)
+ parser.add_option_group(multi)
+
+ return parser
+
+def main(command_line=None):
+ if command_line is None:
+ command_line = sys.argv[1:]
+
+ parser = make_parser()
+ (options, args) = parser.parse_args(command_line)
+
+ if options.inname is None:
+ parser.error("Need eland input file name")
+ return 1
+
+ if options.inname == '-':
+ instream = sys.stdin
+ elif os.path.exists(options.inname):
+ instream = open(options.inname, 'r')
+ else:
+ parser.error('%s was not found' % (options.inname))
+ return 1
+
+ if options.outname is None:
+ # if outname wasn't defined, and we're reading from stdout
+ if instream is sys.stdin:
+ # write to stdout
+ outstream = sys.stdout
+ else:
+ # if there's a name write to name.bde
+ options.outname = os.path.splitext(options.inname)[0]+'.bed'
+ print >>sys.stderr, "defaulting to outputname", options.outname
+ elif options.outname == '-':
+ outstream = sys.stdout
+ elif os.path.exists(options.outname):
+ parser.error("not overwriting %s" % (options.outname))
+ return 1
+ else:
+ outstream = open(options.outname, 'w')
+
+ if options.flowcell is not None and options.lane is not None:
+ # get our name/description out of the database
+ name, description = make_description(
+ options.database, options.flowcell, options.lane
+ )
+ else:
+ name = options.name
+ description = options.description
+
+ if options.multi:
+ make_bed_from_multi_eland_stream(instream, outstream,
+ name, description,
+ options.prefix,
+ options.reads)
+
+ else:
+ make_bed_from_eland_stream(instream, outstream,
+ name, description,
+ options.prefix)
+ return 0
+
+if __name__ == "__main__":
+ sys.exit(main(sys.argv[1:]))
+
--- /dev/null
+#!/bin/sh
+
+if [ -z $HTSW_ROOT ]; then
+ HTSW_ROOT=$(pwd)
+fi
+
+PYTHONPATH=${HTSW_ROOT}/htswcommon:${HTSW_ROOT}/htswdataprod:${HTSW_ROOT}/htswfrontend:$PYTHONPATH
+
+nosetests -w ${HTSW_ROOT}/htswcommon
+nosetests -w ${HTSW_ROOT}/htswdataprod
projects / htsworkflow.git / commitdiff