Major refactoring of auxiliary functions (#10)

* refactor map_habitat.py

delete unnecessary code repetition
fix readability
makes code cleaner

* refactor filter_length.py

minor changes, improve code readability

* fix problem with merging step

* fix format of the table

previously, there was a column separated by multiple spaces
this occasionated one of the values as NA, which breaks the
merge procedure by PANDAS.

* refactor and fix column renaming

column names were wrongly associated before to float valued columns
now it is fixed by always fixing the position of these columns as the first 
two.

* fix variable naming

* refactor map_quality.py

fix some internal functions, reduce overworking
improved readability

* fix error in the loop to include only small orfs

The loop before would include all sequences regardless the result
of the tests before (length and if was already seen), so now it loops
first in the length, if condition met, then goes to the seen test, and if
not seen before, then it is added to our dictionary, avoiding overwritting
and conflicts.

* refactor and enhance the deeplca function

Improved readability of the code, and reduce the complexity and 
time to compute LCA for each smorf.

* Update map_taxonomy.py

* Update map_taxonomy.py

* Update map_taxonomy.py

* fix error in format of one variable

Author: celiosantosjr

gmsc_mapper/filter_length.py

@@ -1,22 +1,40 @@
-def filter_length(queryfile,tmpdirname,N):
+message_error_all_long = '''GMSC-mapper Error: Input sequences are all more than 303nt or 100aa,
+which will be filtered. There are some options:
+1.If you don't want to filter, please use --nofilter flag.
+2.Please check if your input consist of contigs, which should
+use -i not --nt-genes or --aa-genes as input.
+'''
+
+message_error_longer='''GMSC-mapper Warning: Input has seqences more than 303nt or 100aa,
+which will be filtered. If you don't want to filter,
+please use --nofilter flag.
+'''
+
+
+def filter_length(queryfile, tmpdirname, N):
     import sys
+    
     from os import path
     from .fasta import fasta_iter
-    filtered_file = path.join(tmpdirname,"filtered.faa")
+    
+    filtered_file = path.join(tmpdirname, "filtered.faa")
 
-    with open(filtered_file,'wt') as of:
+    with open(filtered_file, 'wt') as of:
         all_longer_flag = 1
         longer_exist_flag = 0
-        for ID,seq in fasta_iter(queryfile):
+        
+        for ID, seq in fasta_iter(queryfile):
             if len(seq) < N:
                 all_longer_flag = 0
                 of.write(f'>{ID}\n{seq}\n')
             else:
                 longer_exist_flag = 1
+        
         if all_longer_flag:
-            sys.stderr.write("GMSC-mapper Error: Input sequences are all more than 303nt or 100aa,which will be filtered. 1.If you don't want to filter,please use --nofilter flag. 2.Please check if your input is contigs,which should use -i not --nt-genes or --aa-genes as input.\n")
+            sys.stderr.write(message_error_all_long)
             sys.exit(1)
         else:
             if longer_exist_flag:
-                print("GMSC-mapper Warning: Input has seqences more than 303nt or 100aa,which will be filtered.If you don't want to filter,please use --nofilter flag.\n")
-    return filtered_file
+                print(message_error_longer)
+    
+    return filtered_file

gmsc_mapper/map_habitat.py

@@ -1,39 +1,75 @@
 import pandas as pd
+
 from os import path
 
-def smorf_habitat(outdir,habitatfile,resultfile):
-    habitat_file = path.join(outdir,"habitat.out.smorfs.tsv")	
 
-    result = pd.read_csv(resultfile,sep='\t',header=None)
-    result = result.rename(columns={0:'qseqid',1:'sseqid'})
-    if habitatfile.endswith('.gz'):
-        reader =  pd.read_csv(habitatfile,compression="gzip",sep="\t",chunksize=5000000,header=None)
-    if habitatfile.endswith('.xz'):
-        reader =  pd.read_csv(habitatfile,compression="xz",sep="\t",chunksize=5000000,header=None)
-    if habitatfile.endswith('.bz2'):
-        reader =  pd.read_csv(habitatfile,compression="bz2",sep="\t",chunksize=5000000,header=None)
-    else:
-        reader =  pd.read_csv(habitatfile,sep="\t",chunksize=5000000,header=None)
+def fixdf(x):
+    x = x.dropna()
+    x = x.drop_duplicates()
+    return ','.join(x)
+    
+    
+def formatlabel(x):
+    x = x.split(',')
+    x = list(set(x))
+    x = sorted(x)
+    return ','.join(x)
+    
+        
+def smorf_habitat(outdir, habitatfile, resultfile):
+    habitat_file = path.join(outdir, "habitat.out.smorfs.tsv")	
+
+    result = pd.read_csv(resultfile,
+                         sep='\t',
+                         header=None)
+                         
+    result.rename({0: 'qseqid', 1: 'sseqid'},
+                  axis=1,
+                  inplace=True)
+                         
+    reader =  pd.read_table(habitatfile,
+                            sep="\t",
+                            chunksize=5_000_000,
+                            header=None,
+                            names=['sseqid', 'habitat'])
 
     output_list = []
     for chunk in reader:
-        chunk.columns = ['sseqid','habitat']
-        output_chunk = pd.merge(result,chunk,how='left')[['qseqid', 'habitat']]
+        output_chunk = result.merge(on='sseqid',
+                                    right=chunk,
+                                    how='left')
+        output_chunk = output_chunk[['qseqid', 'habitat']]
         output_list.append(output_chunk)
-    output = pd.concat(output_list, axis=0).sort_values(by='qseqid')
-    output = output.groupby('qseqid',as_index=False,sort=False).agg({'habitat':lambda x : ','.join(x.dropna().drop_duplicates())})
-    output['habitat'] = output['habitat'].apply(lambda x: ','.join(sorted(list(set(x.split(','))))))
-    output.to_csv(habitat_file,sep='\t',index=False)
+        
+    output = pd.concat(output_list,
+                       axis=0)
+    
+    output = output.sort_values(by='qseqid')
+    
+    output = output.groupby('qseqid',
+                            as_index=False,
+                            sort=False)
+    
+    output = output.agg({'habitat':lambda x : fixdf(x)})
+    output['habitat'] = output['habitat'].apply(lambda x: formatlabel(x))
+    
+    output.to_csv(habitat_file,
+                  sep='\t',
+                  index=False)
+
+    wdf = output['habitat'].apply(lambda x: len(x.split(',')))
+    number_dict = dict(wdf.value_counts())
+    number_dict_normalize = dict(wdf.value_counts(normalize=True))
 
-    number_dict = dict(output['habitat'].apply(lambda x: len(x.split(','))).value_counts())
-    number_dict_normalize = dict(output['habitat'].apply(lambda x: len(x.split(','))).value_counts(normalize=True))
     if 1 in number_dict.keys():
         single_number = number_dict[1]
         single_percentage = number_dict_normalize[1]
     else:
         single_number = 0
         single_percentage = 0
+        
     multi_number = output['habitat'].size - single_number
     multi_percentage = 1 - single_percentage
 
-    return single_number,single_percentage,multi_number,multi_percentage
+    return (single_number, single_percentage, multi_number, multi_percentage)
+    

gmsc_mapper/map_quality.py

@@ -1,35 +1,60 @@
-from os import path
 import pandas as pd
 
-def smorf_quality(outdir,qualityfile,resultfile):
-    quality_file = path.join(outdir,"quality.out.smorfs.tsv")	
+from os import path
+
 
-    result = pd.read_csv(resultfile,sep='\t',header=None)
-    result = result.rename(columns={0:'qseqid',1:'sseqid'})
-    if qualityfile.endswith('.gz'):
-        ref_quality =  pd.read_csv(qualityfile,compression="gzip",sep='\t',header=None)
-    if qualityfile.endswith('.xz'):
-        ref_quality =  pd.read_csv(qualityfile,compression="xz",sep='\t',header=None)
-    if qualityfile.endswith('.bz2'):
-        ref_quality =  pd.read_csv(qualityfile,compression="bz2",sep='\t',header=None)
+def judgefunc(x):
+    x = x.split(',')
+    if 'high quality' in x:
+        return 'high quality'
     else:
-        ref_quality =  pd.read_csv(qualityfile,sep="\t",header=None)
+        return 'low quality'
+    
+
+def smorf_quality(outdir:str, qualityfile:str, resultfile:str) -> tuple:
+    result = pd.read_csv(resultfile,
+                         sep='\t',
+                         header=None)
+    
+    result = result.rename({0:'qseqid', 1:'sseqid'},
+                           axis=1)
+    
+    quality_file = path.join(outdir,
+                             "quality.out.smorfs.tsv")	
+    
+    ref_quality =  pd.read_table(qualityfile,
+                                 sep="\t",
+                                 header=None)
 
     ref_quality.columns = ['sseqid']
     ref_quality['quality'] = 'high quality'
 
-    output = pd.merge(result,ref_quality,how='left')[['qseqid', 'quality']].fillna('low quality')
-    output = output.groupby('qseqid',as_index=False,sort=False).agg({'quality':lambda x : ','.join(x.drop_duplicates())})
-    rule = {"high quality":0,"low quality":1}
-    output['quality'] = output['quality'].apply(lambda x: sorted(x.split(','),key=lambda x:rule[x])[0])
-
+    output = result.merge(on='sseqid',
+                          right=ref_quality,
+                          how='left')
+    
+    output = output[['qseqid', 'quality']]
+    output = output.fillna('low quality')
+    
+    output = output.groupby('qseqid',
+                            as_index=False,
+                            sort=False)
+    
+    output = output.agg({'quality': lambda x: ','.join(x.drop_duplicates())})  
+    output['quality'] = output['quality'].apply(lambda x: judgefunc(x))
+    
     number_dict = dict(output['quality'].value_counts())
     number_dict_normalize = dict(output['quality'].value_counts(normalize=True))
+    
     if "high quality" in number_dict.keys():
         number = number_dict['high quality']
         percentage = number_dict_normalize['high quality']
     else:
         number = 0
         percentage = 0
-    output.to_csv(quality_file,sep='\t',index=False)
-    return number,percentage
+    
+    output.to_csv(quality_file,
+                  sep='\t',
+                  index=False)
+    
+    return (number, percentage)