#!/usr/bin/env python # -*- encoding: utf-8 -*- ''' @file :analysis_pdb.py @Description: : 分析PDB结构 @Date :2023/11/27 09:57:00 @Author :lyzeng @Email :pylyzeng@gmail.com @version :1.0 ''' # micromamba create -n modeller modeller biopython pymol-open-source biopandas requests -y -c conda-forge -c salilab # modeller注册码：MODELIRANJE from dataclasses import dataclass, field from Bio.PDB import PDBParser from Bio.SeqUtils import seq1 from typing import List, Dict, Tuple, Optional from functools import reduce, partial from Bio.PDB import MMCIFIO, PDBIO, Chain, Structure from biopandas.pdb import PandasPdb from pathlib import Path from Bio.SeqRecord import SeqRecord from Bio import SeqIO import requests from copy import deepcopy from pymol import cmd @dataclass class PDBAnalyzer: """ PDBAnalyzer class for analyzing protein structures from PDB files. Attributes: pdb_file (Path): Path to the PDB file to be analyzed. structure (object): The parsed PDB structure, initialized in __post_init__. """ pdb_file: Path pid: Optional[str] = field(default=None, init=False) structure: object = field(init=False) biodf: PandasPdb = field(init=False) protein_state: str = field(init=False) # Apo or Holo chain_id_list: List[str] = field(init=False) def __post_init__(self): """ Initialize the PDB structure after the object is created. """ parser = PDBParser(QUIET=True) self.pid = self.pdb_file.stem.lower() if len(self.pdb_file.stem) == 4 else None self.structure = parser.get_structure('PDB_structure', self.pdb_file.as_posix()) self.biodata = PandasPdb().read_pdb(self.pdb_file.as_posix()) self.protein_state = 'Holo' if 'HETATM' in self.biodata.df.keys() else 'Apo' self.chain_id_list = self.biodata.df['ATOM']['chain_id'].drop_duplicates().to_list() def check_continuity(self, chain, missing_char): """ Check the continuity of residues in a protein chain. Args: chain (Chain): A Chain object from Bio.PDB. missing_char (str): Character to denote missing residues. Returns: List[Tuple[str, int]]: List of tuples containing residue names and their sequence numbers. """ def reducer(acc, residue): # Accumulate residues, filling in gaps with the specified missing character expected_resseq = acc[-1][1] + 1 if acc else residue.id[1] while expected_resseq < residue.id[1]: acc.append((missing_char, expected_resseq)) expected_resseq += 1 acc.append((residue.resname, residue.id[1])) return acc return reduce(reducer, chain, []) def residue_to_single_letter(self, residue_name): """ Convert a three-letter residue name to a single-letter code. Args: residue_name (str): Three-letter residue name. Returns: str: Single-letter residue code or the original string if it's a special character. """ # Handle special characters separately if residue_name in {'-', 'X', ''}: return residue_name # Convert standard residue names to single-letter codes return seq1(residue_name) def extract_sequences(self, missing_char='-') -> Dict[str, str]: """ Extract amino acid sequences from the structure, with gaps denoted by the specified character. Args: missing_char (str): Character to use for missing residues (default is '-'). Returns: Dict[str, str]: Dictionary of chain IDs mapped to their amino acid sequences. """ # Create a continuity check function with the specified missing character check_continuity_new = partial(self.check_continuity, missing_char=missing_char) sequences = {} # Process each chain in the structure for model in self.structure: chains = model.get_list() for chain in chains: # Check continuity and get the sequence of residues chain_sequence = check_continuity_new(chain) # Convert to single-letter sequence single_letter_sequence = ''.join(self.residue_to_single_letter(res[0]) for res in chain_sequence) sequences[chain.id] = single_letter_sequence return sequences def extract_sequences_info(self) -> Dict[str, List[int]]: """ Get missing residues information for each chain in the structure. Utilizes map and filter for efficient processing. Returns: Dict[str, List[int]]: Dictionary mapping each chain ID to a list of missing residue numbers. """ def find_missing(chain): observed = [residue.get_id()[1] for residue in chain if residue.id[0] == ' '] if observed: full_range = range(min(observed), max(observed) + 1) return chain.get_id(), sorted(set(full_range) - set(observed)) return chain.get_id(), [] chains = [chain for model in self.structure for chain in model] missing_info = map(find_missing, chains) return dict(filter(lambda x: x[1], missing_info)) def extract_chain(self, chain_id: str): """ Extract a specific chain from the structure. use biopython Args: chain_id (str): ID of the chain to be extracted. Returns: object: An object containing the extracted chain with a save method. """ target_chain = next((chain for model in self.structure for chain in model if chain.id == chain_id), None) if target_chain is None: raise ValueError(f"Chain {chain_id} not found in the structure.") def save_to_file(filename: str, file_format: str = 'pdb'): """ Save the extracted chain to a file in the specified format. Args: filename (str): Name of the file to save. file_format (str): Format of the file ('pdb' or 'cif'). Default is 'pdb'. """ if file_format not in ['pdb', 'cif']: raise ValueError("Unsupported file format. Use 'pdb' or 'cif'.") if file_format == 'pdb': io = PDBIO() elif file_format == 'cif': io = MMCIFIO() io.set_structure(target_chain) io.save(filename) return type('ChainExtractor', (object,), {'save': save_to_file}) def get_missing_residues_as_loops(self) -> Dict[str, List[Tuple[int, int, int]]]: """ Get missing residues information formatted as loops for PyRosetta, along with chain identifiers. Returns: Dict[str, List[Tuple[int, int, int]]]: Dictionary mapping chain IDs to lists of tuples representing loop regions (start, end, cut). """ missing_residues = self.extract_sequences_info() loops_by_chain = {} for chain_id, missing_res_nums in missing_residues.items(): loops = [] for num in missing_res_nums: # Define the loop region, assuming simple start and end points loops.append((num - 1, num + 1, num)) # Example loop definition loops_by_chain[chain_id] = loops return loops_by_chain def change_chain_identifier(self, old_chain_id: str, new_chain_id: str, split: bool = False) -> PandasPdb: """ Change the identifier of a specific chain. default use biopandas , not test yet. Args: old_chain_id (str): Original identifier of the chain. new_chain_id (str): New identifier to be assigned to the chain. return: object: An new object,which not influce original self.biodata, containing the changed chain with a save method. """ biodata = self.split_chain(old_chain_id) if split else deepcopy(self.biodata) # get ATOM colum biodata.df['ATOM']['chain_id'] = biodata.df['ATOM']['chain_id'].replace(old_chain_id.strip().upper(), new_chain_id.strip().upper()) return biodata # use to_pdb to save def change_chain_identifier_biopython(self, old_chain_id: str, new_chain_id: str, split: bool = False) -> PDBIO: io = PDBIO() def change_id(chain): if chain.id == old_chain_id: chain.id = new_chain_id return chain if split: # 使用 filter 和 map 创建一个只包含已更改的特定链的新结构 new_structure = Structure.Structure("modified_structure") for model in self.structure: filtered_chains = filter(lambda c: c.id == old_chain_id, model) changed_chains = map(change_id, filtered_chains) new_model = reduce(lambda m, c: m.add(c) or m, changed_chains, Structure.Model.Model(model.id)) new_structure.add(new_model) io.set_structure(new_structure) else: # 直接在整个结构上更改ID all_chains = (chain for model in self.structure for chain in model) list(map(change_id, all_chains)) io.set_structure(self.structure) return io def split_chain(self, chain_id: str) -> PandasPdb: biodata = deepcopy(self.biodata) df = biodata.df['ATOM'] # get ATOM colum tmp = PandasPdb() # create empty pdb file object tmp.df['ATOM'] = df[df['chain_id'] == chain_id.strip().upper()] return tmp # use to_pdb method save def save_biodata(self, out_file: Path): self.biodata.to_pdb(out_file.as_posix()) @staticmethod def write_seq_to_fasta_single_line(seq_record: SeqRecord, output_file: Path): with open(output_file, 'w') as output_handle: output_handle.write(f">{seq_record.description}\n") output_handle.write(str(seq_record.seq) + "\n") def download_fasta(self, pdb_id: str = None): if pdb_id == None: pdb_id = self.pid if pdb_id == None: raise ValueError("PDB ID not found.") url = f"https://www.rcsb.org/fasta/entry/{pdb_id.upper()}" response = requests.get(url) if response.status_code == 200: output_file = self.pdb_file.parent / f"{pdb_id}.fasta" with open(output_file, 'w') as file: file.write(response.text) return output_file else: raise Exception(f"Failed to download FASTA file for PDB ID {pdb_id}") def filter_sequences(self, chain_id: str) -> List[SeqRecord]: return list(filter(lambda x: f"Chain {chain_id}" in x.description, self.read_fasta())) def read_fasta(self) -> SeqRecord: fasta_file = self.pdb_file.parent / f"{self.pid}.fasta" if not fasta_file.exists(): self.download_fasta() return SeqIO.parse(fasta_file.as_posix(), "fasta") @dataclass class PDBAlign: template_file: Path target_file: Path out_file: Path = field(default=None) # 输出文件路径 def align(self): cmd.reinitialize() # 首先，加载模板结构 cmd.load(self.template_file.as_posix(), "template") # 加载并对齐所有目标结构 cmd.load(self.target_file.as_posix(), "target") cmd.align("target", "template") return cmd.get_pdbstr('target') def main(PDB_ID, PDB_file_path): # 示例 # PDB_ID = '5sws' # PDB_file_path = Path(f'{PDB_ID}.pdb') analyzer = PDBAnalyzer(PDB_file_path) sequences = analyzer.extract_sequences(missing_char='-') # 或者 'X', 或者 '' print(f'Residues info for {PDB_ID}: \n',sequences) missing_info = analyzer.extract_sequences_info() print(f'Missing residues info for {PDB_ID}:\n {missing_info}') # 示例: 使用biopython提取A链（将会保留HETATM） chain_extractor = analyzer.extract_chain('A') # 假设要提取的链ID是 'A' chain_extractor.save('biopython_extracted_chain_A.pdb') # 保存为PDB文件 # 示例: 使用biopandas提取A链（将不会保留HETATM） chain_extractor = analyzer.split_chain('A') # 假设要提取的链ID是 'A' chain_extractor.to_pdb('biopandas_extracted_chain_A.pdb') # 保存为PDB文件 # A链改B链, 并分割保存为单独文件 analyzer.change_chain_identifier('A', 'B', split=True).to_pdb(f'{PDB_ID}_B.pdb') # 分割所有的链 split_dict = {} for j in analyzer.chain_id_list: fn = Path(f'{PDB_ID}_{j}.pdb') analyzer.split_chain(j).to_pdb(fn.as_posix()) split_dict[j]=fn.read_text() # 修复loop区域 from build_modeller import PDBModeler from modeller import ModellerError mc_dict = {} for mc in missing_info: out_file = f'5sws_{mc}.pdb' analyzer.split_chain(mc).to_pdb(out_file) # get misschain pdb file mc_fasta = analyzer.filter_sequences(mc) # get misschain fasta file if len(mc_fasta) == 1: mc_fasta = mc_fasta[0] out_fasta_file = Path(f'{PDB_ID}_{mc}.fasta') analyzer.write_seq_to_fasta_single_line(mc_fasta, out_fasta_file) print(f'>{mc_fasta.description}') print(mc_fasta.seq) modeller = PDBModeler(PDB_file_path, out_fasta_file, Path('.'), mc, 1, 'refine.very_fast') try: modeller_results = modeller.make_model() except ModellerError: print(f'Failed to build model for chain {mc}') print(f'No loops detected in {out_fasta_file.name}') print(f'may pdb file sequence is not correct') continue except Exception as e: raise e print(f'Model files: {[file.name for file in modeller_results]}') # change id to original for i in modeller_results: manalyzer = PDBAnalyzer(i) manalyzer.change_chain_identifier('A', mc, split=False).to_pdb(i) if len(modeller_results) == 1: # use pymol to align aligner = PDBAlign(PDB_file_path, modeller_results[0],Path(f'{PDB_ID}_merge_model.pdb')) pdbstr = aligner.align() mc_dict[mc] = pdbstr else: print('more than one model file, please set num_loop to 1') else: raise ValueError(f'only can fix one chain content: {mc_fasta}') # 使用示例 split_dict.update(mc_dict) # 更新 split_dict import_and_merge_pdb_strings(split_dict, "merged_object", f'{PDB_ID}.modellerfix.pdb') def import_and_merge_pdb_strings(pdb_strings, merged_object_name, output_file): """ 从 PDB 字符串导入多个对象，将它们合并，并保存为一个 PDB 文件。 :param pdb_strings: 字典，键为对象名称，值为 PDB 数据字符串。 :param merged_object_name: 合并后的对象名称。 :param output_file: 输出文件的路径。 """ # 初始化 PyMOL（如果在脚本或非交互式环境中使用） cmd.reinitialize() # 从字符串导入每个对象 for chain_id, pdb_str in pdb_strings.items(): object_name = f"chain_{chain_id}" # 创建唯一的对象名称 cmd.read_pdbstr(pdb_str, object_name) # 合并所有对象为一个 object_names = [f"chain_{chain_id}" for chain_id in pdb_strings.keys()] cmd.create(merged_object_name, ' or '.join(object_names)) # 保存合并后的对象 cmd.save(output_file, merged_object_name) if __name__ == "__main__": # import argparse # parser = argparse.ArgumentParser(description="Build model by Modeller") # parser.add_argument("-s", "--structure", help="Structure file") # parser.add_argument("-o", "--outdir", help="Output directory") # parser.add_argument("-f", "--fasta", help="Fasta file") # parser.add_argument("-n", "--num_loop", help="Number of loop model") # parser.add_argument("-m", "--md_level", help="MD level") # parser.add_argument("-c", "--chain", help="Chain ID") # args = parser.parse_args() pdbfiles = [i for i in Path('../PDBfile').glob('*.pdb')] for i in pdbfiles: PDB_file_path = i PDB_ID = i.stem analyzer = PDBAnalyzer(PDB_file_path) sequences = analyzer.extract_sequences(missing_char='-') # 或者 'X', 或者 '' print(f'Residues info for {PDB_ID}: \n',sequences) missing_info = analyzer.extract_sequences_info() print(f'Missing residues info for {PDB_ID}:\n {missing_info}') # 示例: 使用biopython提取A链（将会保留HETATM） chain_extractor = analyzer.extract_chain('A') # 假设要提取的链ID是 'A' chain_extractor.save('biopython_extracted_chain_A.pdb') # 保存为PDB文件 # 示例: 使用biopandas提取A链（将不会保留HETATM） chain_extractor = analyzer.split_chain('A') # 假设要提取的链ID是 'A' chain_extractor.to_pdb('biopandas_extracted_chain_A.pdb') # 保存为PDB文件 # A链改B链, 并分割保存为单独文件 analyzer.change_chain_identifier('A', 'B', split=True).to_pdb(f'{PDB_ID}_B.pdb') # 分割所有的链 split_dict = {} for j in analyzer.chain_id_list: fn = Path(f'{PDB_ID}_{j}.pdb') analyzer.split_chain(j).to_pdb(fn.as_posix()) split_dict[j]=fn.read_text() # 修复loop区域 from build_modeller import PDBModeler from modeller import ModellerError mc_dict = {} for mc in missing_info: out_file = f'5sws_{mc}.pdb' analyzer.split_chain(mc).to_pdb(out_file) # get misschain pdb file mc_fasta = analyzer.filter_sequences(mc) # get misschain fasta file if len(mc_fasta) == 1: mc_fasta = mc_fasta[0] out_fasta_file = Path(f'{PDB_ID}_{mc}.fasta') analyzer.write_seq_to_fasta_single_line(mc_fasta, out_fasta_file) print(f'>{mc_fasta.description}') print(mc_fasta.seq) modeller = PDBModeler(PDB_file_path, out_fasta_file, Path('.'), mc, 1, 'refine.very_fast') try: modeller_results = modeller.make_model() except ModellerError: print(f'Failed to build model for chain {mc}') print(f'No loops detected in {out_fasta_file.name}') print(f'may pdb file sequence is not correct') continue except Exception as e: raise e print(f'Model files: {[file.name for file in modeller_results]}') # change id to original for i in modeller_results: manalyzer = PDBAnalyzer(i) manalyzer.change_chain_identifier('A', mc, split=False).to_pdb(i) if len(modeller_results) == 1: # use pymol to align aligner = PDBAlign(PDB_file_path, modeller_results[0],Path(f'{PDB_ID}_merge_model.pdb')) pdbstr = aligner.align() mc_dict[mc] = pdbstr else: print('more than one model file, please set num_loop to 1') elif len(mc_fasta) == 0: continue else: raise ValueError(f'only can fix one chain content: {mc_fasta}') # 使用示例 split_dict.update(mc_dict) # 更新 split_dict import_and_merge_pdb_strings(split_dict, "merged_object", f'{PDB_ID}.modellerfix.pdb')