<jupyter_start><jupyter_text>Parse FDA approved drugs from DrugBank<jupyter_code>if not os.path.isfile('../data/processed/fda.csv'):
    # Parse DrugBank XML file.
    ns = '{http://www.drugbank.ca}'
    tree = ET.parse('../data/external/drugbank_517.xml')
    rows = [(drug.findtext(ns + 'drugbank-id[@primary="true"]'),
             drug.findtext(ns + 'name'),
             '|'.join([group.text for group in
                       drug.findall(f'{ns}groups/{ns}group')]),
             '|'.join([code.get('code') for code in
                       drug.findall(f'{ns}atc-codes/{ns}atc-code')]),
             drug.findtext(f'{ns}calculated-properties/'
                           f'{ns}property[{ns}kind="SMILES"]/{ns}value'),
             drug.findtext(f'{ns}experimental-properties/'
                           f'{ns}property[{ns}kind="logP"]/{ns}value'))
            for drug in tree.getroot()]

    approved_drugs = (pd.DataFrame(rows, columns=['drugbank_id', 'name',
                                                  'groups', 'atc_codes',
                                                  'smiles', 'logP'])
                     .dropna(subset=['smiles']))
    # Filter on FDA approved drugs.
    approved_drugs = approved_drugs[approved_drugs['groups']
                                    .str.contains('approved')]
    # Only retain drugs with valid and unique SMILES.
    smiles = []
    for drug_smiles in approved_drugs['smiles']:
        mol = Chem.MolFromSmiles(drug_smiles)
        smiles.append(Chem.MolToSmiles(mol, False)
                      if mol is not None else None)
    approved_drugs['smiles'] = smiles
    approved_drugs = (approved_drugs.dropna(subset=['smiles'])
                      .drop_duplicates('smiles')
                      .reset_index(drop=True))
    approved_drugs.to_csv('../data/processed/fda.csv', index=False)
else:
    approved_drugs = pd.read_csv('../data/processed/fda.csv')<jupyter_output><empty_output><jupyter_text>Epidermis probability for FDA approved drugs<jupyter_code># Load the trained classifier.
classifier = joblib.load('../data/processed/rf.joblib')
# Calculate features using Mordred.
mordred_calculator = Calculator(descriptors, ignore_3D=True)
# Get the original feature labels used during training.
compounds = pd.read_csv('../data/processed/compound_smiles.csv')
mols = compounds['smiles'].apply(Chem.MolFromSmiles)
features_orig = pd.DataFrame(mordred_calculator.pandas(mols)
                             .select_dtypes(exclude='object')
                             .astype(np.float32))
feature_labels = features_orig.columns
if 'epidermis_prob' not in approved_drugs.columns:
    # Generate features for the FDA approved drugs.
    # (Exclude features not encountered during training.)
    mols = approved_drugs['smiles'].apply(Chem.MolFromSmiles)
    features = pd.DataFrame(mordred_calculator.pandas(mols)[feature_labels]
                            .astype(np.float32))

    approved_drugs['epidermis_prob'] = \
        classifier.predict_proba(features.values)[:, 1]
    approved_drugs.to_csv('../data/processed/fda.csv', index=False)
else:
    approved_drugs = pd.read_csv('../data/processed/fda.csv')
(approved_drugs.sort_values('epidermis_prob', ascending=False)
 [['name', 'epidermis_prob']].head(10))<jupyter_output><empty_output><jupyter_text>Epidermis probability for drug biotransformation products<jupyter_code>data_dir = os.path.abspath(os.path.join(
    os.getcwd(), '../data/processed/biotransformer'))
if not os.path.isdir(data_dir):
    os.makedirs(data_dir)
bin_dir = os.path.abspath(os.path.join(os.getcwd(), '../bin/biotransformer'))

# https://bitbucket.org/djoumbou/biotransformer/
def _biotransform(drugbank_id, smiles):
    cmd = f"""cd "{bin_dir}" && \
        java -jar biotransformer-2.0.1.jar \
        -b allHuman \
        -ismi "{smiles}" \
        -k pred \
        -ocsv "{data_dir}/{drugbank_id}.csv" >/dev/null 2>&1"""
    if not os.path.isfile(os.path.join(data_dir, f'{drugbank_id}.csv')):
        subprocess.run(cmd, shell=True)
if not os.path.isfile('../data/processed/fda_biotransformations.csv'):
    # Generate human biotransformation products.
    joblib.Parallel(n_jobs=-1)(
        joblib.delayed(_biotransform)(drugbank_id, smiles)
        for drugbank_id, smiles in tqdm.tqdm(
            zip(approved_drugs['drugbank_id'], approved_drugs['smiles']),
            desc='Biotransformations predicted', total=len(approved_drugs)))
    # Read all biotransformation files for each drug.
    biotransformations = []
    for drugbank_id in approved_drugs['drugbank_id']:
        filename = os.path.join(data_dir, f'{drugbank_id}.csv')
        if os.path.isfile(filename):
            biotransform = pd.read_csv(filename, usecols=['SMILES'])
            biotransform['drugbank_id'] = drugbank_id
            biotransformations.append(biotransform)
    biotransformations = (pd.merge(pd.concat(biotransformations),
                                   approved_drugs[['drugbank_id',
                                                   'atc_codes']],
                                   on='drugbank_id')
                       .copy().dropna()
                       .rename(columns={'SMILES': 'smiles'}))
    # Only retain biotransformations with valid and unique SMILES.
    smiles = []
    for biotransform_smiles in biotransformations['smiles']:
        mol = Chem.MolFromSmiles(biotransform_smiles)
        smiles.append(Chem.MolToSmiles(mol, False)
                      if mol is not None else None)
    biotransformations['smiles'] = smiles
    biotransformations = (biotransformations.dropna(subset=['smiles'])
                          .drop_duplicates('smiles')
                          .reset_index(drop=True))
    biotransformations.to_csv('../data/processed/fda_biotransformations.csv',
                              index=False)
else:
    biotransformations = pd.read_csv(
        '../data/processed/fda_biotransformations.csv')
if 'epidermis_prob' not in biotransformations.columns:
    # Generate features for the biotransformations.
    # (Exclude features not encountered during training.)
    mols = biotransformations['smiles'].apply(Chem.MolFromSmiles)
    features = pd.DataFrame(mordred_calculator.pandas(mols)[feature_labels]
                            .astype(np.float32))
    
    biotransformations['epidermis_prob'] = \
        classifier.predict_proba(features.values)[:, 1]
    biotransformations.to_csv('../data/processed/fda_biotransformations.csv',
                              index=False)
else:
    biotransformations = pd.read_csv(
        '../data/processed/fda_biotransformations.csv')<jupyter_output><empty_output><jupyter_text>Plotting<jupyter_code>atc_map = {'A': 'Alimentary tract and metabolism',
           'B': 'Blood and blood forming organs',
           'C': 'Cardiovascular system',
           'D': 'Dermatologicals',
           'G': 'Genito-urinary system and sex hormones',
           'H': 'Systemic hormonal preparations, excluding sex hormones and insulins',
           'J': 'Antiinfectives for systemic use',
           'L': 'Antineoplastic and immunomodulating agents',
           'M': 'Musculo-skeletal system',
           'N': 'Nervous system',
           'P': 'Antiparasitic products, insecticides and repellents',
           'R': 'Respiratory system',
           'S': 'Sensory organs',
           'V': 'Various'}

approved_drugs['type'] = 'FDA approved drugs'
biotransformations['type'] = 'Human biotransformations'
combined = pd.concat([approved_drugs, biotransformations], ignore_index=True)
combined['atc_codes'] = combined['atc_codes'].str.split('|')
combined = combined.explode('atc_codes')
combined['atc_level1'] = combined['atc_codes'].str[:1].map(atc_map)
width = 7
height = width / 1.618
fig, ax = plt.subplots(figsize=(width, height))

sns.kdeplot(data=combined, x='epidermis_prob', hue='type', ax=ax,
            common_norm=False, fill=True)

ax.set_xlim(0, 1.01)
ax.xaxis.set_major_formatter(mticker.PercentFormatter(xmax=1))

ax.set_xlabel('Predicted probability to observe on epidermis')

ax.legend_.set_title('')

sns.despine()

plt.savefig('fda_biotransformations_predictions.png', dpi=300,
            bbox_inches='tight')
plt.show()
plt.close()
width = 7
height = width / 1.618
fig, ax = plt.subplots(figsize=(width, height * 2))

combined_atc = combined[~combined['atc_level1'].isin(['', 'Dermatologicals'])]
order = (combined_atc.groupby(['atc_level1']).median()
         .sort_values(by='epidermis_prob', ascending=False).index)

sns.violinplot(x='epidermis_prob', y='atc_level1', hue='type',
               data=combined_atc, order=order, split=True, orient='h', ax=ax)

ax.set_xlim(0, 1.01)
ax.xaxis.set_major_formatter(mticker.PercentFormatter(xmax=1))

ax.set_xlabel('Predicted probability to observe on epidermis')
ax.set_ylabel('ATC level 1')

ax.legend(loc='upper center', bbox_to_anchor=(0.5, 1.06), ncol=2)

sns.despine()

plt.savefig('fda_biotransformations_predictions_atc.png', dpi=300,
            bbox_inches='tight')
plt.show()
plt.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13655902127190567, 'fraction_numerical': 0.006502810536757412, 'mean_word_length': 2.983318700614574, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 19, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '602578', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3146, 'n_tokens_neox': 2880, 'n_words': 477, 'path': 'notebooks/fda_biotransformations.ipynb', 'repo_name': 'bittremieux/drugs_epidermis', 'size': '9073'}