ENH update 04_forzen

Author: cocodyq

General_Scripts/00_Remove_redundancy_and_cluster/01_deduplicate_sort_merge_extract.py

@@ -87,10 +87,9 @@ def extract_seq(infile1,infile2,outfile1,outfile2):
     fastaset = set()
     with gzip.open(infile1,"rt") as f:
         for line in f:
-            line = line.strip()
-            linelist = line.split("\t")
-            if linelist[0] != "1":
-                fastaset.add(linelist[1])
+            count,seq = line.strip().split("\t")
+            if count != "1":
+                fastaset.add(seq)
 
     with gzip.open(outfile1, "wt", compresslevel=1) as out1, \
         gzip.open(outfile2, "wt", compresslevel=1) as out2:

General_Scripts/00_Remove_redundancy_and_cluster/07_identify_clusters.py

@@ -13,7 +13,7 @@ def identify(infile,outfile):
                     continue
                 else:
                     nameset.add(linelist[0])
-                    out.write(linelist[0]+"\t"+linelist[2]+"\n")
+                    out.write(f'{linelist[2]}\t{linelist[0]}\n')
 
 for i in range(24):
     INPUT_FILE_1 = "sub"+str(i)+".faa.gz.tsv"

General_Scripts/04_Frozen/01_rename_list.py

@@ -2,63 +2,65 @@
 Concept:
 Generate original name - 100AA rename list.
 Peptides are named: >GMSC10.100AA.XXX_XXX_XXX
-Numbers were assigned in order of increasing number of copies. 
-So that the lower the number, the lower the number of copies of that peptide was present in the input data. 
-And if the number of copies is same, numbers were assigned in order of letters of peptides.
+Numbers were assigned in order of increasing number of copies. If the number of copies is same, numbers were assigned in order of letters of peptides.
 '''
 
 def sort(infile,outfile,n,prefix):
     from operator import itemgetter   
     import gzip
-    seqnumber_list=[]
+
+    seqnumber_list = []
+
     with gzip.open(infile,"rt") as f1:
         for line in f1 :
-            line = line.strip()
-            linelist = line.split("\t")
-            if linelist[0] != "1":
-                seqnumber_tup = (int(linelist[0]),linelist[1])
+            count,seq = line.strip().split("\t")
+            if count != "1":
+                seqnumber_tup = (int(count),seq)
                 seqnumber_list.append(seqnumber_tup)
-    sortseqnumber_list=sorted(seqnumber_list,key=itemgetter(0,1))
+
+    sortseqnumber_list = sorted(seqnumber_list,key=itemgetter(0,1))
+
     with open(outfile,"wt") as out:
-        for i in range (len(seqnumber_list)):
+        for item in sortseqnumber_list:
             nf = f'{n:09}'
-            out.write(f'{sortseqnumber_list[i][0]}\t{sortseqnumber_list[i][1]}\t{prefix}.{nf[:3]}_{nf[3:6]}_{nf[6:9]}\n')
+            out.write(f'{item[0]}\t{item[1]}\t{prefix}.{nf[:3]}_{nf[3:6]}_{nf[6:9]}\n')
             n += 1
 
 def rename_nonsingleton(infile1,infile2,outfile):
     from fasta import fasta_iter
+
     fastadict={}
+
     for ID,seq in fasta_iter(infile1):
         fastadict[seq] = ID
-    out = open(outfile, "w")
-    with open (infile2) as f1:
-        for line in f1 :
-            line = line.strip()
-            linelist = line.split("\t")
-            name = fastadict[linelist[1]]
-            newname = linelist[2]
-            out.write(f'{name}\t{newname}\n')
-    out.close()   
+
+    with open(outfile,'wt') as out:
+        with open (infile2) as f1:
+            for line in f1 :
+                count,seq,newname = line.strip().split("\t")
+                name = fastadict[seq]
+                out.write(f'{name}\t{newname}\n')  
 
 def rename_singleton(infile1,infile2,outfile,n,prefix):
     from fasta import fasta_iter
+
     name=set()
-    with open (infile1) as f1:
+
+    with open(infile1) as f1:
         for line in f1 :
-            line = line.strip()
-            linelist = line.split("\t")
-            name.add(linelist[0])
-    out = open(outfile, "w")
-    for ID,seq in fasta_iter(infile2):
-        if ID in name:
-            nf = f'{n:09}'
-            out.write(f'{ID}\t{prefix}.{nf[:3]}_{nf[3:6]}_{nf[6:9]}\n')
-            n += 1
-    out.close()          
+            singleton,cluster = line.strip().split("\t")
+            name.add(singleton)
+
+    with open(outfile,'wt') as out:
+        for ID,seq in fasta_iter(infile2):
+            if ID in name:
+                nf = f'{n:09}'
+                out.write(f'{ID}\t{prefix}.{nf[:3]}_{nf[3:6]}_{nf[6:9]}\n')
+                n += 1     
 
 INPUT_FILE_1 = "metag_ProG.raw_number.tsv.gz"
 INPUT_FILE_2 = "metag_ProG_nonsingleton.faa.gz"
-INPUT_FILE_3 = "singleton_0.5_0.9.tsv"
+INPUT_FILE_3 = "singleton_0.9.tsv"
 INPUT_FILE_4 = "metag_ProG_singleton.faa.gz"
 OUTPUT_FILE_1 = "nonsingleton_rename_seq.tsv"
 OUTPUT_FILE_2 = "nonsingleton_rename.tsv"

General_Scripts/04_Frozen/02_100AA_faa_fna.py

@@ -6,27 +6,27 @@
 def getseq(infile1,infile2,outfile):   
     from fasta import fasta_iter
     import lzma
+    
     name = {}
-    out1 = lzma.open(outfile, "wt")
+    out = lzma.open(outfile, "wt")
 
-    with lzma.open(infile1,"rt") as f1:
+    with open(infile1,"rt") as f1:
         for line in f1:
-            line = line.strip()
-            linelist = line.split("\t")
-            name[linelist[0]] = linelist[1] 
+            old,new = line.strip().split("\t")
+            name[old] = new
 
     for ID,seq in fasta_iter(infile2):
         if ID in name.keys():
-            out1.write(f'>{name[ID]}\n{seq}\n')           
-    out1.close()        
+            out.write(f'>{name[ID]}\n{seq}\n')           
+    out.close()        
 
-INPUT_FILE_1 = "./data/100AA_rename.tsv.xz"  
-INPUT_FILE_2 = "./data/metag_ProG_dedup.faa.gz"
+INPUT_FILE_1 = "100AA_rename.tsv"  
+INPUT_FILE_2 = "metag_ProG_dedup.faa.gz"
 INPUT_FILE_3 = "GMSC10.metag_smorfs.fna.xz"
 INPUT_FILE_4 = "GMSC.ProGenomes2.smorfs.fna.xz"
-OUTPUT_FILE_1 = "./data/frozen/100AA_GMSC.faa.xz"
-OUTPUT_FILE_2 = "./data/frozen/100AA_metag.fna.xz"
-OUTPUT_FILE_3 = "./data/frozen/100AA_prog.fna.xz"
+OUTPUT_FILE_1 = "100AA_GMSC.faa.xz"
+OUTPUT_FILE_2 = "100AA_metag.fna.xz"
+OUTPUT_FILE_3 = "100AA_prog.fna.xz"
 
 getseq(INPUT_FILE_1,INPUT_FILE_2,OUTPUT_FILE_1)
 getseq(INPUT_FILE_1,INPUT_FILE_3,OUTPUT_FILE_2)

General_Scripts/04_Frozen/03_90AA_faa_fna.py

@@ -11,105 +11,98 @@
 '''
 Generate original name - 90AA rename list.
 Peptides are named: >GMSC10.90AA.XXX_XXX_XXX
-Numbers were assigned in order of increasing number of copies. 
-So that the lower the number, the lower the number of copies of that peptide was present in the input data. 
-And if the number of copies is same, numbers were assigned in order of letters of peptides.
 '''
 def rename(infile1,infile2,outfile,n,prefix):   
     number = {}
-    seqnumber_list=[]
+    seqnumber_list = []
 
     with gzip.open(infile2,"rt") as f2:
         for line in f2 :
-            line = line.strip()
-            linelist = line.split("\t")
-            if linelist[0] in number.keys():
-                number[linelist[0]] += 1
+            cluster,member = line.strip().split("\t")
+            if cluster in number.keys():
+                number[cluster] += 1
             else:
-                number[linelist[0]] = 1
+                number[cluster] = 1
 
     for ID,seq in fasta_iter(infile1):
         seqnumber_tup = (int(number[ID]),seq,ID) 
         seqnumber_list.append(seqnumber_tup)
 
-    sortseqnumber_list=sorted(seqnumber_list,key=itemgetter(0,1))
+    sortseqnumber_list = sorted(seqnumber_list,key=itemgetter(0,1))
     with lzma.open(outfile,"wt") as out:
-        for i in range (len(sortseqnumber_list)):
+        for item in sortseqnumber_list:
             nf = f'{n:09}'
-            out.write(f'{sortseqnumber_list[i][2]}\t{prefix}.{nf[:3]}_{nf[3:6]}_{nf[6:9]}\n')
-            n += 1
+            out.write(f'{item[2]}\t{prefix}.{nf[:3]}_{nf[3:6]}_{nf[6:9]}\n')
+            n += 1 
 
 '''
-Generate originalname - 100AA - 90AA rename list.
+Generate original name - 100AA - 90AA rename list.
 '''
 def rename_all(infile1,infile2,outfile):   
     name = {}
     out1 = lzma.open(outfile, "wt")
 
     with lzma.open(infile1,"rt") as f1:
         for line in f1:
-            line = line.strip()
-            linelist = line.split("\t")
-            name[linelist[0]] = linelist[1] 
+            old,new = line.strip().split("\t")
+            name[old] = new
 
     with gzip.open(infile2,"rt") as f2:
         for line in f2:
-            line = line.strip().strip(">")
-            linelist = line.split("\t")
-            if linelist[0] in name.keys():
-                out1.write(linelist[0]+"\t"+linelist[1]+"\t"+name[linelist[0]]+"\n")         
-    out1.close()   
+            old,new = line.strip().split("\t")
+            if old in name.keys():
+                out1.write(f'{old}\t{new}\t{name[new]}\n')         
+    out1.close()
 
 '''
 Generate rename and sequence of 90AA faa.
 '''
 def getfaa(infile1,infile2,outfile):   
     name = {}
-    out1 = lzma.open(outfile, "wt")
+    out = lzma.open(outfile, "wt")
 
     with lzma.open(infile1,"rt") as f1:
         for line in f1:
-            line = line.strip()
-            linelist = line.split("\t")
-            name[linelist[0]] = linelist[1] 
-            
-    
+            old,new = line.strip().split("\t")
+            name[old] = new
+                
     for ID,seq in fasta_iter(infile2):
-        out1.write(f'>{name[ID]}\n{seq}\n')            
-    out1.close()  
+        out.write(f'>{name[ID]}\n{seq}\n')            
+    out.close()  
 
 '''
 Generate rename and sequence of 90AA fna.
 '''
 def getfna(infile1,infile2,outfile):       
     fasta = {}
     table = {}
-    out1 = lzma.open(outfile, "wt")
+
+    out = lzma.open(outfile, "wt")
+
     with lzma.open(infile1,"rt") as f1:
         for line in f1:
-            line = line.strip()
-            linelist = line.split("\t")
-            table[linelist[1]] = linelist[2]
+            old,name100,name90 = line.strip().split("\t")
+            table[name100] = name90
 
     for ID,seq in fasta_iter(infile2):
         if ID in table.keys():           
             fasta[table[ID]] = seq
-    table = {}
+
     for ID,seq in sorted(fasta.items()):
-        out1.write(f">{ID}\n{seq}\n")               
-    out1.close()  
+        out.write(f">{ID}\n{seq}\n")               
+    out.close()  
 
-INPUT_FILE_1 = "./clust_result/0.5_result/metag_ProG_nonsingleton_0.9_clu_rep.faa.gz"
-INPUT_FILE_2 = "./clust_result/0.5_result/metag_ProG_nonsingleton_0.9_clu.tsv.gz"
-INPUT_FILE_3 = "./data/100AA_rename.tsv.xz"
-INPUT_FILE_4 = "./data/frozen/100AA_GMSC.fna.xz"
+INPUT_FILE_1 = "metag_ProG_nonsingleton_0.9_clu_rep.faa.gz"
+INPUT_FILE_2 = "metag_ProG_nonsingleton_0.9_clu.tsv.gz"
+INPUT_FILE_3 = "100AA_rename.tsv.xz"
+INPUT_FILE_4 = "100AA_GMSC.fna.xz"
 
-OUTPUT_FILE_1 = "./data/frozen/90AA_rename.tsv.xz"
-OUTPUT_FILE_2 = "./data/frozen/90AA_rename_all.tsv.xz"
-OUTPUT_FILE_3 = "./data/frozen/90AA_GMSC.faa.xz"
-OUTPUT_FILE_4 = "./data/frozen/90AA_GMSC.fna.xz"
+OUTPUT_FILE_1 = "90AA_rename.tsv.xz"
+OUTPUT_FILE_2 = "90AA_rename_all.tsv.xz"
+OUTPUT_FILE_3 = "90AA_GMSC.faa.xz"
+OUTPUT_FILE_4 = "90AA_GMSC.fna.xz"
 
 rename(INPUT_FILE_1,INPUT_FILE_2,OUTPUT_FILE_1,0,'GMSC10.90AA')
 rename_all(OUTPUT_FILE_1,INPUT_FILE_3,OUTPUT_FILE_2)
 getfaa(OUTPUT_FILE_1,INPUT_FILE_1,OUTPUT_FILE_3)
-getfna(OUTPUT_FILE_2,INPUT_FILE_4,OUTPUT_FILE_4 )
+getfna(OUTPUT_FILE_2,INPUT_FILE_4,OUTPUT_FILE_4)

General_Scripts/04_Frozen/04_family.py

@@ -1,43 +1,36 @@
 '''
 Concept:
-Generate the table including the name of smORFs, and the clusters they belong to at 90% identity.
+Generate the table including the 100AA smORFs and 90AA families.
 '''
 
 def generate_family(infile1,infile2,infile3,outfile):   
     import lzma
     import gzip
-    name50 = {}
+
     name90 = {}
-    out1 = lzma.open(outfile, "wt")
-    
-    with lzma.open(infile1,"rt") as f1:
-        for line in f1:
-            line = line.strip()
-            linelist = line.split("\t")
-            name50[linelist[1]] = linelist[2] 
-            
-    with lzma.open(infile2,"rt") as f2:
-        for line in f2:
-            line = line.strip()
-            linelist = line.split("\t")
-            name90[linelist[1]] = linelist[2]   
+    out = open(outfile, "wt")
 
-    with gzip.open(infile3,"rt") as f3:
-        for line in f3:
-            line = line.strip()
-            linelist = line.split("\t")
-            if len(linelist) == 2:
-                out1.write(linelist[0]+"\t"+""+"\t"+name50[linelist[1]]+"\n")
-            elif len(linelist) == 3 and linelist[1] != "":
-                out1.write(linelist[0]+"\t"+name90[linelist[2]]+"\t"+name50[linelist[1]]+"\n")
-            else:
-                out1.write(linelist[0]+"\t"+name90[linelist[2]]+"\n")           
-    out1.close()        
+    with lzma.open(infile1,"rt") as f:
+        for line in f:
+            old,new100,new90 = line.strip().split("\t")
+            name90[new100] = new90
+    
+    name100 = {}
+    with open(infile2,'rt') as f:
+        for line in f:
+            old,new = line.strip().split("\t")
+            name100[old] = new
+
+    with gzip.open(infile3,"rt") as f:
+        for line in f:
+            cluster,member = line.strip().split("\t")
+            out.write(f'{name100[member]}\t{name90[cluster]}\n')
+    out.close()        
 
-INPUT_FILE_1 = "./data/frozen/50AA_rename_all.tsv.xz"  
-INPUT_FILE_2 = "./data/frozen/90AA_rename_all.tsv.xz"
-INPUT_FILE_3 = "./clust_result/result/all_0.5_0.9_rename.tsv.gz"
-OUTPUT_FILE = "./data/frozen/all_0.9_0.5_family.tsv.xz"
+INPUT_FILE_1 = "90AA_rename_all.tsv.xz"
+INPUT_FILE_2 = "100AA_rename.tsv"  
+INPUT_FILE_3 = "all_0.9.tsv.gz"
+OUTPUT_FILE = "GMSC.cluster.tsv"
 
 generate_family(INPUT_FILE_1,INPUT_FILE_2,INPUT_FILE_3,OUTPUT_FILE)
 

General_Scripts/04_Frozen/README.md

@@ -1,7 +1,8 @@
-# 04_Frozen
+## 04_Frozen
+
 | **Code** | **Description** | **Input** | **Output** |
 | :---: | :---: | :---: | :---: |
-| 01_rename_list.py | Rename 100AA sequences | metag_ProG.raw_number.tsv.gz metag_ProG_nonsingleton.faa.gz singleton_0.5_0.9.tsv metag_ProG_singleton.faa.gz| nonsingleton_rename.tsv singleton_rename.tsv |
-| 02_100AA_faa_fna.py | Generate 100AA faa and fna file with new identifier | 100AA_rename.tsv.xz metag_ProG_dedup.faa.gz GMSC10.metag_smorfs.fna.xz GMSC.ProGenomes2.smorfs.fna.xz | 100AA_GMSC.faa.xz 100AA_metag.fna.xz 100AA_prog.fna.xz |
+| 01_rename_list.py | Rename 100AA sequences | metag_ProG.raw_number.tsv.gz metag_ProG_nonsingleton.faa.gz singleton_0.9.tsv metag_ProG_singleton.faa.gz| nonsingleton_rename.tsv singleton_rename.tsv |
+| 02_100AA_faa_fna.py | Generate 100AA faa and fna file with new identifier | 100AA_rename.tsv metag_ProG_dedup.faa.gz GMSC10.metag_smorfs.fna.xz GMSC.ProGenomes2.smorfs.fna.xz | 100AA_GMSC.faa.xz 100AA_metag.fna.xz 100AA_prog.fna.xz |
 | 03_90AA_faa_fna.py | Rename 90AA sequences and generate 90AA faa and fna file with new identifier | metag_ProG_nonsingleton_0.9_clu_rep.faa.gz metag_ProG_nonsingleton_0.9_clu.tsv.gz 100AA_rename.tsv.xz 100AA_GMSC.fna.xz | 90AA_rename.tsv.xz 90AA_rename_all.tsv.xz 90AA_GMSC.faa.xz 90AA_GMSC.fna.xz |
-| 04_family.py | Generate the cluster table | 90AA_rename_all.tsv.xz all_0.5_0.9_rename.tsv.gz | all_0.9_0.5_family.tsv.xz |
+| 04_family.py | Generate the family table | 90AA_rename_all.tsv.xz 100AA_rename.tsv all_0.9.tsv.gz | GMSC.cluster.tsv |