ENH update 08_Conserved_domain_annotation

Author: cocodyq

General_Scripts/02_Habitat_mapping/02_multi_habitat.py

@@ -86,7 +86,7 @@ def general(infile1,infile2,outfile):
 INPUT_FILE_3 = "habitat_general.txt"
 OUTPUT_FILE_1 = "all_cluster_multi_habitat.tsv.xz"
 OUTPUT_FILE_2 = "100AA_multi_habitat.tsv.xz" 
-OUTPUT_FILE_3 = "100AA_multi_general_habitat.tsv.xz"
+OUTPUT_FILE_3 = "GMSC10.100AA.general_habitat.tsv.xz"
 
 multi_habitat(INPUT_FILE_1,OUTPUT_FILE_1)
 extract(INPUT_FILE_2,OUTPUT_FILE_1,OUTPUT_FILE_2)

General_Scripts/02_Habitat_mapping/03_map_cluster_habitat.py

@@ -84,12 +84,12 @@ def general(infile1,infile2,outfile):
 
     out.close()
 
-INPUT_FILE_1 = "100AA_multi_general_habitat.tsv.xz"   
+INPUT_FILE_1 = "GMSC10.100AA.general_habitat.tsv.xz"   
 INPUT_FILE_2 = "GMSC.cluster.tsv.gz"
 INPUT_FILE_3 = "habitat_general.txt"
 OUTPUT_FILE_1 = "cluster_multi_habitat_90.tsv.xz" 
 OUTPUT_FILE_2 = "90AA_multi_habitat.tsv.xz"
-OUTPUT_FILE_3 = "90AA_multi_general_habitat.tsv.xz"
+OUTPUT_FILE_3 = "GMSC10.90AA.general_habitat.tsv.xz"
 
 mapcluster(INPUT_FILE_1,INPUT_FILE_2,OUTPUT_FILE_1)
 multi_habitat(OUTPUT_FILE_1,OUTPUT_FILE_2)

General_Scripts/02_Habitat_mapping/README.md

@@ -3,5 +3,5 @@
 | **Code** | **Description** | **Input** | **Output** |
 | :---: | :---: | :---: | :---: |
 | 01_map_habitat.py | Map habitat for all the smORFs from metaG | metadata.tsv GMSC10.metag_smorfs.rename.txt.xz dedup_cluster.tsv.gz| metag_cluster_habitat.tsv.xz |
-| 02_multi_habitat.py | Combine multiple habitats for each smORF from the same cluster | metag_cluster_habitat.tsv.xz GMSC.cluster.tsv.gz habitat_general.txt| 100AA_multi_general_habitat.tsv.xz |
-| 03_map_cluster_habitat.py | Map multiple habitats to 90% identity smORFs clusters. | GMSC.cluster.tsv.gz 100AA_multi_general_habitat.tsv.xz habitat_general.txt| 90AA_multi_general_habitat.tsv.xz |
+| 02_multi_habitat.py | Combine multiple habitats for each smORF from the same cluster | metag_cluster_habitat.tsv.xz GMSC.cluster.tsv.gz habitat_general.txt| GMSC10.100AA.general_habitat.tsv.xz |
+| 03_map_cluster_habitat.py | Map multiple habitats to 90% identity smORFs clusters. | GMSC.cluster.tsv.gz 100AA_multi_general_habitat.tsv.xz habitat_general.txt| GMSC10.90AA.general_habitat.tsv.xz |

General_Scripts/08_Conserved_domain_annotation/01_Annotation/01_cdd.sh

@@ -3,4 +3,4 @@
 # Concept: 
 # Map to CDD database
 
-rpsblast -query 90AA_GMSC.faa -out 90AA_cdd.tsv -db ~/Cdd -num_threads 20 -evalue 0.01 -outfmt "6 qseqid sseqid qlen score length pident evalue"
+rpsblast -query 90AA_GMSC.faa -out 90AA_cdd.tsv -db ./Cdd -num_threads 20 -evalue 0.01 -outfmt "6 qseqid sseqid qlen score length pident evalue"

General_Scripts/08_Conserved_domain_annotation/01_Annotation/02_add_pssm_length.py

@@ -1,4 +1,5 @@
 '''
+Concept:
 Add length of PSSM and filter with target coverage >80%.
 '''
 def add_length(infile1,infile2,outfile):
@@ -23,6 +24,7 @@ def add_length(infile1,infile2,outfile):
                         out.write(f'{line}\t{cdd_dict[pssm]}\n')
 
 def filter_cov(infile,outfile):
+    import pandas as pd
     result = pd.read_csv(infile,compression='gzip',sep='\t',header=None,names=['smorf','cdd','query_length','score','align_length','identity','evalue','target_length'])
     result['tcov'] = result['align_length']/result['target_length']
     result = result[result['tcov'] >0.8]

General_Scripts/08_Conserved_domain_annotation/02_multi-phylum_analysis/03_multi_specific.py renamed to General_Scripts/08_Conserved_domain_annotation/02_multi-phylum_analysis/01_multi_specific.py

@@ -1,5 +1,6 @@
 '''
-Analyse if smORFs 90AA families are specific or multiple at the taxonomy rank.
+Concept:
+Analyse if 90AA families are specific or multiple at the taxonomy rank.
 '''
 
 import pandas as pd
@@ -51,26 +52,6 @@ def cal(infile,outfile):
         out.write(f'{smorf}\t{flag_multi}\t{flag_specific}\t{number}\n')
     out.close()
 
-# change to new identifier
-
-def store(infile1):
-    import lzma
-    name = {}
-    with lzma.open(infile1,'rt') as f1:
-        for line in f1:
-            old,new = line.strip().split('\t')
-            name[old] = new
-    return name
-
-def map_multi(name,infile,outfile):
-    out1 = open(outfile,'wt')
-    with open(infile,'rt') as f2:
-        for line in f2:
-            family,anno = line.strip().split('\t',1)
-            if family in name.keys():
-                out1.write(f'{name[family]}\t{anno}\n')
-    out1.close()
-
 # Calculate number of taxonomy specific
 def merge(infile,outfile):
     km = 0
@@ -158,12 +139,8 @@ def merge(infile,outfile):
         out.write(f'kingdom-specific\t{ok}\t{pm}\t{ps}\n')
 
 infile1 = 'metag_cluster_tax_90.tsv'
-infile2 = '90AA_rename.tsv.xz'
 outfile1 = '90AA_taxa_multi_specific.tsv'
-outfile2 = '90AA_multi_newname.tsv'
-outfile3 = '90AA_specific_multi.tsv'
+outfile2 = '90AA_specific_multi.tsv'
 
 cal(infile1,outfile1)
-name = store(infile2)
-map_multi(name,outfile1,outfile2)
-merge(outfile1,outfile3)
+merge(outfile1,outfile2)

General_Scripts/08_Conserved_domain_annotation/02_multi-phylum_analysis/01_get_all_habitat.py renamed to General_Scripts/08_Conserved_domain_annotation/02_multi-phylum_analysis/02_get_all_habitat.py

@@ -1,4 +1,5 @@
 '''
+Concept:
 Select 90AA smORF families across all 8 habitat categories.
 Add cdd annotation to these 90AA smORF families.
 '''
@@ -113,7 +114,7 @@ def merge_habitat_motif(infile1,infile2,infile3,outfile):
 infile2 = '1_cdd_tcov_90AA.tsv.gz'
 infile3 = 'cddid_all.tbl.gz'
 outfile1 = 'all_habitat_smorf.tsv'
-outfile2 = 'all_habitat_motif_right.tsv'
+outfile2 = 'all_habitat_smorf_motif.tsv'
 
 map_high(infile1,outfile1)
 merge_habitat_motif(infile2,outfile1,infile3,outfile2)

General_Scripts/08_Conserved_domain_annotation/02_multi-phylum_analysis/02_species_number.py renamed to General_Scripts/08_Conserved_domain_annotation/02_multi-phylum_analysis/03_species_number.py

@@ -1,4 +1,5 @@
 '''
+Concept:
 Calculate the species number in 90AA families.
 '''
 def tax_format(infile,outfile):
@@ -27,42 +28,26 @@ def cal_species(infile,outfile):
         for key,value in species_number.items():
             out.write(f'{key}\t{len(value)}\n')
 
-def store(infile):
-    import lzma
-    name = {}
-    with lzma.open(infile,'rt') as f:
-        for line in f:
-            old,new = line.strip().split('\t')
-            name[old] = new
-    return name
-
-def select(infile2):
+def count(infile1,infile2,outfile):
     housekeeping = set()
-    with open(infile2,'rt') as f2:
-        for line in f2:
+    with open(infile1,'rt') as f:
+        for line in f:
             smorf,habitat = line.strip().split('\t')
             housekeeping.add(smorf)
-    return housekeeping
 
-def count(name,housekeeping,infile,outfile):
-    out = open(outfile,'wt')
-    with open(infile,'rt') as f1:
-        for line in f1:
-            smorf,number = line.strip().split('\t')
-            if smorf in name.keys():
-                if name[smorf] in housekeeping:
-                    out.write(f'{name[smorf]}\t{number}\n')
-    out.close()
+    with open(outfile,'wt') as out:
+        with open(infile2,'rt') as f:
+            for line in f:
+                smorf,number = line.strip().split('\t')
+                if smorf in housekeeping:
+                    out.write(f'{smorf}\t{number}\n')
 
 infile1 = "metag_cluster_tax_90.tsv.xz"
-infile2 = '90AA_rename.tsv.xz'
-infile3 = 'all_habitat_smorf.tsv'
+infile2 = 'all_habitat_smorf.tsv'
 outfile1 = "metag_cluster_tax_90.tsv"
 outfile2 = '90AA_species_number.tsv'
 outfile3 = 'housekeeping_species.tsv'
 
 tax_format(infile1,outfile1)
 cal_species(outfile1,outfile2)
-name = store(infile2)
-housekeeping = select(infile3)
-count(name,housekeeping,outfile2,outfile3)
+count(infile2,outfile2,outfile3)

General_Scripts/08_Conserved_domain_annotation/02_multi-phylum_analysis/04_merge_multi-phylum.py

@@ -31,7 +31,7 @@ def merge():
     df.to_csv('housekeeping_motif_species_multi_phylum_all.tsv',sep='\t',index=None)
 
 infile1 = 'housekeeping_species.tsv'
-infile2 = '90AA_multi_newname.tsv'
+infile2 = '90AA_taxa_multi_specific.tsv'
 outfile = 'housekeeping_multi.tsv'
 
 map_multi(infile1,infile2,outfile)

General_Scripts/08_Conserved_domain_annotation/02_multi-phylum_analysis/05_map_sample_taxonomy.py

@@ -19,12 +19,11 @@ def map_90(sample_100,infile2,outfile):
     sample_90 = {}
     with lzma.open(infile2,'rt') as f2:
         for line in f2:
-            linelist = line.strip().split('\t')
-            if linelist[1] != '':
-                if linelist[0] in sample_100.keys():
-                    if linelist[1] not in sample_90.keys():
-                        sample_90[linelist[1]] = []
-                    sample_90[linelist[1]].append(sample_100[linelist[0]])
+            member,cluster = line.strip().split('\t')
+            if member in sample_100.keys():
+                if cluster not in sample_90.keys():
+                    sample_90[cluster] = []
+                sample_90[cluster].append(sample_100[member])
     sample_100 = {}
     with open(outfile,'wt') as out:
         for key,value in sample_90.items():
@@ -35,47 +34,44 @@ def map_sample(infile1,infile2,outfile):
     new = set()
     with open(infile1,'rt') as f1:
         for line in f1:
-            new.add(line.strip())
+            smorf,number = line.strip().split('\t')
+            if number > 100:
+                new.add(smorf)
 
     with open(outfile,'wt') as out:
         with open(infile2,'rt') as f2:
             for line in f2:
-                linelist = line.strip().split('\t')
-                if linelist[0] in new:
+                cluster,sample = line.strip().split('\t')
+                if cluster in new:
                     out.write(line)
 
-def map_taxonomy(infile1,infile2,infile3,outfile):
+def map_taxonomy(infile1,infile2,outfile):
     import lzma
-    new = set()
-    old = {}
-    with open(infile1,'rt') as f1:
-        for line in f1:
-            new.add(line.strip())
+    seqs = set()
 
-    with lzma.open(infile2,'rt') as f2:
-        for line in f2:
-            linelist = line.strip().split('\t')
-            if linelist[1] in new:
-                old[linelist[0]] = linelist[1]
+    with open(infile1,'rt') as f:
+        for line in f:
+            smorf,number = line.strip().split('\t')
+            if number > 100:
+                seqs.add(smorf)
 
     with open(outfile,'wt') as out:
-        with lzma.open(infile3,'rt') as f3:
-            for line in f3:
-                line = line.strip()
-                linelist = line.split('\t',2)
-                if linelist[0] in old:
-                    out.write(f'{old[linelist[0]]}\t{line}\n')
+        with lzma.open(infile2,'rt') as f:
+            for line in f:
+                cluster,member,taxonomy = line.strip().split('\t',2)
+                if cluster in seqs:
+                    out.write(line)
 
 infile1 = '100AA_sample.tsv.xz'
-infile2 = 'all_0.9_0.5_family_sort.tsv.xz'
-infile3 = 'housekeeping.txt'
-infile4 = '90AA_rename.tsv.xz'
-infile5 = 'metag_cluster_tax_90.tsv.xz'
+infile2 = 'GMSC.cluster.tsv.gz'
+infile3 = 'housekeeping_species.tsv'
+infile4 = 'metag_cluster_tax_90.tsv.xz'
+
 outfile1 = '90AA_sample.tsv'
 outfile2 = 'housekeeping_sample.txt'
 outfile3 = 'housekeeping_taxonomy.txt'
 
 sample_100 = store_100(infile1)
 map_90(sample_100,infile2,outfile1)
 map_sample(infile3,outfile1,outfile2)
-map_taxonomy(infile3,infile4,infile5,outfile3)
+map_taxonomy(infile3,infile4,outfile3)