from Bio import PDB
from Bio.PDB.Polypeptide import PPBuilder
import sys

def parse_pdb(pdb_file):
    """PDBファイルを解析して構造オブジェクトを取得"""
    parser = PDB.PDBParser(QUIET=True)
    structure = parser.get_structure('model', pdb_file)
    return structure

def get_chain_info(structure):
    """各チェインの情報（配列、長さ、残基番号）を取得"""
    ppb = PPBuilder()
    chain_info = {}
    for model in structure:
        for chain in model:
            seq = ""
            res_nums = []
            for pp in ppb.build_peptides(chain):
                seq += pp.get_sequence()
                res_nums += [res.get_id()[1] for res in pp]
            if not seq:
                # PPBuilderが失敗した場合、手動で配列を取得
                seq = ''.join([PDB.Polypeptide.three_to_one(res.get_resname()) 
                              for res in chain if PDB.is_aa(res)])
                res_nums = [res.get_id()[1] for res in chain if PDB.is_aa(res)]
            if seq:
                chain_id = chain.get_id()
                chain_info[chain_id] = {
                    'sequence': str(seq),
                    'length': len(seq),
                    'res_nums': res_nums,
                    'start': min(res_nums) if res_nums else None,
                    'end': max(res_nums) if res_nums else None,
                }
    return chain_info

def filter_chains(chain_info):
    """短いチェインの削除と重複チェインの削除"""
    # 長さが10以下のチェインを削除
    filtered = {cid: info for cid, info in chain_info.items() if info['length'] > 10}
    
    # チェインをグループ化（例えば、Aa, Ab, Ac は同じグループA）
    groups = {}
    for cid, info in filtered.items():
        group_key = cid[0]  # 最初の文字でグループ化
        if group_key not in groups:
            groups[group_key] = []
        groups[group_key].append((cid, info))
    
    # 各グループ内で重複チェインを削除
    final_chains = set()
    for group, chains in groups.items():
        # チェインを長さの降順でソート
        chains_sorted = sorted(chains, key=lambda x: x[1]['length'], reverse=True)
        kept = []
        for cid, info in chains_sorted:
            overlap = False
            for kept_cid, kept_info in kept:
                # 残基番号の範囲が重複するかチェック
                if info['start'] <= kept_info['end'] and info['end'] >= kept_info['start']:
                    overlap = True
                    break
            if not overlap:
                kept.append((cid, info))
        # 保持するチェインIDを追加
        for cid, info in kept:
            final_chains.add(cid)
    return final_chains

def write_clean_pdb(pdb_file, output_file, kept_chains):
    """クリーンアップされたチェインのみを含む新しいPDBファイルを作成"""
    parser = PDB.PDBParser(QUIET=True)
    io = PDB.PDBIO()
    structure = parser.get_structure('model', pdb_file)
    
    class ChainSelect(PDB.Select):
        def accept_chain(self, chain):
            return chain.get_id() in kept_chains
    
    io.set_structure(structure)
    io.save(output_file, select=ChainSelect())

def main():
    if len(sys.argv) != 3:
        print("使用法: python filter_modelangelo_chains.py input.pdb output_cleaned.pdb")
        sys.exit(1)
    
    input_pdb = sys.argv[1]
    output_pdb = sys.argv[2]
    
    print("PDBファイルを解析中...")
    structure = parse_pdb(input_pdb)
    
    print("チェイン情報を取得中...")
    chain_info = get_chain_info(structure)
    
    print("チェインをフィルタリング中...")
    kept_chains = filter_chains(chain_info)
    
    print(f"保持するチェイン: {', '.join(kept_chains)}")
    
    print("クリーンアップされたPDBファイルを作成中...")
    write_clean_pdb(input_pdb, output_pdb, kept_chains)
    
    print(f"フィルタリング完了。クリーンアップされたPDBは {output_pdb} に保存されました。")

if __name__ == "__main__":
    main()