<jupyter_start><jupyter_text>想定1. pdbをダウンロード2. correct resnumで残基番号を修正3. 残基番号をもとに配列を生成．このときmissing residueは'-'か何かにする4. fastaと生成した配列をアライメント<jupyter_code>target = '4PNE_B'
pdb_id = target.split('_')[0]
chain = target.split('_')[1]
fasta_path = (fasta_dir / target).with_suffix('.fasta')
pdb_path = (native_dir / pdb_id).with_suffix('.pdb')
fasta_seq = ReadSeq.fasta2seq(fasta_path)
correct_pdb_path = (correct_pdb_dir / pdb_id).with_suffix('.pdb')
mol = parsePDB(pdb_path, chain=chain)
new_resnum = correct_renums(mol)
mol.setResnums(new_resnum)
writePDB(str(correct_pdb_path), mol)
fasta_seq
pdb_seq, pdb_resnums = ReadSeq.pdb2seq(correct_pdb_path, chain, read_HETATM=True, insert_gap=False)
pdb_resnums
pdb_seq
seq_len = pdb_resnums[-1] - pdb_resnums[0] + 1
seq_array = np.array(['-'] * seq_len)
seq_array
seq_indices = pdb_resnums - pdb_resnums[0]
seq_indices
seq_array[seq_indices] = list(pdb_seq)
seq_array
pdb_seq_correct = ''.join(seq_array)
pdb_seq_correct
pairwise2.align.globalms(fasta_seq, pdb_seq_correct, 5, -1, -2, -1)
def align_fasta_pdb(target):
    pdb_id = target.split('_')[0]
    chain = target.split('_')[1]
    fasta_path = (fasta_dir / target).with_suffix('.fasta')
    pdb_path = (native_dir / pdb_id).with_suffix('.pdb')
    correct_pdb_path = (correct_pdb_dir / target).with_suffix('.pdb')
    fasta_seq = ReadSeq.fasta2seq(fasta_path)
    mol = parsePDB(pdb_path, chain=chain)
    new_resnums = correct_renums(mol)
    mol.setResnums(new_resnums)
    pdb_seq, pdb_resnum = ReadSeq.mol2seq(mol, insert_gap=False)
    seq_len = np.max(pdb_resnum) - np.min(pdb_resnum) + 1
    seq_array = np.array(['-'] * seq_len)
    seq_indices = pdb_resnum - pdb_resnum[0]
    try:
        seq_array[seq_indices] = list(pdb_seq)
    except IndexError as e:
        print(target, e)
        print(fasta_seq)
        print(pdb_seq)
        print(pdb_resnum)
        print(seq_indices)
        return None
    pdb_seq_correct = ''.join(seq_array)
    alignment_list = pairwise2.align.globalms(fasta_seq, pdb_seq_correct, 5, -1, -2, -1)
    return alignment_list

for fasta in fasta_dir.glob('*.fasta'):
    target = fasta.stem
    alignment_list = align_fasta_pdb(target)
    # try:
    #     alignment_list = align_fasta_pdb(target)
    # except IndexError as e:
    #     print('target:', target, e)
    if len(alignment_list) > 1:
        print(target)
        for alignment in alignment_list:
            print(alignment)
        print()
from functools import reduce

def get_native_target(target):
    pdb_id = target.split('_')[0]
    chain = target.split('_')[1]
    fasta_path = (fasta_dir / target).with_suffix('.fasta')
    pdb_path = (native_dir / pdb_id).with_suffix('.pdb')
    correct_pdb_path = (correct_pdb_dir / target).with_suffix('.pdb')
    fasta_seq = ReadSeq.fasta2seq(fasta_path)
    mol = parsePDB(pdb_path, chain=chain)
    new_resnums = correct_renums(mol)
    mol.setResnums(new_resnums)
    # TODO ここでresnumsをmolのresnumsの順番から変更してしまっているのでバグが出る
    pdb_seq, pdb_resnums = ReadSeq.mol2seq(mol)
    # pdb_seq, pdb_resnum = ReadSeq.mol2seq(mol, insert_gap=False)
    # if pdb_resnum[0] != np.min(pdb_resnum):
    #     print('different!!!', pdb_resnum[0], np.min(pdb_resnum))
    #     del_indices = np.where(pdb_resnum < pdb_resnum[0])[0]
    #     pdb_seq = ''.join(np.delete(list(pdb_seq), del_indices, 0))
    #     pdb_resnum = np.delete(pdb_resnum, del_indices, 0)
    # if pdb_resnum[0] != np.min(pdb_resnum):
    #     print('different again!!!', pdb_resnum[0], np.min(pdb_resnum))
    # seq_len = np.max(pdb_resnum) - np.min(pdb_resnum) + 1
    # seq_array = np.array(['-'] * seq_len)
    # seq_indices = pdb_resnum - pdb_resnum[0]
    # seq_array[seq_indices] = list(pdb_seq)
    # pdb_seq_correct = ''.join(seq_array)
    align_fseq, align_pseq, align_findices, align_pindices = AlignSeq.align_fasta_and_pdb(fasta_seq, pdb_seq, alignment_percentage_threshold=0.7)
    # TODO 単純にpdb_seq_correctとmol.getResindices()の順序が一部変更されているのでalign_pindicesをそのまま使うとずれが生じる
    # sel_mol = mol.select('resindex {}'.format(reduce(lambda a, b: str(a) + ' ' + str(b), align_pindices)))
    sel_pdb_resnums = pdb_resnums[align_pindices]
    sel_mol = mol.select('resnum {}'.format(reduce(lambda a, b: str(a) + ' ' + str(b), sel_pdb_resnums)))
    assert sel_mol is not None
    # pdb_resnum = sel_mol.getResnums()
    pdb_resnum_start, pdb_resnum_end = str(sel_pdb_resnums[0]), str(sel_pdb_resnums[-1])
    _, scop_resnum_start, scop_resnum_end = read_fasta(fasta_path)
    if pdb_resnum_start != scop_resnum_start or pdb_resnum_end != scop_resnum_end:
        print('The residue number of pdb and scop are different')
        print('pdb start: {}, scop start: {}, pdb end: {}, scop end:{}'
                .format(pdb_resnum_start, scop_resnum_start, pdb_resnum_end, scop_resnum_end))
    fasta_resnum = align_findices + 1
    # TODO ここもalign_pindicesをそのまま使うとずれが生じる
    convert_resnum_dict = dict(zip(sel_pdb_resnums, fasta_resnum))
    # TODO ResindicesではなくResnumsを使う
    new_resnum = [convert_resnum_dict[resnum] for resnum in sel_mol.getResnums()]
    sel_mol.setResnums(new_resnum)
    writePDB(str(correct_pdb_path), sel_mol)
for fasta in fasta_dir.glob('*.fasta'):
    target = fasta.stem
    print(target)
    get_native_target(target)
    # try:
    #     get_native_target(target)
    # except ValueError as e:
    #     print(e)<jupyter_output>4EFH_A
2UUB_E
4WE5_A
5V93_o
1HV5_A
2V2F_F
4HRR_H
5XX9_B
1IFV_A
1W3B_A
3DMI_A
2F2H_A
1DE4_A
4ZPL_A
6KNA_A
3RHH_C
2WP2_B
5CZJ_B
1F2E_A
2BWJ_A
1AYZ_A
1I0D_A
3JAI_SS
4PNE_B
2MR9_A
1P3C_A
2I2Y_A
6K9F_B
6O8W_g
2D9Y_A
2VKE_A
2DS4_A
1OHE_A
3O1D_A
1RSY_A
6CD6_A
1VLK_A
5GJH_A
1I4K_A
4IOP_B
2DKF_A
5DE0_C
1WLU_A
3GFO_A
5JP5_A
1EA9_C
4RM4_A
6SPF_T
5JT8_A
2CNW_D
2HSZ_A
1HTA_A
4G12_A
1TW4_A
1ST9_A
4NBU_B
6KWQ_A
1YR0_A
3M0G_A
1MZE_A
2RKG_A
1V2F_A
1XWN_A
1XZP_A
3JB4_B
1S3J_A
1KWF_A
1OH0_A
1UTI_A
1IRU_I
2GH9_A
5Y4E_A
4NYN_A
1EXR_A
2P6R_A
2UUB_Q
3HRX_A
1UB7_A
2CYB_B
4V52_BC
1KA9_F
2B5A_A
1W8O_A
2V90_C
The residue number of pdb and scop are different
pdb start: 247, scop start: 247, pdb end: 338, scop end:339
4KE6_E
1RE5_A
3CRN_A
3TJ3_A
2YV3_A
2C21_A
1V70_A
1BQU_A
1G4I_A
1QNI_A
1K28_A
4M1Q_B
3GM8_A
2ZJR_H
1RT8_A
1RHF_A<jupyter_text>Test the extracted pdb sequence matches to the fasta sequence for all targets<jupyter_code>cut_pdb_dir = Path('correct_pdb')

def test_match_from_seq(fasta_seq, pdb_seq, pdb_resnum):
    fasta_seq_array = np.array(list(fasta_seq))
    pdb_seq_array = np.copy(fasta_seq_array)
    pdb_seq_array[pdb_resnum - 1] = list(pdb_seq)
    num_diff = np.count_nonzero(fasta_seq_array != pdb_seq_array)
    num_missing = len(fasta_seq) - len(pdb_seq)
    return num_diff, num_missing

def test_match(target):
    chain = target.split('_')[1]
    pdb_path = (cut_pdb_dir / target).with_suffix('.pdb')
    fasta_path = (fasta_dir / target).with_suffix('.fasta')
    fasta_seq = ReadSeq.fasta2seq(fasta_path)
    pdb_seq, pdb_resnum = ReadSeq.pdb2seq(pdb_path, chain=chain, read_HETATM=True, insert_gap=False)

    # pdb_resnum_resname_dict = dict(zip(pdb_resnum, pdb_seq))
    # missing = 0
    # diff = 0
    # for i, r in enumerate(fasta_seq):
    #     fasta_resnum = i + 1
    #     if fasta_resnum in pdb_resnum_resname_dict:
    #         if r != pdb_resnum_resname_dict[fasta_resnum]:
    #             diff += 1
    #     else:
    #         missing += 1
    print('length:', len(fasta_seq))
    # print('diff:', diff)
    # print('missing:', missing)

    diff, missing = test_match_from_seq(fasta_seq, pdb_seq, pdb_resnum)
    print('diff:', diff)
    print('missing:', missing)


for fasta in fasta_dir.glob('*.fasta'):
    target = fasta.stem
    print(target)
    test_match(target)
    print()
target = '2V2F_F'
chain = target.split('_')[1]
pdb_path = (cut_pdb_dir / target).with_suffix('.pdb')
fasta_path = (fasta_dir / target).with_suffix('.fasta')
fasta_seq = ReadSeq.fasta2seq(fasta_path)
pdb_seq, pdb_resnum = ReadSeq.pdb2seq(pdb_path, chain=chain, read_HETATM=True, insert_gap=False)
fasta_seq
pdb_seq
pdb_resnum
pdb_seq_insert = np.array(['-'] * len(fasta_seq))
pdb_seq_insert[pdb_resnum - 1] = list(pdb_seq)
pdb_seq_insert
def test_match(fasta_seq, pdb_seq, pdb_resnum):
    fasta_seq_array = np.array(list(fasta_seq))
    pdb_seq_array = np.copy(fasta_seq_array)
    pdb_seq_array[pdb_resnum - 1] = list(pdb_seq)
    num_diff = np.count_nonzero(fasta_seq_array != pdb_seq_array)
    num_missing = len(fasta_seq) - len(pdb_seq)
    return num_diff, num_missing
target = '3M0G_A'
chain = target.split('_')[1]
pdb_path = (cut_pdb_dir / target).with_suffix('.pdb')
fasta_path = (fasta_dir / target).with_suffix('.fasta')
fasta_seq = ReadSeq.fasta2seq(fasta_path)
pdb_seq, pdb_resnum = ReadSeq.pdb2seq(pdb_path, chain=chain, read_HETATM=True, insert_gap=True)
pdb_seq
fasta_seq
AlignSeq.align_seq(fasta_seq, pdb_seq)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11516726943942134, 'fraction_numerical': 0.027689873417721517, 'mean_word_length': 4.486050836949783, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '208248', 'license': 'MIT', 'n_tokens_mistral': 4137, 'n_tokens_neox': 4019, 'n_words': 622, 'path': 'src/test/align_fasta_pdb/cut_native.ipynb', 'repo_name': 'yutake27/HMDM', 'size': '8848'}