From 125440bc659b5bb329d75c08c4430c4da877dfa1 Mon Sep 17 00:00:00 2001 From: Godwin Ani Date: Tue, 11 Mar 2025 10:11:32 -0400 Subject: [PATCH] Add files via upload --- PTL2/Scripts/guidance.py | 157 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100644 PTL2/Scripts/guidance.py diff --git a/PTL2/Scripts/guidance.py b/PTL2/Scripts/guidance.py new file mode 100644 index 0000000..6794a99 --- /dev/null +++ b/PTL2/Scripts/guidance.py @@ -0,0 +1,157 @@ +# 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