From 0a0433278442f348e375602bb3078bbb38f1bc82 Mon Sep 17 00:00:00 2001 From: Godwin Ani Date: Tue, 11 Mar 2025 13:23:07 -0400 Subject: [PATCH] Delete PTL2/Scripts/guidance.py --- PTL2/Scripts/guidance.py | 157 --------------------------------------- 1 file changed, 157 deletions(-) delete mode 100644 PTL2/Scripts/guidance.py diff --git a/PTL2/Scripts/guidance.py b/PTL2/Scripts/guidance.py deleted file mode 100644 index 6794a99..0000000 --- a/PTL2/Scripts/guidance.py +++ /dev/null @@ -1,157 +0,0 @@ -# Last updated Apr 2 2024 -# Authors: Auden Cote-L'Heureux and Mario Ceron-Romero - -# This script runs Guidance in an iterative fashion for more both MSA construction -# and more rigorous homology assessment than what is offered in EukPhylo part 1. -# Guidance runs until the input number of iterations (--guidance_iters, default = 5) -# has been reached, or until there are no sequences below the sequence score cutoff. -# All sequences below the score cutoff (--seq_cutoff, default = 0.3) are removed at -# each iteration. By default, EukPhylo does not remove residues that fall below the -# given residue cutoff (--res_cutoff) and columns that fall below the given column -# cutoff (--col_cutoff, defaults are 0), though this can be turned on by adjusting -# these parameters. Outputs at this point are found in the “Guidance_NotGapTrimmed” -# output folder. We then run MSAs through TrimAl to remove all sites in the alignment -# that are at least 95% gaps (or --gap_trim_cutoff) generating files in the “Guidance” -# output folder. - -# This step is either intended to be run starting with --start = unaligned (but not raw) -# inputs, meaning one amino acid alignment per OG. It can also be run directly after the -# preguidance step. The run() function is called in two places: in eukphylo.py generally, -# and in contamination.py if the contamination loop is using Guidance as the re-alignment -# method. - -#Dependencies -import os, sys, re -import subprocess -import logging -from Bio import SeqIO - -# Set up logging configuration -logging.basicConfig(filename='guidance_run.log', level=logging.DEBUG, - format='%(asctime)s - %(levelname)s - %(message)s') - -# Function to run system commands and capture logs -def run_command(command, log_file): - try: - logging.info(f"Running command: {command}") - result = subprocess.run(command, shell=True, check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True) - logging.info(f"Command output: {result.stdout}") - return result - except subprocess.CalledProcessError as e: - logging.error(f"Command failed with error: {e.stderr}") - with open(log_file, 'a') as log: - log.write(f"Command failed: {command}\nError: {e.stderr}\n") - return None - -# Called in eukphylo.py and contamination.py -def run(params): - if params.start == 'raw' or params.start == 'unaligned': - # Checking that pre-Guidance has been run or that unaligned files per OG are provided. - if params.start == 'raw': - preguidance_path = params.output + '/Output/Pre-Guidance' - else: - preguidance_path = params.data - - if not os.path.isdir(preguidance_path): - logging.error(f"ERROR: The path {preguidance_path} could not be found when trying to locate pre-Guidance (unaligned) files.") - exit() - - if len([f for f in os.listdir(preguidance_path) if f.endswith('.fa') or f.endswith('.faa') or f.endswith('.fasta')]) == 0: - logging.error(f"ERROR: No pre-Guidance (unaligned) files found at {preguidance_path}.") - exit() - - # Creating intermediate folders - os.mkdir(params.output + '/Output/Intermediate/Guidance') - os.mkdir(params.output + '/Output/Intermediate/Guidance/Input') - os.mkdir(params.output + '/Output/Intermediate/Guidance/Output') - - guidance_input = params.output + '/Output/Intermediate/Guidance/Input/' - subprocess.run(f'cp -r {preguidance_path}/* {guidance_input}', shell=True) - - guidance_removed_file = open(params.output + '/Output/GuidanceRemovedSeqs.txt', 'w') - guidance_removed_file.write('Sequence\tScore\n') - - too_many_seqs = False - - # Check if any files have more than 2000 sequences - for file in [f for f in os.listdir(guidance_input) if f.endswith('.fa') or f.endswith('.faa') or f.endswith('.fasta')]: - nseqs = len([rec for rec in SeqIO.parse(guidance_input + '/' + file, 'fasta')]) - if nseqs > 2000: - too_many_seqs = True - break - - if too_many_seqs and not params.allow_large_files: - return False - - # Process each unaligned AA fasta file - for file in [f for f in os.listdir(guidance_input) if f.endswith('.fa') or f.endswith('.faa') or f.endswith('.fasta')]: - tax_guidance_outdir = params.output + '/Output/Intermediate/Guidance/Output/' + file.split('.')[0].split('_preguidance')[0] - os.mkdir(tax_guidance_outdir) - - fail = False - # Iterative process - for i in range(params.guidance_iters): - n_recs = len([r for r in SeqIO.parse(guidance_input + '/' + file, 'fasta')]) - - if n_recs < 4: - logging.warning(f"Gene family {file.split('.')[0].split('_preguidance')[0]} contains fewer than 4 sequences after {i} iterations.") - subprocess.run(f'rm -rf {tax_guidance_outdir}', shell=True) - if i == 0: - fail = True - break - - # Determining MAFFT algorithm - mafft_alg = 'genafpair' if n_recs < 200 else 'auto' - - # Running Guidance - command = f"perl /guidance/guidance.v2.02/www/Guidance/guidance.pl --seqFile {guidance_input}/{file} --msaProgram MAFFT --seqType aa --outDir {tax_guidance_outdir} --seqCutoff {params.seq_cutoff} --colCutoff {params.col_cutoff} --outOrder as_input --bootstraps 10 --MSA_Param '\\--{mafft_alg} --maxiterate 1000 --thread {params.guidance_threads} --bl 62 --anysymbol' > {params.output}/Output/Intermediate/Guidance/Output/{file[:10]}/log.txt" - result = run_command(command, f'{params.output}/Output/Intermediate/Guidance/Output/{file[:10]}/log.txt') - - if result is None: - fail = True - break - - if os.path.isfile(tax_guidance_outdir + '/MSA.MAFFT.Guidance2_res_pair_seq.scr_with_Names'): - seqs_below = len([line for line in open(tax_guidance_outdir + '/MSA.MAFFT.Guidance2_res_pair_seq.scr_with_Names').readlines()[1:-1] if float(line.split()[-1]) < params.seq_cutoff]) - - if n_recs - seqs_below < 4: - logging.warning(f"Gene family {file.split('.')[0].split('_preguidance')[0]} contains fewer than 4 sequences after {i + 1} iterations.") - subprocess.run(f'rm -rf {tax_guidance_outdir}', shell=True) - break - - if seqs_below == 0 or i == params.guidance_iters - 1: - logging.info(f"Guidance complete after {i + 1} iterations for gene family {file.split('.')[0].split('_preguidance')[0]}") - break - - for line in [line for line in open(tax_guidance_outdir + '/MSA.MAFFT.Guidance2_res_pair_seq.scr_with_Names').readlines()[1:-1] if float(line.split()[-1]) < params.seq_cutoff]: - guidance_removed_file.write(line) - - subprocess.run(f'cp {tax_guidance_outdir}/Seqs.Orig.fas.FIXED.Without_low_SP_Seq.With_Names {guidance_input}/{file}', shell=True) - - # Handling intermediate files - if params.keep_iter: - iteration_folder = params.output + '/Output/Intermediate/Guidance/Iterations/' + str(i + 1) + '/' + file.split('.')[0].split('_preguidance')[0] - os.makedirs(iteration_folder, exist_ok=True) - subprocess.run(f'cp -r {tax_guidance_outdir}/* {iteration_folder}', shell=True) - - if not params.keep_temp: - for gdir_file in os.listdir(iteration_folder): - if gdir_file not in ('MSA.MAFFT.Guidance2_res_pair_seq.scr_with_Names', 'MSA.MAFFT.aln.With_Names', 'MSA.MAFFT.Guidance2_res_pair_col.scr', 'log', 'postGuidance_preTrimAl_unaligned.fasta'): - subprocess.run(f'rm -r {iteration_folder}/{gdir_file}', shell=True) - else: - if gdir_file == 'MSA.MAFFT.aln.With_Names': - subprocess.run(f'mv {iteration_folder}/{gdir_file} {iteration_folder}/{file.split(".")[0].split("_preguidance")[0]}_{gdir_file}.aln', shell=True) - else: - subprocess.run(f'mv {iteration_folder}/{gdir_file} {iteration_folder}/{file.split(".")[0].split("_preguidance")[0]}_{gdir_file}', shell=True) - - subprocess.run(f'rm -r {tax_guidance_outdir}/*', shell=True) - else: - fail = True - break - - if not fail: - logging.info(f"Residue and column cutoffs applied for {file.split('.')[0].split('_preguidance')[0]}") - # Additional steps to apply residue/column cutoffs and finalize the alignment can go here - - guidance_removed_file.close() - return True