Using PaDELPy to Generate Molecular Fingerprints for Machine Learning-Based QSAR

A fingerprint model is a matrix plus a label. Each row is a molecule, each column is a bit or a count from a defined fingerprint, and one extra column says what the molecule did in the assay. PaDEL builds that matrix. The number you can defend is the score on rows the forest did not train on.

This note computes Substructure fingerprints for a curated HCV NS5B table, drops bits with almost no variance, and fits a random forest. It reports three Matthews correlations: training, five-fold cross-validation on the training split, and a held-out test split. It does not claim that those values are high. Run the notebook and read them.

What PaDELPy is calling

PaDEL-Descriptor is a Java program that calculates descriptors and fingerprints. PaDELPy is a Python wrapper around that program. It saves you from writing the Java command yourself. It does not remove the dependency: a Java runtime still has to be available for the wrapper to call.

The fingerprint catalog used here is the XML set from the dataprofessor/padel tutorial. The install lines below use the notebook shell shortcut. In a terminal, drop the exclamation mark.

# Colab or Jupyter shell shortcuts. Locally, run these without the "!".
!pip install padelpy
!wget https://github.com/dataprofessor/padel/raw/main/fingerprints_xml.zip
!unzip -o fingerprints_xml.zip

Unzipping that archive produces twelve XML files. The names line up with the labels as follows. The map is written out on purpose. Zipping a hand-made label list to glob("*.xml") works only while the directory contains exactly these files, in the order the labels assume. One extra XML file would silently point “Substructure” at the wrong definition.

import glob

xml_files = sorted(glob.glob("*.xml"))
xml_files

fp = {
    "AtomPairs2DCount": "AtomPairs2DFingerprintCount.xml",
    "AtomPairs2D": "AtomPairs2DFingerprinter.xml",
    "EState": "EStateFingerprinter.xml",
    "CDKextended": "ExtendedFingerprinter.xml",
    "CDK": "Fingerprinter.xml",
    "CDKgraphonly": "GraphOnlyFingerprinter.xml",
    "KlekotaRothCount": "KlekotaRothFingerprintCount.xml",
    "KlekotaRoth": "KlekotaRothFingerprinter.xml",
    "MACCS": "MACCSFingerprinter.xml",
    "PubChem": "PubchemFingerprinter.xml",
    "SubstructureCount": "SubstructureFingerprintCount.xml",
    "Substructure": "SubstructureFingerprinter.xml",
}

Count variants (AtomPairs2DCount, KlekotaRothCount, SubstructureCount) store how often a pattern occurs. The others store presence or absence. This tutorial uses Substructure, the presence/absence version.

The table

The file read here is the CSV hosted with that same tutorial:

https://raw.githubusercontent.com/dataprofessor/data/master/HCV_NS5B_Curated.csv

It is a curated HCV NS5B bioactivity table. The same curated set is the one distributed with the chaninlab/hcvpred project. The notebook uses the raw CSV URL because that is the file the code actually loads. The label column is Activity. It is already a class label in the file. This note does not apply a second cutoff.

import pandas as pd

df = pd.read_csv(
    "https://raw.githubusercontent.com/dataprofessor/data/master/HCV_NS5B_Curated.csv"
)
df[["CANONICAL_SMILES", "CMPD_CHEMBLID", "Activity"]].head()

PaDEL wants a structure file, not a data frame. The two-column SMILES file below has no header: SMILES, then the ChEMBL ID, separated by a tab.

df2 = pd.concat([df["CANONICAL_SMILES"], df["CMPD_CHEMBLID"]], axis=1)
df2.to_csv("molecule.smi", sep="\t", index=False, header=False)

From SMILES to the matrix

from padelpy import padeldescriptor

fingerprint = "Substructure"
fingerprint_output_file = fingerprint + ".csv"
fingerprint_descriptortypes = fp[fingerprint]

padeldescriptor(
    mol_dir="molecule.smi",
    d_file=fingerprint_output_file,
    descriptortypes=fingerprint_descriptortypes,
    detectaromaticity=True,
    standardizenitro=True,
    standardizetautomers=True,
    threads=2,
    removesalt=True,
    log=True,
    fingerprints=True,
)

descriptors = pd.read_csv(fingerprint_output_file)

PaDEL writes the molecule name into a Name column and one column per fingerprint bit. The next block refuses to continue if the row count changed. Pairing Activity by position is valid only when every input SMILES produced one output row, in the same order. A skipped molecule would shift every later label.

if len(descriptors) != len(df):
    raise ValueError(
        "PaDEL returned a different number of rows than the input table. "
        "Align Activity on molecule ID before training."
    )

X = descriptors.drop(columns=["Name"])
y = df["Activity"].reset_index(drop=True)

Bits that are on for almost every molecule, or off for almost every molecule, do not separate the classes. VarianceThreshold with a cutoff of 0.1 drops those columns. For a bit that is either 0 or 1, a variance above 0.1 means the bit is on in a substantial fraction of the rows and off in a substantial fraction. The cutoff is the one used in this tutorial. It is aggressive for a sparse fingerprint, and Substructure fingerprints are being asked to survive it.

from sklearn.feature_selection import VarianceThreshold

def remove_low_variance(input_data, threshold=0.1):
    selection = VarianceThreshold(threshold)
    selection.fit(input_data)
    return input_data[input_data.columns[selection.get_support(indices=True)]]

X = remove_low_variance(X, threshold=0.1)

A forest, and three different scores

The split is 80 percent training and 20 percent test, with random_state fixed so the split can be repeated. The forest has 500 trees.

from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import matthews_corrcoef
from sklearn.model_selection import cross_val_score, train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

model = RandomForestClassifier(n_estimators=500, random_state=42)
model.fit(X_train, y_train)

mcc_train = matthews_corrcoef(y_train, model.predict(X_train))
mcc_test = matthews_corrcoef(y_test, model.predict(X_test))

# Refit on five folds of the training split only. The test rows stay untouched.
# cross_val_score clones the classifier, so this does not reuse the fitted trees.
mcc_cv = cross_val_score(
    RandomForestClassifier(n_estimators=500, random_state=42),
    X_train,
    y_train,
    cv=5,
    scoring="matthews_corrcoef",
).mean()

performance_metrics = pd.DataFrame(
    {
        "Model": ["Random Forest"],
        "MCC_Train": [mcc_train],
        "MCC_CV": [mcc_cv],
        "MCC_Test": [mcc_test],
    }
)
performance_metrics

An earlier version of this note built that table twice and referred to mcc_cv without ever computing it. The cross-validation line above is the missing definition. It is the mean Matthews correlation across five folds of the training set.

Read the three numbers in that order, and do not collapse them into “the MCC.”

• Training MCC is the forest scored on rows it has already seen. A large forest can look excellent here and still fail on new molecules.
• Cross-validated MCC refits the same settings five times inside the training split. It is the honest internal number. It is still not an external test.
• Test MCC is the 20 percent that train_test_split held out. One split is a single draw. It is the right place to look, and it is not a confidence interval.

If a fold contains only one class, scikit-learn will warn and that fold’s score is not interpretable. That is a reason to look at the class balance before trusting the mean.

Bottom line

PaDELPy is the step that turns the SMILES column into a fingerprint matrix, provided Java is installed and the XML file is the fingerprint you think it is. The random forest is then an ordinary classifier on that matrix.

The earlier conclusion that “high” Matthews correlations show the model is effective was not backed by a computed cross-validation score. Compute the three numbers. Keep the test set out of the fitting. If you try another fingerprint from the same XML list, change fingerprint and run the same split again, so the comparison is the representation rather than a new random partition.