<jupyter_start><jupyter_text>Create commands for running NetMHCpan4.0<jupyter_code>cmd_template = '-BA -p peptides_8-10.pep -a {allele} -xls -xlsfile {allele_file}'
cmds = []
for allele in hla_alleles['allele'].values:
    allele_file = 'netmhc_preds/%s_preds.xls' % (allele.replace(':', ''))
    # Check if allele file exists.
    if os.path.exists(allele_file):
        continue
    cmd = cmd_template.format(
        allele=allele,
        allele_file=allele_file,
    )
    cmds.append(cmd)
print('# Commands: ', len(cmds))

with open('netmhc_args.txt', 'w') as f:
    for cmd in cmds:
        f.write(cmd+'\n')
# For NetMHCpan4.1:

cmd_template = '-BA -p peptides_8-10.pep -a {allele} -xls -xlsfile {allele_file}'
cmds = []
for allele in hla_alleles['allele'].values:
    allele_file = 'netmhc-4.1_preds/%s_preds.xls' % (allele.replace(':', ''))
    # Check if allele file exists.
    if os.path.exists(allele_file):
        continue
    cmd = cmd_template.format(
        allele=allele,
        allele_file=allele_file,
    )
    cmds.append(cmd)
print('# Commands: ', len(cmds))

with open('netmhc-4.1_args.txt', 'w') as f:
    for cmd in cmds:
        f.write(cmd+'\n')<jupyter_output># Commands:  230<jupyter_text>Command to run:    ```cat netmhc_args.txt | xargs -P 40 -d '\n' -n 1 ./netMHCpan-4.0/netMHCpancat netmhc-4.1_args.txt | xargs -P 35 -d '\n' -n 1 ./netMHCpan-4.1/netMHCpan```  Load NetMHCpan4.0 predictions<jupyter_code>dfs = []
for allele in hla_alleles['allele'].values:
    df = pd.read_csv(
        './netmhc_preds/%s_preds.xls' % allele.replace(':', ''),
        delimiter='\t',
        skiprows=[0],
    )
    df['genotype'] = allele
    df = df.drop(columns=['Pos', 'ID', 'core', 'icore', 'Rank', 'Ave', 'NB'])
    dfs.append(df)
netmhc1_data = pd.concat(dfs)
netmhc1_data['sequence_length'] = [len(x) for x in netmhc1_data['Peptide'].values]
netmhc1_data['loci'] = [x[:5] for x in netmhc1_data['genotype'].values]
netmhc1_data
a74 = netmhc1_data.loc[netmhc1_data['genotype'].str.contains('HLA-A74')].groupby('Peptide').agg('mean').reset_index() #['mean', 'count'])
a74['loci'] = 'HLA-A'
a74['genotype'] = 'HLA-A74'
#a74

c17 = netmhc1_data.loc[netmhc1_data['genotype'].str.contains('HLA-C17')].groupby('Peptide').agg('mean').reset_index()
c17['loci'] = 'HLA-C'
c17['genotype'] = 'HLA-C17'
#c17

c18 = netmhc1_data.loc[netmhc1_data['genotype'].str.contains('HLA-C18')].groupby('Peptide').agg('mean').reset_index()
c18['loci'] = 'HLA-C'
c18['genotype'] = 'HLA-C18'
#c18
pd.concat([netmhc1_data, a74, c17, c18], sort=False)
data_pivot = pd.concat([netmhc1_data, a74, c17, c18], sort=False).pivot_table(
    index='Peptide',
    columns=['loci', 'genotype'],
    values='1-log50k',
)
data_pivot.to_pickle('mhc1_haplotype_netmhc_pred_affinity_pivot.pkl.gz', protocol=2)
data_pivot<jupyter_output><empty_output><jupyter_text>Load NetMHCpan4.1 Predictions<jupyter_code>dfs = []
for allele in hla_alleles['allele'].values:
    df = pd.read_csv(
        './netmhc-4.1_preds/%s_preds.xls' % allele.replace(':', ''),
        delimiter='\t',
        skiprows=[0],
    )
    df['genotype'] = allele
    df = df.drop(columns=['Pos', 'ID', 'core', 'icore', 'Ave', 'NB'])
    dfs.append(df)
netmhc41_data = pd.concat(dfs)
netmhc41_data['sequence_length'] = [len(x) for x in netmhc41_data['Peptide'].values]
netmhc41_data['loci'] = [x[:5] for x in netmhc41_data['genotype'].values]
netmhc41_data['BA_nM'] = 50000**(1-netmhc41_data['BA-score'])
netmhc41_data
a74 = netmhc41_data.loc[netmhc41_data['genotype'].str.contains('HLA-A74')].groupby('Peptide').agg('mean').reset_index() #['mean', 'count'])
a74['loci'] = 'HLA-A'
a74['genotype'] = 'HLA-A74'
#a74

c17 = netmhc41_data.loc[netmhc41_data['genotype'].str.contains('HLA-C17')].groupby('Peptide').agg('mean').reset_index()
c17['loci'] = 'HLA-C'
c17['genotype'] = 'HLA-C17'
#c17

c18 = netmhc41_data.loc[netmhc41_data['genotype'].str.contains('HLA-C18')].groupby('Peptide').agg('mean').reset_index()
c18['loci'] = 'HLA-C'
c18['genotype'] = 'HLA-C18'
#c18
data_pivot = pd.concat([netmhc41_data, a74, c17, c18], sort=False).pivot_table(
    index='Peptide',
    columns=['loci', 'genotype'],
    values='BA-score',
)
data_pivot.to_pickle('mhc1_haplotype_netmhc-4.1_pred_affinity_pivot.pkl.gz', protocol=2)
data_pivot<jupyter_output><empty_output><jupyter_text>Create data for MHCflurry<jupyter_code># Create dataframe with all MHC/peptide pairs.
a = mhc1_data.copy()
b = hla_alleles.copy()
a['key'] = 0
b['key'] = 0
pmhc_pairs = a.merge(b, how='outer')
pmhc_pairs = pmhc_pairs.drop(columns=['key'])
pmhc_pairs
pmhc_pairs[['allele', 'sequence']].to_csv('mhc1_8-10_haplotype_pairs.csv', index=False, header=False)
pmhc_pairs.rename(columns={'sequence': 'peptide'})[['allele', 'peptide']].to_csv(
    'mhc1_8-10_haplotype_pairs_withheader.csv', index=False, header=True)<jupyter_output><empty_output><jupyter_text>Command to run MHCflurry:```mhcflurry-predict \    --out mhc1_haplotype_preds_mhcflurry.csv \    mhc1_8-10_haplotype_pairs_withheader.csvmhcflurry-predict \    --out mhc1_haplotype_preds_mhcflurry2.0.csv \    mhc1_8-10_haplotype_pairs_withheader.csv```  Load MHCflurry Predictions<jupyter_code>def transform_affinity(x):
    x = np.clip(x, a_min=None, a_max=50000)
    return 1 - np.log(x) / np.log(50000)

# print(transform_affinity(500))
mhcflurry_preds = pd.read_csv('mhc1_haplotype_preds_mhcflurry.csv')

mhcflurry_preds = mhcflurry_preds.rename(columns={'peptide': 'sequence'})
mhcflurry_preds['sequence_length'] = [len(x) for x in mhcflurry_preds['sequence'].values]

# Filter dataframe to alleles in hla_alleles.
# mhcflurry_preds = mhcflurry_preds.merge(hla_alleles, on='allele')

# Add epitope protein data.
# mhcflurry_preds = mhcflurry_preds.merge(mhc1_data[['sequence', 'Protein']], on='sequence')

# Compute logistic-transformed binding affinity.
mhcflurry_preds['transformed_aff'] = [transform_affinity(x) for x in mhcflurry_preds['mhcflurry_affinity'].values]

mhcflurry_preds
df = mhcflurry_preds.copy()
df['loci'] = [x[:5] for x in df['allele'].values]
df = df.drop(columns=['sequence_length', 'mhcflurry_affinity_percentile', 'mhcflurry_processing_score', 'mhcflurry_presentation_score'])
df = df.rename(columns={'sequence': 'Peptide', 'allele': 'genotype'})
df
a74 = df.loc[df['genotype'].str.contains('HLA-A74')].groupby('Peptide').agg('mean').reset_index() #['mean', 'count'])
a74['loci'] = 'HLA-A'
a74['genotype'] = 'HLA-A74'
#a74

c17 = df.loc[df['genotype'].str.contains('HLA-C17')].groupby('Peptide').agg('mean').reset_index()
c17['loci'] = 'HLA-C'
c17['genotype'] = 'HLA-C17'
#c17

c18 = df.loc[df['genotype'].str.contains('HLA-C18')].groupby('Peptide').agg('mean').reset_index()
c18['loci'] = 'HLA-C'
c18['genotype'] = 'HLA-C18'
#c18
pd.concat([df, a74, c17, c18], sort=False)
data_pivot = pd.concat([df, a74, c17, c18], sort=False).pivot_table(
    index='Peptide',
    columns=['loci', 'genotype'],
    values='transformed_aff',
)
data_pivot.to_pickle('mhc1_haplotype_mhcflurry_pred_affinity_pivot.pkl.gz', protocol=2)
data_pivot
df = mhcflurry_preds.copy()
df['loci'] = [x[:5] for x in df['allele'].values]
df = df.drop(columns=['sequence_length', 'mhcflurry_affinity_percentile', 'mhcflurry_processing_score', 'mhcflurry_presentation_score'])
df = df.rename(columns={'sequence': 'Peptide', 'allele': 'genotype'})

# df2 = df.groupby(['Peptide', 'loci']).count().reset_index()[['Peptide', 'loci']]
# df2['genotype'] = 'unknown'
# df2['transformed_aff'] = 0.
# df2['mhcflurry_affinity'] = 0.

# df_with_unknown = pd.concat([df, df2], sort=False)

# data_pivot = df_with_unknown.pivot_table(
#     index='Peptide',
#     columns=['loci', 'genotype'],
#     values='transformed_aff',
# )

data_pivot = df.pivot_table(
    index='Peptide',
    columns=['loci', 'genotype'],
    values='transformed_aff',
)
# data_pivot.to_pickle('mhc1_haplotype_mhcflurry_pred_affinity_pivot.pkl.gz', protocol=2)
data_pivot<jupyter_output><empty_output><jupyter_text>Load MHCflurry 2.0 Predictions<jupyter_code>mhcflurry_preds = pd.read_csv('mhc1_haplotype_preds_mhcflurry2.0.csv')

mhcflurry_preds = mhcflurry_preds.rename(columns={'peptide': 'sequence'})
mhcflurry_preds['sequence_length'] = [len(x) for x in mhcflurry_preds['sequence'].values]

# Filter dataframe to alleles in hla_alleles.
# mhcflurry_preds = mhcflurry_preds.merge(hla_alleles, on='allele')

# Add epitope protein data.
# mhcflurry_preds = mhcflurry_preds.merge(mhc1_data[['sequence', 'Protein']], on='sequence')

# Compute logistic-transformed binding affinity.
mhcflurry_preds['transformed_aff'] = 1 - np.log(mhcflurry_preds['mhcflurry_affinity']) / np.log(50000)

mhcflurry_preds['loci'] = [x[:5] for x in mhcflurry_preds['allele'].values]
mhcflurry_preds.drop(
    columns=['sequence_length', 'mhcflurry_affinity_percentile', 'mhcflurry_processing_score', 'mhcflurry_presentation_score'],
    inplace=True,
)

mhcflurry_preds.rename(columns={'allele': 'genotype', 'sequence': 'Peptide'}, inplace=True)

mhcflurry_preds
df = mhcflurry_preds

a74 = df.loc[df['genotype'].str.contains('HLA-A74')].groupby('Peptide').agg('mean').reset_index() #['mean', 'count'])
a74['loci'] = 'HLA-A'
a74['genotype'] = 'HLA-A74'
#a74

c17 = df.loc[df['genotype'].str.contains('HLA-C17')].groupby('Peptide').agg('mean').reset_index()
c17['loci'] = 'HLA-C'
c17['genotype'] = 'HLA-C17'
#c17

c18 = df.loc[df['genotype'].str.contains('HLA-C18')].groupby('Peptide').agg('mean').reset_index()
c18['loci'] = 'HLA-C'
c18['genotype'] = 'HLA-C18'
#c18

data_pivot = pd.concat([df, a74, c17, c18], sort=False).pivot_table(
    index='Peptide',
    columns=['loci', 'genotype'],
    values='transformed_aff',
)
data_pivot.to_pickle('mhc1_haplotype_mhcflurry2.0_pred_affinity_pivot.pkl.gz', protocol=2)
data_pivot<jupyter_output><empty_output><jupyter_text>Mean Ensemble from NetMHC and MHCflurry<jupyter_code>ens_netmhc1 = netmhc1_data.copy().drop(columns=['sequence_length', 'loci', '1-log50k'])
ens_netmhc1 = ens_netmhc1.rename(columns={'nM': 'netmhc_nM'})
ens_netmhc1
ens_mhcflurry = mhcflurry_preds.copy().drop(
    columns=['mhcflurry_affinity_percentile', 'mhcflurry_processing_score', 'mhcflurry_presentation_score',
             'transformed_aff', 'sequence_length'])
ens_mhcflurry = ens_mhcflurry.rename(
    columns={'allele': 'genotype', 'sequence': 'Peptide', 'mhcflurry_affinity': 'mhcflurry_nM'})
ens_mhcflurry
ens_combined = ens_netmhc1.merge(ens_mhcflurry, on=['genotype', 'Peptide'], how='inner')
ens_combined['loci'] = [x[:5] for x in ens_combined['genotype'].values]
ens_combined['mean_nM'] = (ens_combined['netmhc_nM'] + ens_combined['mhcflurry_nM']) / 2
ens_combined['max_nM'] = ens_combined[['netmhc_nM', 'mhcflurry_nM']].max(axis=1)
ens_combined['mean_transformed'] = [transform_affinity(x) for x in ens_combined['mean_nM'].values]
ens_combined['max_transformed'] = [transform_affinity(x) for x in ens_combined['max_nM'].values]
ens_combined
print('Total pMHC pairs: ', len(ens_combined))
print('# NetMHC binders: ', (ens_combined['netmhc_nM'].values <= 50).sum())
print('# MHCflurry binders: ', (ens_combined['mhcflurry_nM'].values <= 50).sum())
print('# Mean Ens. binders: ', (ens_combined['mean_nM'].values <= 50).sum())
print('# Max Ens. binders: ', (ens_combined['max_nM'].values <= 50).sum())
print(transform_affinity(50))
transform_affinity(31838)
df = ens_combined

a74 = df.loc[df['genotype'].str.contains('HLA-A74')].groupby('Peptide').agg('mean').reset_index() #['mean', 'count'])
a74['loci'] = 'HLA-A'
a74['genotype'] = 'HLA-A74'
#a74

c17 = df.loc[df['genotype'].str.contains('HLA-C17')].groupby('Peptide').agg('mean').reset_index()
c17['loci'] = 'HLA-C'
c17['genotype'] = 'HLA-C17'
#c17

c18 = df.loc[df['genotype'].str.contains('HLA-C18')].groupby('Peptide').agg('mean').reset_index()
c18['loci'] = 'HLA-C'
c18['genotype'] = 'HLA-C18'
#c18

data_pivot = pd.concat([df, a74, c17, c18], sort=False).pivot_table(
    index='Peptide',
    columns=['loci', 'genotype'],
    values='mean_transformed',
)
data_pivot.to_pickle('mhc1_haplotype_mean-ensemble_pivot.pkl.gz', protocol=2)
data_pivot
df = ens_combined

a74 = df.loc[df['genotype'].str.contains('HLA-A74')].groupby('Peptide').agg('mean').reset_index() #['mean', 'count'])
a74['loci'] = 'HLA-A'
a74['genotype'] = 'HLA-A74'
#a74

c17 = df.loc[df['genotype'].str.contains('HLA-C17')].groupby('Peptide').agg('mean').reset_index()
c17['loci'] = 'HLA-C'
c17['genotype'] = 'HLA-C17'
#c17

c18 = df.loc[df['genotype'].str.contains('HLA-C18')].groupby('Peptide').agg('mean').reset_index()
c18['loci'] = 'HLA-C'
c18['genotype'] = 'HLA-C18'
#c18

data_pivot = pd.concat([df, a74, c17, c18], sort=False).pivot_table(
    index='Peptide',
    columns=['loci', 'genotype'],
    values='max_transformed',
)
data_pivot.to_pickle('mhc1_haplotype_max-ensemble_pivot.pkl.gz', protocol=2)
data_pivot<jupyter_output><empty_output><jupyter_text>Mean ensemble of NetMHCpan4.0 and MHCflurry2.0<jupyter_code>netmhc_pivot = pd.read_pickle('mhc1_haplotype_netmhc_pred_affinity_pivot.pkl.gz')
mhcflurry_pivot = pd.read_pickle('mhc1_haplotype_mhcflurry2.0_pred_affinity_pivot.pkl.gz')
assert netmhc_pivot.shape == mhcflurry_pivot.shape
assert set(netmhc_pivot.T.index.values.tolist()) == set(mhcflurry_pivot.T.index.values.tolist())
assert set(netmhc_pivot.index.values.tolist()) == set(mhcflurry_pivot.index.values.tolist())

netmhc_pivot_nm = 50000**(1-netmhc_pivot)
mhcflurry_pivot_nm = 50000**(1-mhcflurry_pivot)

ens_pivot_nm = (netmhc_pivot_nm + mhcflurry_pivot_nm) / 2
ens_pivot = 1 - np.log(ens_pivot_nm) / np.log(50000)

ens_pivot.to_pickle('mhc1_netmhcpan4.0_mhcflurry2.0_ensemble_affinity.pkl.gz', protocol=2)
ens_pivot<jupyter_output><empty_output><jupyter_text>PUFFIN<jupyter_code># Create dataframe with all MHC/peptide pairs.
a = mhc1_data.copy()
b = hla_alleles.copy()
a['key'] = 0
b['key'] = 0
pmhc_pairs = a.merge(b, how='outer')
pmhc_pairs = pmhc_pairs.drop(columns=['key'])
pmhc_pairs
pmhc_pairs['placeholder'] = -1
pmhc_pairs
pmhc_pairs[['sequence', 'placeholder', 'allele']].to_csv(
    'puffin_mhc1_preds/datafile.tsv', sep='\t', index=False, header=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18938480096501809, 'fraction_numerical': 0.0354786064003406, 'mean_word_length': 6.5127931769722816, 'pattern_counts': {'":': 0, '<': 40, '<?xml version=': 0, '>': 36, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '498489', 'license': 'MIT', 'n_tokens_mistral': 6360, 'n_tokens_neox': 5910, 'n_words': 790, 'path': 'scoring/peptide_scoring_mhc1.ipynb', 'repo_name': 'gifford-lab/optivax', 'size': '14093'}