Molfeat Molecular Featurization

uv pip install molfeat

# Optional extras for specific featurizer types
uv pip install "molfeat[transformer]"   # ChemBERTa, ChemGPT, MolT5
uv pip install "molfeat[dgl]"           # GIN graph neural networks
uv pip install "molfeat[graphormer]"    # Graphormer models
uv pip install "molfeat[fcd]"           # FCD descriptors
uv pip install "molfeat[map4]"          # MAP4 fingerprints
uv pip install "molfeat[all]"           # All dependencies

from molfeat.calc import FPCalculator
from molfeat.trans import MoleculeTransformer

smiles = ["CCO", "CC(=O)O", "c1ccccc1", "CC(C)O"]

# Create fingerprint calculator + transformer
calc = FPCalculator("ecfp", radius=3, fpSize=2048)
transformer = MoleculeTransformer(calc, n_jobs=-1)

# Featurize batch in parallel
features = transformer(smiles)
print(f"Shape: {features.shape}")  # (4, 2048)

# Save configuration for reproducibility
transformer.to_state_yaml_file("featurizer_config.yml")

# Save
transformer.to_state_yaml_file("config.yml")
transformer.to_state_json_file("config.json")

# Reload
loaded = MoleculeTransformer.from_state_yaml_file("config.yml")

from molfeat.calc import FPCalculator

# ECFP — most popular, general-purpose
ecfp = FPCalculator("ecfp", radius=3, fpSize=2048)
fp = ecfp("CCO")
print(f"ECFP shape: {fp.shape}")  # (2048,)

# MACCS keys — 167-bit structural keys, fast scaffold similarity
maccs = FPCalculator("maccs")
fp = maccs("c1ccccc1")
print(f"MACCS shape: {fp.shape}")  # (167,)

# Count-based fingerprints (non-binary)
ecfp_count = FPCalculator("ecfp-count", radius=3, fpSize=2048)

# MAP4 — MinHashed atom-pair, efficient for large databases
map4 = FPCalculator("map4")
print(f"MAP4 shape: {map4('CCO').shape}")  # (1024,)

from molfeat.calc import RDKitDescriptors2D, MordredDescriptors

# RDKit 2D — 200+ named properties (MW, logP, TPSA, etc.)
desc2d = RDKitDescriptors2D()
descriptors = desc2d("CCO")
print(f"2D descriptors: {len(descriptors)}")  # 200+
print(f"Feature names: {desc2d.columns[:5]}")

# Mordred — 1800+ comprehensive descriptors
mordred = MordredDescriptors()
descriptors = mordred("c1ccccc1O")
print(f"Mordred descriptors: {len(descriptors)}")  # 1800+

from molfeat.calc import CATSCalculator, USRDescriptors

# CATS — pharmacophore point pair distributions
cats = CATSCalculator(mode="2D", scale="raw")
descriptors = cats("CC(C)Cc1ccc(C)cc1C")
print(f"CATS shape: {descriptors.shape}")  # (21,)

# USR — ultrafast shape recognition
usr = USRDescriptors()
shape = usr("CC(=O)Oc1ccccc1C(=O)O")
print(f"USR shape: {shape.shape}")  # (12,)

from molfeat.trans import MoleculeTransformer, FeatConcat
from molfeat.calc import FPCalculator

smiles = ["CCO", "CC(=O)O", "c1ccccc1", "CC(C)O", "CCCC"]

# Parallel batch processing
transformer = MoleculeTransformer(FPCalculator("ecfp"), n_jobs=-1)
features = transformer(smiles)
print(f"Batch shape: {features.shape}")  # (5, 2048)

# Concatenate multiple featurizers
concat = FeatConcat([
    FPCalculator("maccs"),      # 167 dims
    FPCalculator("ecfp")        # 2048 dims
])
combo_transformer = MoleculeTransformer(concat, n_jobs=-1)
combo_features = combo_transformer(smiles)
print(f"Combined shape: {combo_features.shape}")  # (5, 2215)

# Error-tolerant processing
safe_transformer = MoleculeTransformer(
    FPCalculator("ecfp"), n_jobs=-1,
    ignore_errors=True, verbose=True
)
features = safe_transformer(["CCO", "invalid", "c1ccccc1"])
# Returns None for failed molecules

from molfeat.trans.pretrained import PretrainedMolTransformer

# ChemBERTa — RoBERTa trained on 77M PubChem compounds
chemberta = PretrainedMolTransformer("ChemBERTa-77M-MLM", n_jobs=-1)
embeddings = chemberta(["CCO", "CC(=O)O", "c1ccccc1"])
print(f"ChemBERTa shape: {embeddings.shape}")  # (3, 768)

# GIN — graph neural network pretrained on ChEMBL
gin = PretrainedMolTransformer("gin-supervised-masking", n_jobs=-1)
graph_emb = gin(["CCO", "CC(=O)O"])
print(f"GIN shape: {graph_emb.shape}")

from molfeat.store.modelstore import ModelStore

store = ModelStore()
print(f"Total available: {len(store.available_models)}")

# Search for specific model
results = store.search(name="ChemBERTa")
for model in results:
    print(f"  {model.name}: {model.description}")

# View usage and load
card = store.search(name="ChemBERTa-77M-MLM")[0]
card.usage()
transformer = store.load("ChemBERTa-77M-MLM")

from molfeat.calc import FPCalculator
from molfeat.trans import MoleculeTransformer
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score

# Featurize molecules
transformer = MoleculeTransformer(FPCalculator("ecfp", radius=3), n_jobs=-1)
X = transformer(smiles_train)
print(f"Features shape: {X.shape}")

# Train and evaluate
model = RandomForestRegressor(n_estimators=100)
scores = cross_val_score(model, X, y_train, cv=5, scoring='r2')
print(f"R² = {scores.mean():.3f} ± {scores.std():.3f}")

# Save for deployment
transformer.to_state_yaml_file("production_featurizer.yml")

from sklearn.ensemble import RandomForestClassifier

# Step 1: Featurize known actives/inactives
transformer = MoleculeTransformer(FPCalculator("ecfp"), n_jobs=-1)
X_train = transformer(train_smiles)

# Step 2: Train classifier
clf = RandomForestClassifier(n_estimators=500, n_jobs=-1)
clf.fit(X_train, train_labels)

# Step 3: Screen library (e.g., 1M compounds)
X_screen = transformer(screening_smiles)
predictions = clf.predict_proba(X_screen)[:, 1]

# Step 4: Rank and select top hits
top_indices = predictions.argsort()[::-1][:1000]
top_hits = [screening_smiles[i] for i in top_indices]
print(f"Top 1000 hits selected from {len(screening_smiles)} compounds")

from molfeat.calc import FPCalculator, RDKitDescriptors2D
from sklearn.metrics import roc_auc_score

featurizers = {
    'ECFP': FPCalculator("ecfp"),
    'MACCS': FPCalculator("maccs"),
    'Descriptors': RDKitDescriptors2D(),
}

for name, calc in featurizers.items():
    transformer = MoleculeTransformer(calc, n_jobs=-1)
    X_train = transformer(smiles_train)
    X_test = transformer(smiles_test)
    clf = RandomForestClassifier(n_estimators=100)
    clf.fit(X_train, y_train)
    auc = roc_auc_score(y_test, clf.predict_proba(X_test)[:, 1])
    print(f"{name}: AUC = {auc:.3f}")

from sklearn.pipeline import Pipeline
from sklearn.ensemble import RandomForestClassifier

pipeline = Pipeline([
    ('featurizer', MoleculeTransformer(FPCalculator("ecfp"), n_jobs=-1)),
    ('classifier', RandomForestClassifier(n_estimators=100))
])
pipeline.fit(smiles_train, y_train)
predictions = pipeline.predict(smiles_test)

from sklearn.metrics.pairwise import cosine_similarity

calc = FPCalculator("ecfp")
query_fp = calc("CC(=O)Oc1ccccc1C(=O)O").reshape(1, -1)  # Aspirin

transformer = MoleculeTransformer(calc, n_jobs=-1)
db_fps = transformer(database_smiles)

similarities = cosine_similarity(query_fp, db_fps)[0]
top_k = similarities.argsort()[-10:][::-1]
for i in top_k:
    print(f"  {database_smiles[i]}: {similarities[i]:.3f}")

import numpy as np

def featurize_chunks(smiles_list, transformer, chunk_size=10000):
    all_features = []
    for i in range(0, len(smiles_list), chunk_size):
        chunk = smiles_list[i:i+chunk_size]
        features = transformer(chunk)
        all_features.append(features)
        print(f"Processed {min(i+chunk_size, len(smiles_list))}/{len(smiles_list)}")
    return np.vstack(all_features)