Update Cluster_v2.0.py

Utilities/for_fastas/Cluster_v2.0.py

@@ -1,101 +1,66 @@
-#Professor L. Katz and Godwin Ani
+'''
-#9th- Feb- 2023
+#Author, date: Godwin Ani and Laura Katz, 9th- Feb - 2023.
-#A clustering program that accepts and validates users input
+#Dependencies: Python3, CD-Hit
+#Inputs: A folder of containing  Amino acid or DNA fasta files.
+#Outputs: A folder of clustered files.
+#Example: python Cluster_v2.0.py --type dna --identity 0.95 --overlap 0.67 --input input_folder_dna --output output_folder_dna
+'''
 
-import os, sys, re
+
+
+import os
+import argparse
 from tqdm import tqdm
+import subprocess
 
-
+def input_validation(value, error_message):
-def cluster():
-    ''' This function takes an amino acid or DNA sequence stored in a fasta format and clusters it.
-It uses two nested functions (cluster_AA and cluster_DNA) to perform this operation.
-Ensure that fasta files to be clustered are stored in a root directory folder named "to_cluster". This folder will automatically be created if none exists.
-The result of the clustering process is saved into a root directory folder named "Clustered" which is created automatically if none exists. '''
-  
-    if not os.path.isdir('to_cluster'):
-        os.mkdir('to_cluster')
-        
-    if not os.path.isdir('Clustered'):
-        os.mkdir('Clustered')
-    # Nested function for amino acids clustering.    
-    def cluster_AA():
-    #input validation for the sequence identity threshold.
-        while True:
     try:
-                val1 =input( 'Amino Acids Sequence Identity Threshold (e.g. 0.99, 0.95)? : ')
+        integer, fractional = value.split('.')
-                integer, fractional = val1.split('.')
+        value = float(value)
-                val1 = float(val1)
+        if int(integer) == 0 and len(fractional) == 2:
-                if int(integer)== 0 and len(fractional) == 2:
+            return value
-                    break
     except ValueError:
         pass
-            print('ERROR! Use format 0.## for Amino acids sequence identity threshold.')
+    print(error_message)
+    exit(1)
 
-      #Input validation for the overlap value.
+def cluster_sequences(program, threshold, overlap, input_folder, output_folder):
-        while True:
+    for file in tqdm(os.listdir(input_folder)):
-            try:
-                val2 =input( 'Amino Acids Alignment Overlap Value (e.g. 0.67, 0.75)? : ')
-                integer, fractional = val2.split('.')
-                val2 = float(val2)
-                if int(integer)== 0 and len(fractional) == 2:
-                    break
-            except ValueError:
-                pass
-            print('ERROR! Use format 0.## for Amino acids sequence alignment overlap value')
-            
-  #Selects amino acids fasta files in "to_cluster" folder and clusters them with CD-HIT.
-        for file in tqdm(os.listdir('to_cluster')):
         if file.endswith('.fasta'):
-                os.system('cd-hit -i to_cluster/' + file + ' -o Clustered/' + file + ' -c %s -d 0 -aS %s' %( val1, val2))
+            subprocess.run([f'{program}', '-i', f'{input_folder}/{file}', '-o', f'{output_folder}/{file}', '-c', f'{threshold}', '-d', '0', '-aS', f'{overlap}'])
-     #Renaming the result of the clustering.
+
-        for file in os.listdir('Clustered'):
+    for file in os.listdir(output_folder):
         if file.endswith('.clstr'):
-                os.rename('Clustered/' + file, 'Clustered/' + file.split('FILE')[0] + 'Clustered.txt')
+            os.rename(f'{output_folder}/{file}', f'{output_folder}/{file.split("FILE")[0]}Clustered.txt')
 
+def main():
+    parser = argparse.ArgumentParser(description='Cluster amino acid or DNA sequences using CD-HIT.')
+    parser.add_argument('--type', choices=['aa', 'dna'], required=True, help='Type of sequences (aa for Amino Acids, dna for DNA)')
+    parser.add_argument('--identity', type=str, required=True, help='Sequence Identity Threshold (e.g., 0.99, 0.95)')
+    parser.add_argument('--overlap', type=str, required=True, help='Sequence Alignment Overlap Value (e.g., 0.67, 0.75)')
+    parser.add_argument('--input', type=str, required=True, help='Input folder containing sequences in fasta format')
+    parser.add_argument('--output', type=str, required=True, help='Output folder for clustered sequences')
 
-    #Nested function for DNA clustering.
+    args = parser.parse_args()
-    def cluster_DNA():
-        #Input validation for the sequence identity threshold.
-        while True:
-            try:
-                val1 =input( 'DNA Sequence Identity Threshold (e.g. 0.99, 0.95)? : ')
-                integer, fractional = val1.split('.')
-                val1 = float(val1)
-                if int(integer)== 0 and len(fractional) == 2:
-                    break
-            except ValueError:
-                pass
-            print('ERROR! Use format 0.## for DNA sequence identity threshold.')
-      #Input validation for the overlap value. 
-        while True:
-            try:
-                val2 =input( 'DNA Sequence Alignment Overlap Value (e.g. 0.67, 0.75)? : ')
-                integer, fractional = val2.split('.')
-                val2 = float(val2)
-                if int(integer)== 0 and len(fractional) == 2:
-                    break
-            except ValueError:
-                pass
-            print('ERROR! Use format 0.## for DNA sequence alignment overlap value.')
 
-            #Selects DNA fasta files in "to_cluster" folder and clusters them with CD-HIT.
+    if not os.path.isdir(args.input):
-        for file in tqdm(os.listdir('to_cluster')):
+        print(f'Error: Input folder "{args.input}" does not exist.')
-            if file.endswith('.fasta'):
+        exit(1)
-                os.system('cd-hit-est -i to_cluster/' + file + ' -o Clustered/' + file + ' -c %s -d 0 -aS %s' %( val1, val2))
 
-        #Renaming the result of the clustering.
+    if not os.path.isdir(args.output):
-        for file in os.listdir('Clustered'):
+        os.mkdir(args.output)
-            if file.endswith('.clstr'):
-                os.rename('Clustered/' + file, 'Clustered/' + file.split('FILE')[0] + 'Clustered.txt')
 
-    # Prompts for user input and function call.
+    if args.type == 'aa':
-    choice_1 = input('Are you clustering Amino Acids or DNA? (AA or DNA) : ')
+        threshold = input_validation(args.threshold, 'ERROR! Use format 0.## for Amino acids sequence identity threshold.')
-    if choice_1 in ['AA', 'Aa', 'aa']:
+        overlap = input_validation(args.overlap, 'ERROR! Use format 0.## for Amino acids sequence alignment overlap value.')
-        cluster_AA()
+        cluster_sequences('cd-hit', threshold, overlap, args.input, args.output)
-    elif choice_1 in ['DNA', 'Dna', 'dna']:
+    elif args.type == 'dna':
-        cluster_DNA()
+        threshold = input_validation(args.threshold, 'ERROR! Use format 0.## for DNA sequence identity threshold.')
+        overlap = input_validation(args.overlap, 'ERROR! Use format 0.## for DNA sequence alignment overlap value.')
+        cluster_sequences('cd-hit-est', threshold, overlap, args.input, args.output)
     else:
-        print('Sorry. This program can only cluster Amino Acids and DNA')
+        print('Invalid sequence type. Choose "aa" for Amino Acids or "dna" for DNA.')
-
+        exit(1)
-cluster()
 
+if __name__ == "__main__":
+    main()