Comment classer le sous-type de cancer du poumon à partir des numéros de copies d’ADN à l’aide de PyTorch

Je vais démontrer comment construire un réseau neuronal à convolution capable de distinguer les types de cancer à l’aide d’un simple classificateur PyTorch. Les données et le code utilisés pour la formation sont accessibles au public et la formation peut être effectuée sur un ordinateur personnel, potentiellement même sur un processeur.

Le cancer est un effet secondaire malheureux de l’accumulation d’erreurs d’information par nos cellules au cours de notre vie, conduisant à une croissance incontrôlée. Dans le cadre de nos recherches, nous étudions les schémas de ces erreurs afin de mieux comprendre la maladie. Vu du point de vue d’un data scientist, le génome humain est une chaîne d’environ trois milliards de lettres avec les lettres A, C, G, T (soit 2 bits d’information par lettre). Une erreur de copie ou un événement externe peut potentiellement supprimer/insérer/modifier une lettre, provoquant une mutation et potentiellement une perturbation de la fonction génomique.

Cependant, les erreurs individuelles ne conduisent presque jamais au développement d’un cancer. Le corps humain dispose de multiples mécanismes pour prévenir le développement du cancer, notamment des protéines dédiées, appelées suppresseurs de tumeurs. Une liste de conditions nécessaires, appelées « les caractéristiques du cancer » doivent être remplies pour qu’une cellule puisse créer une croissance soutenue.

Par conséquent, les modifications apportées aux lettres individuelles de l’ADN sont généralement insuffisantes pour provoquer une croissance proliférative auto-entretenue. La grande majorité des cancers à médiation mutationnelle (par opposition à d’autres sources de cancer, par exemple le virus HPV) présentent également des modifications du numéro de copie (CN). Il s’agit d’événements à grande échelle, qui ajoutent ou suppriment souvent des millions de bases d’ADN à la fois.

Ces vastes changements dans la structure du génome entraînent la perte de gènes qui empêcheraient la formation du cancer, tout en accumulant des gènes favorisant la croissance cellulaire. En séquençant l’ADN de ces cellules, nous pouvons identifier ces changements, qui se produisent assez souvent dans des régions spécifiques au type de cancer. Les valeurs du nombre de copies pour chaque allèle peuvent être dérivées des données de séquençage en utilisant copier le numéro des appelants.

Traitement des profils de numéros de copie

L’un des avantages de travailler avec des profils Copy Number (CN) est qu’ils ne sont pas biométriques et peuvent donc être publiés sans qu’il soit nécessaire de restreindre l’accès. Cela nous permet d’accumuler des données au fil du temps provenant de plusieurs études pour créer des ensembles de données de taille suffisante. Cependant, les données provenant de différentes études ne sont pas toujours directement comparables, car elles peuvent être générées à l’aide de différentes technologies, avoir des résolutions différentes ou être prétraitées de différentes manières.

Pour obtenir les données et les traiter et visualiser conjointement, nous utiliserons l’outil CNSistantélaboré dans le cadre des travaux du Institut de biologie computationnelle du cancer de la Clinique universitaire, Cologne, Allemagne.

Nous clonons d’abord le référentiel et les données et définissons la version utilisée dans ce texte :

git clone [email protected]:schwarzlab/cnsistent.git
cd cnsistent
git checkout v0.9.0

Étant donné que les données que nous utiliserons se trouvent dans le référentiel (~ 1 Go de données), le téléchargement prend quelques minutes. Pour le clonage, Git et Git LFS doivent être présents sur le système.

À l’intérieur du référentiel se trouve un fichier exigences.txt qui répertorie toutes les dépendances pouvant être installées à l’aide de pip install -r requirements.txt.

(Il est recommandé de créer d’abord un environnement virtuel). Une fois la configuration requise installée, CNSistent peut être installé en exécutant pip install -e . dans le même dossier. Le -e flag installe le package à partir de son répertoire source, ce qui est nécessaire pour accéder aux données via l’API.

Le référentiel contient des données brutes provenant de trois ensembles de données : TCGA, PCAWGet TRACEURx. Ceux-ci doivent d’abord être prétraités. Cela peut être fait en exécutant le script bash ./scripts/data_process.sh.

Maintenant, nous avons traité les ensembles de données et pouvons les charger à l’aide de la bibliothèque d’utilitaires de données CNSitent :

import cns.data_utils as cdu
samples_df, cns_df = cdu.main_load("imp")
print(cns_df.head())

Produire le résultat suivant :

|    | sample_id   | chrom   |    start |      end |   major_cn |   minor_cn |
|---:|:------------|:--------|---------:|---------:|-----------:|-----------:|
|  0 | SP101724    | chr1    |        0 | 27256755 |          2 |          2 |
|  1 | SP101724    | chr1    | 27256755 | 28028200 |          3 |          2 |
|  2 | SP101724    | chr1    | 28028200 | 32976095 |          2 |          2 |
|  3 | SP101724    | chr1    | 32976095 | 33354394 |          5 |          2 |
|  4 | SP101724    | chr1    | 33354394 | 33554783 |          3 |          2 |

Ce tableau montre les données du numéro de copie avec les colonnes suivantes :

• sample_id: l’identifiant de l’échantillon,
• chrom: le chromosome,
• start: la position de départ du segment (indexé 0 inclus),
• end: la position finale du segment (exclusif indexé 0),
• major_cn: le nombre de copies de l’allèle majeur (le plus gros des deux),
• minor_cn: le nombre de copies de l’allèle mineur (le plus petit des deux).

Sur la première ligne on peut donc voir un segment indiquant que l’échantillon SP101724 possède 2 copies de l’allèle majeur et 2 copies de l’allèle mineur (4 au total) dans la région du chromosome 1 de 0 à 27,26 mégabase.

La deuxième trame de données que nous avons chargée, samples_df, contient les métadonnées des échantillons. Pour nos besoins, seul le type est important. Nous pouvons étudier les types disponibles en exécutant :

import matplotlib.pyplot as plt
type_counts = samples_df["type"].value_counts()
plt.figure(figsize=(10, 6))
type_counts.plot(kind='bar')
plt.ylabel('Count')
plt.xticks(rotation=90)

Dans l’exemple ci-dessus, nous pouvons observer un problème potentiel avec les données : les longueurs des segments individuels ne sont pas uniformes. Le premier segment mesure 27,26 mégabases, tandis que le second ne mesure que 0,77 mégabases. C’est un problème pour le réseau neuronal, qui s’attend à ce que l’entrée soit d’une taille fixe.

Nous pourrions techniquement prendre tous les points d’arrêt existants et créer des segments entre tous les points d’arrêt de l’ensemble de données, ce qu’on appelle segmentation cohérente minimale. Cela entraînerait cependant un très grand nombre de segments — une vérification rapide à l’aide de len(cns_df[“end”].unique()) montre qu’il y a 823652 points d’arrêt uniques.

Alternativement, nous pouvons utiliser CNSistent pour créer une nouvelle segmentation à l’aide d’un algorithme de binning. Cela créera des segments de taille fixe, qui pourront être utilisés comme entrée du réseau neuronal. Dans notre travail, nous avons déterminé 1 à 3 segments de mégabase pour offrir le meilleur compromis entre précision et surajustement. Nous créons d’abord la segmentation puis l’appliquons pour obtenir de nouveaux fichiers CNS à l’aide du script Bash suivant :

threads=8
cns segment whole --out "./out/segs_3MB.bed" --split 3000000 --remove gaps - filter 300000 
for dataset in TRACERx PCAWG TCGA_hg19; 
do 
 cns aggregate ./out/${dataset}_cns_imp.tsv - segments ./out/segs_3MB.bed - out ./out/${dataset}_bin_3MB.tsv - samples ./out/${dataset}_samples.tsv - threads $threads
done

La boucle traite chaque ensemble de données séparément, tout en conservant la même segmentation. Le--threadsL’indicateur est utilisé pour accélérer le processus en exécutant l’agrégation en parallèle, en ajustant la valeur en fonction du nombre de cœurs disponibles.

Le--remove gaps  --filter 300000les arguments supprimeront les régions de faible mappabilité (c’est-à-dire lacunes) et filtrez les segments de moins de 300 Ko. Le --split 3000000 L’argument créera des segments de 3 Mo.

Carcinome pulmonaire non à petites cellules

Dans ce texte, nous nous concentrerons sur la classification des carcinomes pulmonaires non à petites cellules, qui représentent environ 85 % de tous les cancers du poumon, en particulier la distinction entre adénocarcinome et carcinome épidermoïde. Il est important de faire la différence entre les deux car leurs régimes de traitement seront différents et de nouvelles méthodes laissent espérer une détection non invasive à partir d’échantillons de sang ou d’écouvillons nasaux.

Nous utiliserons les segments produits ci-dessus et les chargerons à l’aide d’une fonction utilitaire fournie à l’aide d’une fonction utilitaire. Puisque nous classons deux types de cancer, nous pouvons filtrer les échantillons pour inclure uniquement les types pertinents, LUAD (adénocarcinome) et LUSC (carcinome épidermoïde) et tracez le premier échantillon :

import cns
samples_df, cns_df = cdu.main_load("3MB")
samples_df = samples_df.query("type in ['LUAD', 'LUSC']")
cns_df = cns.select_CNS_samples(cns_df, samples_df)
cns_df = cns.only_aut(cns_df)
cns.fig_lines(cns.cns_head(cns_df, n=3))

Segments de numéros de copie majeurs et mineurs dans des bacs de 3 Mo pour les trois premiers échantillons. Dans ce cas, les trois échantillons proviennent d’un séquençage multirégional du même patient, démontrant à quel point les résultats sont hétérogènes. les cellules cancéreuses peuvent même se trouver dans une seule tumeur.

Modèle de réseau neuronal à convolution

L’exécution du code nécessite Python 3 avec PyTorch 2+ à installer et un shell compatible Bash. Le GPU NVIDIA est recommandé pour un entraînement plus rapide, mais pas nécessaire.

Nous définissons d’abord un réseau de neurones convolutifs à trois couches :

import torch.nn as nn

class CNSConvNet(nn.Module):
    def __init__(self, num_classes):
        super(CNSConvNet, self).__init__()

        self.conv_layers = nn.Sequential(
            nn.Conv1d(in_channels=2, out_channels=16, kernel_size=3, padding=1),
            nn.ReLU(),
            nn.MaxPool1d(kernel_size=2),
            nn.Conv1d(in_channels=16, out_channels=32, kernel_size=3, padding=1),
            nn.ReLU(),
            nn.MaxPool1d(kernel_size=2),
            nn.Conv1d(in_channels=32, out_channels=64, kernel_size=3, padding=1),
            nn.ReLU(),
            nn.MaxPool1d(kernel_size=2)
        )

        self.fc_layers = nn.Sequential(
            nn.LazyLinear(128),
            nn.ReLU(),
            nn.Dropout(0.5),
            nn.Linear(128, num_classes)
        )

    def forward(self, x):
        x = self.conv_layers(x)
        x = x.view(x.size(0), -1) 
        x = self.fc_layers(x)
        return x

Il s’agit d’un CNN profond passe-partout avec 2 canaux d’entrée – un pour chaque allèle – et 3 couches convolutives utilisant un noyau 1D de taille 3 et la fonction d’activation ReLU. Les couches convolutives sont suivies par des couches de pooling maximales avec une taille de noyau de 2. La convolution est traditionnellement utilisée pour détection de bordce qui nous est utile car nous nous intéressons aux changements dans le nombre de copies, c’est-à-dire les bords des segments.

La sortie des couches convolutives est ensuite aplatie et passée à travers deux couches entièrement connectées avec décrochage. Le LazyLinearLa couche connecte la sortie de 64 canaux empilés en une couche de 128 nœuds, sans avoir besoin de calculer le nombre de nœuds à la fin de la convolution. C’est là que se trouvent la plupart de nos paramètres, c’est pourquoi nous appliquons également un abandon pour éviter le surajustement.

Entraîner le modèle

Nous devons d’abord convertir les dataframes en tenseurs Torch. Nous utilisons une fonction utilitaire bins_to_features, qui crée un tableau de fonctionnalités 3D du format (échantillons, allèles, segments). Au cours du processus, nous avons également divisé les données en ensembles de formation et de test dans un rapport de 4 : 1 :

import torch
from torch.utils.data import TensorDataset, DataLoader
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split

# convert data to features and labels
features, samples_list, columns_df = cns.bins_to_features(cns_df)

# convert data to Torch tensors
X = torch.FloatTensor(features)
label_encoder = LabelEncoder()
y = torch.LongTensor(label_encoder.fit_transform(samples_df.loc[samples_list]["type"])) 

# Test/train split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

# Create dataloaders
train_loader = DataLoader(TensorDataset(X_train, y_train), batch_size=32, shuffle=True)
test_loader = DataLoader(TensorDataset(X_test, y_test), batch_size=32, shuffle=False)

Nous pouvons maintenant entraîner le modèle en utilisant la boucle d’entraînement suivante avec 20 époques. L’optimiseur Adam et la perte CrossEntropy sont généralement utilisés pour les tâches de classification, nous les utilisons donc également ici :

# setup the model, loss, and optimizer
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
model = CNSConvNet(num_classes=len(label_encoder.classes_)).to(device)
criterion = nn.CrossEntropyLoss()
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)

# Training loop
num_epochs = 20
for epoch in range(num_epochs):
    model.train()
    running_loss = 0.0
    
    for inputs, labels in train_loader:
        inputs, labels = inputs.to(device), labels.to(device)
        
        # Clear gradients
        optimizer.zero_grad()
        
        # Forward pass
        outputs = model(inputs)
        loss = criterion(outputs, labels)
        
        # Backward pass and optimize
        loss.backward()
        optimizer.step()
        
        running_loss += loss.item()
    
    # Print statistics
    print(f'Epoch {epoch+1}/{num_epochs}, Loss: {running_loss/len(train_loader):.4f}')

Ceci conclut la formation. Ensuite, nous pouvons évaluer le modèle et imprimer la matrice de confusion :

import numpy as np
from sklearn.metrics import confusion_matrix
import seaborn as sns

# Loop over batches in the test set and collect predictions
model.eval()
y_true = []
y_pred = []
with torch.no_grad():
    for inputs, labels in test_loader:
        inputs, labels = inputs.to(device), labels.to(device)
        outputs = model(inputs)
        y_true.extend(labels.cpu().numpy())
        y_pred.extend(outputs.argmax(dim=1).cpu().numpy())
        _, predicted = torch.max(outputs.data, 1)

# Calculate accuracy and confusion matrix
accuracy = (np.array(y_true) == np.array(y_pred)).mean()
cm = confusion_matrix(y_true, y_pred)

# Plot the confusion matrix
plt.figure(figsize=(3, 3), dpi=200)
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', xticklabels=label_encoder.classes_, yticklabels=label_encoder.classes_)
plt.xlabel('Predicted')
plt.ylabel('True')
plt.title('Confusion Matrix, accuracy={:.2f}'.format(accuracy))
plt.savefig("confusion_matrix.png", bbox_inches='tight')

Le processus de formation prend environ 7 secondes au total sur un GPU NVIDIA RTX 4090.

Conclusion

Nous avons développé un classificateur efficace et précis du sous-type de cancer du poumon à partir des données sur le nombre de copies. Comme nous l’avons montréces modèles se transfèrent bien aux nouvelles études et sources de données de séquence.

L’IA à grande échelle est souvent justifiée, entre autres, comme un «solution au cancer« . Cependant, comme dans cet article, les petits modèles avec des approches classiques remplissent généralement bien leur objectif. Certains prétendent même que le véritable obstacle de l’apprentissage automatique en biologie et en médecine n’est pas la résolution de problèmes, mais l’impact réel sur les patients.

Pourtant, l’apprentissage automatique a permis de résoudre au moins une énigme majeure de la biologie computationnelle, revigorer l’accent mis sur l’apprentissage automatique pour le cancer. Avec un peu de chance, nous pourrons peut-être envisager la prochaine décennie comme celle où nous aurons enfin « résolu » le cancer.

Bonus : Cell2Sentence

Les modèles fondamentaux d’aujourd’hui contiennent souvent des informations, telles que les gènes qui ont un nombre élevé de copies dans quels cancers, dans le cadre de leur ensemble de formation. Cela a cependant conduit à la création d’approches basées sur le LLM, comme Cell2Phraseoù un ensemble de données est converti en langage naturel et soumis à un LLM.

Nous pouvons utiliser cette approche ici. Nous convertissons d’abord un échantillon aléatoire de LUSC en une liste classée de gènes, en fonction de leur valeur de nombre de copies :

import cns
import cns.data_utils as cdu
samples_df, cns_df = cdu.main_load("COSMIC")
cns_df = cns.add_total_cn(cns_df)
lusc_df = cns.select_cns_by_type(cns_df, samples_df, "LUSC")
lusc_one_df = cns.cns_head(lusc_df, 1)
rank_list = lusc_one_df.sort_values(by="total_cn", ascending=False)["name"].to_list()

Ensuite, nous convertissons la liste des gènes en une requête en langage naturel :

query = "Below is a list of COSMIC {gene_names} ordered by descending total copy number in a cell. Your task is to identify the cancer type of that sample.\n\n"
query += "{gene_names}: " + " ".join(rank_list) + "\n\n"
query += "Answer:"
print(query)

Soumis à ChatGPT 5, le modèle identifie correctement le cancer LUSC :

Based on the gene copy-number profile you listed, the strongest clue is the massive amplification of 11q13 (CCND1, FADD) together with SOX2, PIK3CA, TP63, EGFR, FGFR1, MYC, CCND1/CCND2/CCND3, CDK6, CDK4. This is a very characteristic pattern of squamous cell carcinoma of the lung (LUSC):
CCND1 amplification is one of the most frequent events in LUSC.
SOX2 amplification (3q26) is a defining driver of lung squamous carcinoma, rarely seen in adenocarcinoma.
TP63 amplification is again highly specific to squamous tumors.
Co-amplifications of PIK3CA, FGFR1, EGFR, MYC fit this profile.
The broad set of genes amplified (many regulators of squamous differentiation and RTK/RAS/PI3K signaling) reflects typical copy-number–driven LUSC biology.
By contrast, lung adenocarcinomas usually show KRAS, EGFR, ALK, ROS1, MET exon 14, ERBB2 point mutations/fusions and have fewer widespread squamous-type CNAs.
Answer: Lung squamous cell carcinoma (LUSC).

Cependant, générer cette sortie pour un échantillon prend plus de temps que classer l’ensemble de nos données avec notre modèle et cela coûterait environ 200 $ en frais d’API pour classer l’ensemble de notre ensemble de données.