CASP16 system

Author: jl4mc

gate/feature/align_models_by_sequence.py

@@ -124,7 +124,7 @@ def get_sequence_by_chain(chain_id, sequence_id_map):
             return sequence
     return ''
 
-def merge_chain_pdbs(chain_mapping, outfile):
+def merge_chain_pdbs(chain_mapping, outfile, check_format=False):
     # reorder chains based on the stoichiometry
     # e.g., A2B2: AB CD
     chain_idx = 0
@@ -143,25 +143,36 @@ def merge_chain_pdbs(chain_mapping, outfile):
         #print(contents)
         fw.write('\n'.join(contents))
 
+    if check_format:
+        parser = PDBParser(QUIET=True)
+        try:
+            structure2 = parser.get_structure('', outfile)
+        except Exception as e:
+            os.system("rm " + outfile)
+
 
 def filter_single_model(inparams):
 
-    clustalw_program, sequence_id_map, inpdb, pdbdir, outpdb = inparams
+    clustalw_program, sequence_id_map, inpdb, pdbdir, outpdb, check_format = inparams
 
     # print(f"Filtering {inpdb}")
 
     # get chain mapping from pdb to fasta file
-    chain_mapping = get_chain_mapping(clustalw_program=clustalw_program,
-                                      sequence_id_map=sequence_id_map, 
-                                      inpdb=inpdb,
-                                      pdbdir=pdbdir)
-    # print(chain_mapping)
-    merge_chain_pdbs(chain_mapping, outpdb)
+    try:
+        chain_mapping = get_chain_mapping(clustalw_program=clustalw_program,
+                                        sequence_id_map=sequence_id_map, 
+                                        inpdb=inpdb,
+                                        pdbdir=pdbdir)
+        # print(chain_mapping)
+        merge_chain_pdbs(chain_mapping, outpdb, check_format)
+    except Exception as e:
+        print(f"Filtering {inpdb} failed!")
+        print(e)
 
     os.system(f"rm -rf {pdbdir}")
 
 
-def align_models(clustalw_program, fasta_path, outdir, input_model_dir):
+def align_models(clustalw_program, fasta_path, outdir, input_model_dir, check_format):
 
     # read sequences from fasta file
     sequences, descriptions = parse_fasta(open(fasta_path).read())
@@ -193,7 +204,7 @@ def align_models(clustalw_program, fasta_path, outdir, input_model_dir):
 
         makedir_if_not_exists(workdir)
 
-        process_list.append([clustalw_program, sequence_id_map, input_model_dir + '/' + model, workdir, outdir + '/' + model.replace('.pdb', '')])
+        process_list.append([clustalw_program, sequence_id_map, input_model_dir + '/' + model, workdir, outdir + '/' + model.replace('.pdb', ''), check_format])
 
     pool = Pool(processes=40)
     results = pool.map(filter_single_model, process_list)
@@ -208,8 +219,9 @@ def align_models(clustalw_program, fasta_path, outdir, input_model_dir):
     parser.add_argument('--outdir', type=str, required=True)
     parser.add_argument('--modeldir', type=str, required=True)
     parser.add_argument('--clustalw_program', type=str, required=True)
-    
+    parser.add_argument('--check_format', default=False, type=lambda x: (str(x).lower() == 'true'))
+
     args = parser.parse_args()
 
-    align_models(args.clustalw_program, args.fasta_path, args.outdir, args.modeldir)
+    align_models(args.clustalw_program, args.fasta_path, args.outdir, args.modeldir, args.check_format)
 

gate/feature/config.py

@@ -17,6 +17,7 @@
         'interface_pairwise_ts_script': os.path.join(ROOTDIR, 'gate', 'feature', 'interface_pairwise_ts.py'),
 
         'icps_script': os.path.join(ROOTDIR, 'gate', 'feature', 'generate_icps_scores.py'),
+        'model_size_script': os.path.join(ROOTDIR, 'gate', 'feature', 'generate_model_size.py'),
         'plddt_script': os.path.join(ROOTDIR, 'gate', 'feature', 'generate_plddt_scores.py'),
         'enqa_script': os.path.join(ROOTDIR, 'gate', 'feature', 'generate_enqa_scores.py'),
         'dproqa_script': os.path.join(ROOTDIR, 'gate', 'feature', 'generate_dproqa_scores.py'),

gate/feature/feature_generation.py

@@ -39,7 +39,7 @@ def generate_enqa_scores(enqa_env_path:str,
         model_size_ratio = dict(zip(list(model_info_df['model']), list(model_info_df['model_size_norm'])))
         target_dict['score_norm'] = []
 
-    max_length_threshold = 2500
+    max_length_threshold = 2600
     # read sequences from fasta file
     sequences, descriptions = parse_fasta(open(fasta_path).read())
     target_length = np.sum(np.array([len(sequence) for sequence in sequences]))

gate/feature/generate_enqa_scores.py

@@ -34,7 +34,7 @@ def generate_gcpnet_scores(gcpnet_ema_env_path: str,
                 f"data.ablate_esm_embeddings=false " \
                 f"model.ablate_af2_plddt=false " \
                 f"data.predict_output_dir={outdir}/workdir"
-    
+    print(cmd)    
     resultfile = f"{outdir}/result.csv"
 
     if not os.path.exists(resultfile):
@@ -52,7 +52,7 @@ def generate_gcpnet_scores(gcpnet_ema_env_path: str,
     model_size_ratio = {}
     if model_csv is not None and os.path.exists(model_csv):                
         model_info_df = pd.read_csv(model_csv)
-        model_size_ratio = dict(zip(list(model_info_df['model']), list(model_info_df['model_size_norm'])))
+        model_size_ratio = dict(zip([str(modelname) for modelname in list(model_info_df['model'])], list(model_info_df['model_size_norm'])))
         data_dict['score_norm'] = []
 
     pred_model_out_dir = os.path.join(outdir, 'pred_pdbs')
@@ -66,7 +66,7 @@ def generate_gcpnet_scores(gcpnet_ema_env_path: str,
         data_dict['score'] += [global_score / 100]
 
         if 'score_norm' in data_dict:
-            data_dict['score_norm'] += [global_score / 100 * float(model_size_ratio[modelname])]
+            data_dict['score_norm'] += [global_score / 100 * float(model_size_ratio[str(modelname)])]
 
         os.system(f"cp {pred_model} {pred_model_out_dir}/{modelname}")
 

gate/feature/generate_gcpnet_ema_scores.py

@@ -353,6 +353,10 @@ def generate_icps_scores(fasta_path, outdir, pairwise_score_csv, input_model_dir
     pairwise_df = pd.read_csv(pairwise_score_csv, index_col=[0])
     models = pairwise_df.columns
     tmscores = np.array([np.mean(np.array(pairwise_df[model])) for model in models])
+    # pairwise_df = pd.read_csv(pairwise_score_csv)
+    # models = pairwise_df['model']
+    # tmscores = pairwise_df['MMalign score']
+
     chain_pdbs = {}
     while True:
         select_model_idx = np.argmax(tmscores)
@@ -432,7 +436,7 @@ def generate_icps_scores(fasta_path, outdir, pairwise_score_csv, input_model_dir
 
             cal_list = [[cdpred_cmap_file, cmap_file] for cmap_file in cmap_files]
             # print(cal_list)
-            pool = Pool(processes=60)
+            pool = Pool(processes=150)
             results = pool.map(icps_recall_wrappeer, cal_list)
             pool.close()
             pool.join()