Update single vs paired guessing code
[htsworkflow.git] / extra / ucsc_encode_submission / ucsc_gather.py
index 4cea711b96bafa3233361c3cbe44310f572a6bff..38d8deefc7a8a1b2dc07b3f445534d24fc4fef0d 100755 (executable)
@@ -4,13 +4,16 @@ import fnmatch
 from glob import glob
 import json
 import logging
+import netrc
 from optparse import OptionParser
 import os
 from pprint import pprint, pformat
 import shlex
 from StringIO import StringIO
-import time
+import stat
+from subprocess import Popen, PIPE
 import sys
+import time
 import types
 import urllib
 import urllib2
@@ -34,10 +37,10 @@ def main(cmdline=None):
     else:
         logging.basicConfig(level = logging.WARNING )        
     
-    apidata = {'apiid': opts.apiid, 'apikey': opts.apikey }
+    apidata = api.make_auth_from_opts(opts, parser)
 
-    if opts.host is None or opts.apiid is None or opts.apikey is None:
-        parser.error("Please specify host url, apiid, apikey")
+    if opts.makeddf and opts.daf is None:
+        parser.error("Please specify your daf when making ddf files")
 
     if len(args) == 0:
         parser.error("I need at least one library submission-dir input file")
@@ -67,22 +70,6 @@ def main(cmdline=None):
 
 
 def make_parser():
-    # Load defaults from the config files
-    config = SafeConfigParser()
-    config.read([os.path.expanduser('~/.htsworkflow.ini'), '/etc/htsworkflow.ini'])
-    
-    sequence_archive = None
-    apiid = None
-    apikey = None
-    apihost = None
-    SECTION = 'sequence_archive'
-    if config.has_section(SECTION):
-        sequence_archive = config.get(SECTION, 'sequence_archive',sequence_archive)
-        sequence_archive = os.path.expanduser(sequence_archive)
-        apiid = config.get(SECTION, 'apiid', apiid)
-        apikey = config.get(SECTION, 'apikey', apikey)
-        apihost = config.get(SECTION, 'host', apihost)
-
     parser = OptionParser()
 
     # commands
@@ -102,20 +89,14 @@ def make_parser():
     parser.add_option('--force', default=False, action="store_true",
                       help="Force regenerating fastqs")
 
-    # configuration options
-    parser.add_option('--apiid', default=apiid, help="Specify API ID")
-    parser.add_option('--apikey', default=apikey, help="Specify API KEY")
-    parser.add_option('--host',  default=apihost,
-                      help="specify HTSWorkflow host",)
-    parser.add_option('--sequence', default=sequence_archive,
-                      help="sequence repository")
-
     # debugging
     parser.add_option('--verbose', default=False, action="store_true",
                       help='verbose logging')
     parser.add_option('--debug', default=False, action="store_true",
                       help='debug logging')
 
+    api.add_auth_options(parser)
+    
     return parser
 
 
@@ -271,8 +252,12 @@ def link_daf(daf_path, library_result_map):
     base_daf = os.path.basename(daf_path)
     
     for lib_id, result_dir in library_result_map:
+        if not os.path.exists(result_dir):
+            raise RuntimeError("Couldn't find target directory %s" %(result_dir,))
         submission_daf = os.path.join(result_dir, base_daf)
         if not os.path.exists(submission_daf):
+            if not os.path.exists(daf_path):
+                raise RuntimeError("Couldn't find daf: %s" %(daf_path,))
             os.link(daf_path, submission_daf)
 
 
@@ -284,7 +269,7 @@ def make_submission_ini(host, apidata, library_result_map, paired=True):
 
     for lib_id, result_dir in library_result_map:
         order_by = ['order_by=files', 'view', 'replicate', 'cell', 
-                    'readType', 'mapAlgorithm', 'insertLength' ]
+                    'readType', 'mapAlgorithm', 'insertLength', 'md5sum' ]
         inifile =  ['[config]']
         inifile += [" ".join(order_by)]
         inifile += ['']
@@ -299,6 +284,7 @@ def make_submission_ini(host, apidata, library_result_map, paired=True):
             attributes = view_map.find_attributes(f, lib_id)
             if attributes is None:
                 raise ValueError("Unrecognized file: %s" % (f,))
+            attributes['md5sum'] = "None"
             
             ext = attributes["extension"]
             if attributes['view'] is None:                   
@@ -307,6 +293,9 @@ def make_submission_ini(host, apidata, library_result_map, paired=True):
                 fastqs.setdefault(ext, set()).add(f)
                 fastq_attributes[ext] = attributes
             else:
+                md5sum = make_md5sum(os.path.join(result_dir,f))
+                if md5sum is not None:
+                    attributes['md5sum']=md5sum
                 inifile.extend(
                     make_submission_section(line_counter,
                                             [f],
@@ -378,6 +367,9 @@ def make_ddf(ininame,  daf_name, guess_ddf=False, make_condor=False, outdir=None
         output = open(ddf_name,'w')
 
     file_list = read_ddf_ini(ininame, output)
+    logging.info(
+        "Read config {0}, found files: {1}".format(
+            ininame, ", ".join(file_list)))
 
     file_list.append(daf_name)
     if ddf_name is not None:
@@ -480,7 +472,7 @@ def make_condor_upload_script(ininame):
     script = """Universe = vanilla
 
 Executable = /usr/bin/lftp
-arguments = -c put ../%(archivename)s -o ftp://detrout@encodeftp.cse.ucsc.edu/%(archivename)s
+arguments = -c put ../%(archivename)s -o ftp://%(ftpuser)s:%(ftppassword)s@%(ftphost)s/%(archivename)s
 
 Error = upload.err.$(Process).log
 Output = upload.out.$(Process).log
@@ -489,14 +481,24 @@ initialdir = %(initialdir)s
 
 queue 
 """
+    auth = netrc.netrc(os.path.expanduser("~diane/.netrc"))
+    
+    encodeftp = 'encodeftp.cse.ucsc.edu'
+    ftpuser = auth.hosts[encodeftp][0]
+    ftppassword = auth.hosts[encodeftp][2]
     context = {'archivename': make_submission_name(ininame),
                'initialdir': os.getcwd(),
-               'user': os.getlogin()}
+               'user': os.getlogin(),
+               'ftpuser': ftpuser,
+               'ftppassword': ftppassword,
+               'ftphost': encodeftp}
 
     condor_script = make_condor_name(ininame, 'upload')
     condor_stream = open(condor_script,'w')
     condor_stream.write(script % context)
     condor_stream.close()
+    os.chmod(condor_script, stat.S_IREAD|stat.S_IWRITE)
+
     return condor_script
 
 
@@ -623,26 +625,50 @@ class NameToViewMap(object):
         ma = 'TH1014'
 
         self.patterns = [
+            # for 2011 Feb 18 elements submission
+            ('final_Cufflinks_genes_*gtf',       'GeneDeNovo'),
+            ('final_Cufflinks_transcripts_*gtf', 'TranscriptDeNovo'),
+            ('final_exonFPKM-Cufflinks-0.9.3-GENCODE-v3c-*.gtf',       
+             'ExonsGencV3c'),
+            ('final_GENCODE-v3-Cufflinks-0.9.3.genes-*gtf',          
+             'GeneGencV3c'),
+            ('final_GENCODE-v3-Cufflinks-0.9.3.transcripts-*gtf',    
+             'TranscriptGencV3c'),
+            ('final_TSS-Cufflinks-0.9.3-GENCODE-v3c-*.gtf', 'TSS'),
+            ('final_junctions-*.bed6+3',                    'Junctions'),
+            
             ('*.bai',                   None),
             ('*.splices.bam',           'Splices'),
             ('*.bam',                   self._guess_bam_view),
             ('junctions.bed',           'Junctions'),
             ('*.jnct',                  'Junctions'),
-            ('*.plus.bigwig',           'PlusSignal'),
-            ('*.minus.bigwig',          'MinusSignal'),
+            ('*unique.bigwig',         None),
+            ('*plus.bigwig',           'PlusSignal'),
+            ('*minus.bigwig',          'MinusSignal'),
             ('*.bigwig',                'Signal'),
             ('*.tar.bz2',               None),
             ('*.condor',                None),
             ('*.daf',                   None),
             ('*.ddf',                   None),
+
+            ('*ufflinks?0.9.3.genes.gtf',       'GeneDeNovo'),
+            ('*ufflinks?0.9.3.transcripts.gtf', 'TranscriptDeNovo'),
+            ('*GENCODE-v3c.exonFPKM.gtf',        'ExonsGencV3c'),
+            ('*GENCODE-v3c.genes.gtf',           'GeneGencV3c'),
+            ('*GENCODE-v3c.transcripts.gtf',     'TranscriptGencV3c'),
+            ('*GENCODE-v3c.TSS.gtf',             'TSS'),
+            ('*.junctions.bed6+3',                'Junctions'),
+            
             ('*.?ufflinks-0.9.0?genes.expr',       'GeneDeNovo'),
             ('*.?ufflinks-0.9.0?transcripts.expr', 'TranscriptDeNovo'),
             ('*.?ufflinks-0.9.0?transcripts.gtf',  'GeneModel'),
-            ('*.GENCODE-v3c?genes.expr',       'GeneGencV3c'),
-            ('*.GENCODE-v3c?transcripts.expr', 'TranscriptGencV3c'),
-            ('*.GENCODE-v4?genes.expr',       'GeneGencV4'),
-            ('*.GENCODE-v4?transcripts.expr', 'TranscriptGencV4'),
-            ('*.GENCODE-v4?transcript.expr', 'TranscriptGencV4'),
+
+            ('*.GENCODE-v3c?genes.expr',       'GeneGCV3c'),
+            ('*.GENCODE-v3c?transcript*.expr', 'TranscriptGCV3c'),
+            ('*.GENCODE-v3c?transcript*.gtf',  'TranscriptGencV3c'),
+            ('*.GENCODE-v4?genes.expr',        None), #'GeneGCV4'),
+            ('*.GENCODE-v4?transcript*.expr',  None), #'TranscriptGCV4'),
+            ('*.GENCODE-v4?transcript*.gtf',   None), #'TranscriptGencV4'),
             ('*_1.75mers.fastq',              'FastqRd1'),
             ('*_2.75mers.fastq',              'FastqRd2'),
             ('*_r1.fastq',              'FastqRd1'),
@@ -651,11 +677,13 @@ class NameToViewMap(object):
             ('*.gtf',                   'GeneModel'),
             ('*.ini',                   None),
             ('*.log',                   None),
+            ('*.md5',                   None),
             ('paired-end-distribution*', 'InsLength'),
             ('*.stats.txt',              'InsLength'),
             ('*.srf',                   None),
             ('*.wig',                   None),
             ('*.zip',                   None),
+            ('transfer_log',            None),
             ]
 
         self.views = {
@@ -671,14 +699,17 @@ class NameToViewMap(object):
             "GeneModel": {"MapAlgorithm": ma},
             "GeneDeNovo": {"MapAlgorithm": ma},
             "TranscriptDeNovo": {"MapAlgorithm": ma},
+            "ExonsGencV3c": {"MapAlgorithm": ma},
             "GeneGencV3c": {"MapAlgorithm": ma},
+            "TSS": {"MapAlgorithm": ma},
+            "GeneGCV3c": {"MapAlgorithm": ma},
+            "TranscriptGCV3c": {"MapAlgorithm": ma},
             "TranscriptGencV3c": {"MapAlgorithm": ma},
-            "GeneGencV4": {"MapAlgorithm": ma},
-            "TranscriptGencV4": {"MapAlgorithm": ma},
+            "GeneGCV4": {"MapAlgorithm": ma},
+            "TranscriptGCV4": {"MapAlgorithm": ma},
             "FastqRd1": {"MapAlgorithm": "NA", "type": "fastq"},
             "FastqRd2": {"MapAlgorithm": "NA", "type": "fastq"},
             "Fastq": {"MapAlgorithm": "NA", "type": "fastq" },
-            "GeneModel": {"MapAlgorithm": ma},
             "InsLength": {"MapAlgorithm": ma},
             }
         # view name is one of the attributes
@@ -732,8 +763,14 @@ class NameToViewMap(object):
 
     def _is_paired(self, lib_id, lib_info):
         """Determine if a library is paired end"""
+        # TODO: encode this information in the library type page.
+        single = (1,3,6)
         if len(lib_info["lane_set"]) == 0:
-            return False
+            # we haven't sequenced anything so guess based on library type
+            if lib_info['library_type_id'] in single:
+                return False
+            else:
+                return True
 
         if not self.lib_paired.has_key(lib_id):
             is_paired = 0
@@ -856,6 +893,46 @@ def validate_filelist(files):
         if not os.path.exists(f):
             raise RuntimeError("%s does not exist" % (f,))
 
+def make_md5sum(filename):
+    """Quickly find the md5sum of a file
+    """
+    md5_cache = os.path.join(filename+".md5")
+    print md5_cache
+    if os.path.exists(md5_cache):
+        logging.debug("Found md5sum in {0}".format(md5_cache))
+        stream = open(md5_cache,'r')
+        lines = stream.readlines()
+        md5sum = parse_md5sum_line(lines, filename)
+    else:
+        md5sum = make_md5sum_unix(filename, md5_cache)
+    return md5sum
+    
+def make_md5sum_unix(filename, md5_cache):
+    cmd = ["md5sum", filename]
+    logging.debug("Running {0}".format(" ".join(cmd)))
+    p = Popen(cmd, stdout=PIPE)
+    stdin, stdout = p.communicate()
+    retcode = p.wait()
+    logging.debug("Finished {0} retcode {1}".format(" ".join(cmd), retcode))
+    if retcode != 0:
+        logging.error("Trouble with md5sum for {0}".format(filename))
+        return None
+    lines = stdin.split(os.linesep)
+    md5sum = parse_md5sum_line(lines, filename)
+    if md5sum is not None:
+        logging.debug("Caching sum in {0}".format(md5_cache))
+        stream = open(md5_cache, "w")
+        stream.write(stdin)
+        stream.close()
+    return md5sum
+
+def parse_md5sum_line(lines, filename):
+    md5sum, md5sum_filename = lines[0].split()
+    if md5sum_filename != filename:
+        errmsg = "MD5sum and I disagre about filename. {0} != {1}"
+        logging.error(errmsg.format(filename, md5sum_filename))
+        return None
+    return md5sum
 
 if __name__ == "__main__":
     main()