Merge pull request #96 from dib-lab/fix_ainput_err

fix assemblyinput for ppl specifying out_path

Author: bluegenes

ep_utils/pipeline_defaults.yaml

@@ -6,7 +6,7 @@ basename: eelpond
 # these outdirectories will be available to all workflows and rules
 eelpond_directories:
   rules: rules
-  animals: animals
+  animals: ep_utils/animals
   logs: logs
   outdirs:
     input_data: input_data

ep_utils/utils.py

@@ -72,6 +72,31 @@ def update_nested_dict(d, other):
 def is_single_end(sample, unit, end = '', assembly = ''):
     return pd.isnull(samples.loc[(sample, unit), "fq2"])
 
+def handle_assemblyinput(assembInput, config):
+    extensions= {}
+    program_params = config['assemblyinput'].get('program_params')
+    assemblyfile = program_params.get('assembly', None)
+    if assemblyfile:
+        assert os.path.exists(assemblyfile), 'Error: cannot find input assembly at {}\n'.format(assemblyfile)
+        sys.stderr.write('\tFound input assembly at {}\n'.format(assemblyfile))
+        assemblyfile = os.path.realpath(assemblyfile)
+    else:
+        sys.stderr.write("\n\tError: trying to run `assemblyinput` workflow, but there's no assembly file specified in your configfile. Please fix.\n\n")
+    gtmap = program_params.get('gene_trans_map', '')
+    if gtmap:
+        assert os.path.exists(gtmap), 'Error: cannot find assembly gene_trans_map at {}\n'.format(gtmap)
+        sys.stderr.write('\tFound input assembly gene-transcript map at {}\n'.format(gtmap))
+        extensions = {'base': ['.fasta', '.fasta.gene_trans_map']}
+        gtmap = os.path.realpath(gtmap)
+    else:
+        program_params['gene_trans_map'] = ''
+        config['no_gene_trans_map']= True
+    # grab user-input assembly extension
+    input_assembly_extension = program_params.get('assembly_extension', '_input')
+    extensions['assembly_extensions'] = [input_assembly_extension]
+    config['assemblyinput'] = {'program_params': program_params, 'eelpond_params': {'extensions':extensions}}
+    return config, input_assembly_extension
+
 def generate_targs(outdir, basename, samples, assembly_exts=[''], base_exts = None, read_exts = None, other_exts = None, contrasts = []):
     base_targets, read_targets, other_targs = [],[],[]
     # handle read targets

rules/deseq2/deseq2.rule

@@ -7,6 +7,10 @@ deseq_params = config['deseq2']['program_params']
 ep_params = config['deseq2']['eelpond_params']
 deseq_dir = config['deseq2']['eelpond_params']['outdir']
 
+# this is handled in `run_eelpond`, but redundancy doesn't hurt
+if config.get('no_gene_trans_map', False):
+    deseq_params['gene_trans_map'] = False
+
 def get_deseq2_threads(wildcards=None):
     # https://github.com/snakemake-workflows/rna-seq-star-deseq2
     # https://twitter.com/mikelove/status/918770188568363008

rules/utils/assemblyinput.rule

@@ -8,15 +8,26 @@ assembly_dir = ep_params['outdir']
 basename = config['basename']
 assembly_extension = aInput_params['assembly_extension']
 
+
 rule cp_input_assembly:
     input: aInput_params['assembly']
     output: join(assembly_dir, basename + assembly_extension + '.fasta')
     log: join(LOGS_DIR, 'assemblyinput/cp_assembly.log')
-    shell: ("cp {input} {output}")
+    params: assembly_dir = ep_params['outdir']
+    shell: 
+        """
+        mkdir -p {params.assembly_dir}
+        cp {input} {output}
+        """
 
 if aInput_params.get('gene_trans_map', None):
     rule cp_input_gene_trans_map:
         input: aInput_params['gene_trans_map'] 
         output: join(assembly_dir, basename + assembly_extension + '.fasta.gene_trans_map')
         log: join(LOGS_DIR, 'assemblyinput/cp_gt_map.log')
-        shell: ("cp {input} {output}")
+        params: assembly_dir = ep_params['outdir']
+        shell: 
+            """
+            mkdir -p {params.assembly_dir}
+            cp {input} {output}
+            """

run_eelpond

@@ -60,42 +60,42 @@ def build_default_params(workdir, targets):
     defaultParams['assembly_extensions'] = list(set(assembly_extensions))
     return defaultParams
 
+
 def build_dirs(ep_dir, params):
+    ''' function to build full paths for all directories '''
     # build main eelpond dir info
-    rules_dir = params['eelpond_directories']['rules']
     params['eelpond_directories']['base_dir'] = ep_dir
-    params['eelpond_directories']['rules'] = os.path.join(ep_dir, rules_dir)
-    animals_dir =  params['eelpond_directories']['animals']
-    params['eelpond_directories']['animals'] = os.path.join(ep_dir, 'ep_utils', animals_dir)
-    
+    params['eelpond_directories']['rules'] = os.path.join(ep_dir, params['eelpond_directories']['rules'])
+    params['eelpond_directories']['animals'] = os.path.join(ep_dir, params['eelpond_directories']['animals'])
     # if desired, user can also provide out_path, and all dirs will be built under there
-    # outdirectory, build all directories relative to that path 
     out_path = params.get('out_path', ep_dir)
-    # if user inputs an absolute path:
-    if os.path.isabs(out_path): # if not absolute, just assume subdirectory of eelpond.  
+    out_path = os.path.expanduser(out_path) # expand any `~` on unix
+    if os.path.isabs(out_path):  # if user inputs an absolute path, check that it exists!
         assert os.path.exists(out_path) and os.path.isdir(out_path), f"Error: provided output path {out_path} is not an existing directory. Please fix.\n\n"
-
+    else: # if not absolute, assume subdirectory of base eelpond dir
+        out_path = os.path.join(ep_dir, out_path)
     # allow user to define basename, and experiment, build outdir name
     basename = params['basename']
     expt = params.get('experiment', '')
     if expt and not expt.startswith('_'):
         expt= '_' + expt
     outdir = basename + "_out" + expt
-   
     # Now join out_path, outdir name
+    out_path = os.path.realpath(out_path)
     outdir = os.path.join(out_path, outdir)  
-    
-    params['eelpond_directories']['out_dir'] = outdir # outdir NAME
-    # if we want logs to be *within* individual outdirs, change logsdir here (else logs in eelpond dir)
-    params['eelpond_directories']['logs'] = join(outdir, os.path.basename(params['eelpond_directories']['logs']))
+    # when using out_path, need to manually build the outdir (snakemake will not automatically create it)
+    if not os.path.exists(outdir):
+        os.makedirs(outdir)
+    # add full path info to the config
+    params['eelpond_directories']['out_dir'] = outdir
+    params['eelpond_directories']['logs'] = join(outdir, params['eelpond_directories']['logs'])
     # build dirs for main eelpond output directories
     outDirs = params['eelpond_directories']['outdirs']
     for targ, outD in outDirs.items():
         outDirs[targ] = os.path.join(outdir, outD)
     # put joined paths back in params file
     params['eelpond_directories']['outdirs'] = outDirs
     # build dirs for included rules
-
     included_rules = params['include_rules']
     for rule in included_rules:
         prog = os.path.basename(rule).split('.rule')[0]
@@ -104,28 +104,6 @@ def build_dirs(ep_dir, params):
         params[prog]['eelpond_params']['outdir'] = os.path.join(outdir, prog_dir) 
     return params
 
-def handle_assemblyInput(assembInput, config):
-    # find files
-    program_params = config['assemblyinput'].get('program_params')
-    assemblyfile = program_params.get('assembly', None)
-    assert os.path.exists(assemblyfile), 'Error: cannot find input assembly at {}\n'.format(assemblyfile)
-    sys.stderr.write('\tFound input assembly at {}\n'.format(assemblyfile))
-    gtmap = program_params.get('gene_trans_map', '')
-    extensions= {}
-    if gtmap:
-        assert os.path.exists(gtmap), 'Error: cannot find assembly gene_trans_map at {}\n'.format(gtmap)
-        sys.stderr.write('\tFound input assembly gene-transcript map at {}\n'.format(gtmap))
-        extensions = {'base': ['.fasta', '.fasta.gene_trans_map']} # kinda hacky. do this better
-    else:
-        program_params['gene_trans_map'] = ''
-    # grab user-input assembly extension
-    input_assembly_extension = program_params.get('assembly_extension', '')
-    config['assemblyinput'] = {}
-    config['assemblyinput']['program_params'] = program_params 
-    extensions['assembly_extensions'] = [input_assembly_extension]
-    config['assemblyinput']['eelpond_params'] = {'extensions': extensions} 
-    return config, input_assembly_extension
-
 def main(args):
 
     thisdir = os.path.abspath(os.path.dirname(__file__))
@@ -188,7 +166,7 @@ def main(args):
         assembInput = configD.get('assemblyinput', None)
         if assembInput:
             targs+=['assemblyinput']
-            configD, assemblyinput_ext = handle_assemblyInput(assembInput, configD)
+            configD, assemblyinput_ext = handle_assemblyinput(assembInput, configD)
         else:
             assemblyinput_ext = None
         if 'assemblyinput' in targs and not assembInput:
@@ -225,6 +203,12 @@ def main(args):
         # add extension to overall assembly_extensions info
         if assemblyinput_ext: # note, need to do it here to prevent override with defaults
             paramsD['assembly_extensions'] = list(set(paramsD.get('assembly_extensions', []) + [assemblyinput_ext]))
+            
+            # NOTE: I think this should be moved from here into an `input_checks` function that checks all appropriate inputs/outputs are specified for the rules in use, and provides appropriate links to docs to help guide the user.
+            if paramsD.get('no_gene_trans_map', False):
+                if paramsD.get('deseq2'):
+                    paramsD['deseq2']['program_params']['gene_trans_map'] = False
+                    sys.stderr.write("\tYou're using `assemblyinput` without specifying a gene-trans-map. Setting differential expression to transcript-level only. See https://dib-lab.github.io/eelpond/deseq2/for details.\n")
 
         # use params to build directory structure
         paramsD = build_dirs(thisdir, paramsD)
@@ -247,6 +231,7 @@ def main(args):
             print('\tconfig: {}'.format(configfile))
             print('\tparams: {}'.format(paramsfile))
             print('\ttargets: {}'.format(repr(targs)))
+            print('\toutput: {}'.format(repr(paramsD['eelpond_directories']['out_dir'])))
             print('\treport: {}'.format(repr(reportfile)))
             print('--------')