В этой статье показано, как создать сверточную нейронную сеть, способную различать типы рака на основе простого классификатора PyTorch. Данные и код для обучения доступны публично, а процесс обучения можно выполнить на обычном персональном компьютере, в том числе на процессоре.
Рак возникает как побочный эффект накопления ошибок в информации клеток на протяжении жизни, что приводит к неконтролируемому росту. Ученые изучают закономерности этих ошибок, чтобы глубже понять природу заболевания. С точки зрения специалиста по данным, геном человека представляет собой строку длиной около трех миллиардов букв A, C, G, T (то есть 2 бита информации на букву). Ошибка копирования или внешнее воздействие может удалить, вставить или изменить букву, вызвав мутацию и возможное нарушение функции генома.
Тем не менее, отдельные ошибки редко приводят к развитию рака. Организм обладает несколькими механизмами защиты от рака, включая специальные белки — так называемые супрессоры опухолей. Для того чтобы клетка могла обеспечить устойчивый рост, необходимо выполнение ряда условий — так называемых «признаков рака».
Поэтому изменения в отдельных буквах ДНК обычно недостаточны для вызова самоподдерживающегося пролиферативного роста. Большинство раков, опосредованных мутациями (в отличие от других причин, например вируса HPV), также демонстрируют изменения в числе копий (CN). Это крупномасштабные события, часто добавляющие или удаляющие миллионы баз ДНК за раз.
Эти обширные изменения в структуре генома приводят к потере генов, которые препятствовали бы формированию рака, одновременно накапливая гены, стимулирующие рост клеток. Секвенируя ДНК этих клеток, мы можем выявить такие изменения, которые часто происходят в регионах, специфичных для типа рака. Значения числа копий для каждого аллеля можно получить из данных секвенирования с помощью инструментов для определения числа копий.
Обработка профилей числа копий
Одно из преимуществ работы с профилями числа копий (CN) заключается в том, что они не являются биометрическими данными и поэтому могут публиковаться без ограничений доступа. Это позволяет накапливать данные со временем из различных исследований для создания наборов данных достаточного объема. Однако данные из разных исследований не всегда напрямую сопоставимы, поскольку они могут быть получены с использованием различных технологий, иметь разное разрешение или быть предварительно обработанными по-разному.
Для получения данных, их совместной обработки и визуализации мы будем использовать инструмент CNSistent, разработанный в рамках работ Института вычислительной биологии рака при Университетской клинике Кёльна, Германия.
Сначала клонируем репозиторий и данные, переключаясь на версию, использованную в этой статье:
git clone [email protected]:schwarzlab/cnsistent.git cd cnsistent git checkout v0.9.0Поскольку данные находятся внутри репозитория (около 1 ГБ), загрузка занимает несколько минут. Для клонирования требуются Git и Git LFS.
В репозитории есть файл requirements.txt, перечисляющий все зависимости, которые можно установить с помощью pip install -r requirements.txt.
(Рекомендуется сначала создать виртуальную среду). После установки зависимостей CNSistent устанавливается запуском pip install -e . в той же папке. Флаг -e устанавливает пакет из исходной директории, что необходимо для доступа к данным через API.
Репозиторий содержит сырые данные из трех наборов: TCGA, PCAWG и TRACERx. Их нужно предварительно обработать, запустив скрипт bash ./scripts/data_process.sh.
Теперь у нас есть обработанные наборы данных, которые можно загрузить с помощью библиотеки утилит CNSistent:
import cns.data_utils as cdu samples_df, cns_df = cdu.main_load("imp") print(cns_df.head())Это дает следующий результат:
| | 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 |Эта таблица отображает данные о числе копий с колонками:
sample_id: идентификатор образца,chrom: хромосома,start: начальная позиция сегмента (0-индексированная, включительная),end: конечная позиция сегмента (0-индексированная, исключительная),major_cn: число копий основного аллеля (большего из двух),minor_cn: число копий минорного аллеля (меньшего из двух).
В первой строке видно, что образец SP101724 имеет 2 копии основного аллеля и 2 копии минорного (всего 4) в регионе хромосомы 1 от 0 до 27,26 мегабаз.
Второй датафрейм, samples_df, содержит метаданные образцов. Для наших целей важен только тип. Доступные типы можно исследовать так:
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)В приведенном примере видно потенциальную проблему с данными — длины отдельных сегментов неравномерны. Первый сегмент составляет 27,26 мегабаз, а второй — всего 0,77 мегабазы. Это проблема для нейронной сети, которая ожидает входные данные фиксированного размера.
Можно было бы взять все существующие точки разрыва и создать сегменты между всеми точками разрыва в наборе данных, так называемую минимально согласованную сегментацию. Однако это привело бы к огромному числу сегментов — быстрая проверка с помощью len(cns_df["end"].unique()) показывает 823652 уникальных точек разрыва.
В качестве альтернативы CNSistent можно использовать для создания новой сегментации с помощью алгоритма биннинга. Это создаст сегменты фиксированного размера, подходящие для ввода в нейронную сеть. В наших исследованиях сегменты 1–3 мегабаз обеспечивают лучший баланс между точностью и переобучением. Сначала создаем сегментацию, затем применяем ее для получения новых файлов CNS с помощью следующего Bash-скрипта:
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Цикл обрабатывает каждый набор данных отдельно, сохраняя одну и ту же сегментацию. Флаг --threads ускоряет процесс за счет параллельного выполнения агрегации, подстраивая значение под количество доступных ядер.
Аргументы --remove gaps --filter 300000 удаляют регионы низкой маппабельности (так называемые гэпы) и отфильтровывают сегменты короче 300 Кб. Аргумент --split 3000000 создает сегменты по 3 Мб.
Недифференцированный рак лёгких
В этой статье мы сосредоточимся на классификации недифференцированного рака лёгких, который составляет около 85% всех случаев рака лёгких, в частности на различии между аденокарциномой и плоскоклеточным раком. Важно различать эти два типа, поскольку их режимы лечения отличаются, а новые методы дают надежду на неинвазивное выявление из образцов крови или мазков из носа.
Мы используем сегменты, созданные выше, и загружаем их с помощью предоставленной утилиты. Поскольку мы классифицируем между двумя типами рака, фильтруем образцы, чтобы включить только релевантные типы: LUAD (аденокарцинома) и LUSC (плоскоклеточный рак), и строим график первого образца:
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))Сегменты основного и минорного числа копий в бинах по 3 Мб для первых трех образцов. В этом случае все три образца получены из мультирегионного секвенирования одного пациента, демонстрируя, насколько гетерогенны могут быть раковые клетки даже внутри одной опухоли.
Модель сверточной нейронной сети
Для запуска кода требуется Python 3 с установленным PyTorch 2+ и оболочка, совместимая с Bash. Рекомендуется NVIDIA GPU для ускорения обучения, но это не обязательно.
Сначала определяем сверточную нейронную сеть с тремя слоями:
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Это стандартная глубокая сверточная сеть с 2 входными каналами — по одному на каждый аллель — и 3 сверточными слоями с 1D-ядром размером 3 и функцией активации ReLU. За сверточными слоями следуют слои максимального пулинга с ядром размером 2. Свертка традиционно применяется для обнаружения краев, что полезно для нас, поскольку мы интересуемся изменениями в числе копий, то есть краями сегментов.
Выход сверточных слоев уплощается и передается через два полносвязных слоя с дропаутом. Слой LazyLinear соединяет выход 64 накопленных каналов в один слой из 128 узлов, без необходимости рассчитывать количество узлов в конце свертки. Здесь сосредоточено большинство параметров, поэтому мы применяем дропаут для предотвращения переобучения.
Обучение модели
Сначала преобразуем датафреймы в тензоры Torch. Используем утилиту bins_to_features, которая создает 3D-массив признаков в формате (образцы, аллели, сегменты). При этом данные делятся на обучающую и тестовую выборки в соотношении 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)Теперь можно обучать модель с помощью следующего цикла обучения на 20 эпох. Оптимизатор Adam и функция потерь CrossEntropy обычно используются для задач классификации, поэтому мы применяем их здесь:
# 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}') Это завершает обучение. Затем оцениваем модель и выводим матрицу путаницы:
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')Процесс обучения занимает около 7 секунд на GPU NVIDIA RTX 4090.
Заключение
Мы разработали эффективный и точный классификатор подтипа рака лёгких на основе данных о числе копий. Как показано, такие модели хорошо переносятся на новые исследования и источники последовательностей данных.
Масштабный ИИ часто оправдывается, в том числе как «решение проблемы рака». Однако, как в этой статье, небольшие модели с классическими подходами обычно справляются с задачей отлично. Некоторые утверждают, что реальное препятствие для машинного обучения в биологии и медицине — не в решении проблем, а в создании реального воздействия на пациентов.
Тем не менее, машинное обучение уже в основном решило одну крупную загадку вычислительной биологии, возродив интерес к машинному обучению для рака. С некоторым везением, следующее десятилетие может стать периодом, когда мы наконец «решим» рак.
Бонус: Cell2Sentence
Современные базовые модели часто содержат информацию, такую как гены с высоким числом копий в различных раках, как часть их обучающего набора. Это привело к созданию подходов на основе LLM, таких как Cell2Sentence, где набор данных преобразуется в естественный язык и подается в LLM.
Мы можем применить этот подход здесь. Сначала преобразуем случайный образец LUSC в ранжированный список генов на основе их значения числа копий:
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()Затем преобразуем список генов в запрос на естественном языке:
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)Поданный в ChatGPT 5, модель правильно идентифицирует рак 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).Однако генерация такого вывода для одного образца занимает больше времени, чем классификация всего набора нашей моделью, и стоила бы около 200 долларов на API-расходах для всего нашего набора данных.