Видимо, астрологи Хабра объявили день Python в биологии. Сегодня мы приготовили для вас материал, в котором аспирант Эрин Уилсон делится кодом умозрительного примера прогнозирования количественной оценки последовательностей ДНК. За подробностями приглашаем под кат — к старту нашего флагманского курса по Data Science.

ДНК — поток комплексных данных, который можно представить в виде нуклеотидной последовательности ACTG. Тем не менее многие сложные паттерны и структурные нюансы трудно понять человеку, который смотрит на необработанную последовательность нуклеотидов.

В последние годы в моделировании ДНК достигли заметного прогресса. Этот прогресс стал возможным благодаря технологиям глубокого обучения.

Для моделирования ДНК исследователи применяли различные методы. В их числе свёрточные нейросети (CNN), сети с долговременной и кратковременной памятью (LSTM) и даже трансформеры для прогноза различных геномных измерений непосредственно по последовательностям ДНК. Самое полезное в этих моделях — автоматический поиск паттернов и мотивов последовательностей ДНК при достаточном количестве качественных обучающих данных. При этом исследователь не должен заранее указывать искомые паттерны. Нейросеть может их прогнозировать. Полезность глубокого обучения в геномике прибавляет энтузиазма всем, кому нужно сопоставлять последовательности ДНК с биологическими функциями!

Как аспирант, увлечённая математическими методами решения задач рационального природопользования и синтетической биологии, я с интересом изучила возможности PyTorch для исследования паттернов последовательностей ДНК.

Обучающих материалов по PyTorch немало, но многие из них посвящены работе с графикой и языковыми данными. ДНК использовалась в качестве входных данных в целом ряде интересных проектов, которые разработали в PyTorch фреймворки для моделирования различных биологических явлений 1, 2, 3, но они могут быть довольно трудными для новичка.

Мне было нелегко найти примеры в PyTorch для новичков, которые также были бы связаны с ДНК, поэтому я составила это краткое руководство на случай, если будущие архитекторы моделей ДНК найдут его полезным для начала работы.

Само руководство можно запустить в интерактивном режиме блокнота Jupyter. Проследить его план можно по краткому описанию ключевых понятий и Github gists в остальной части этой статьи.

Научим решать интересные задачи, чтобы вы прокачали карьеру или стали востребованным IT-специалистом:

• Полный курс по Data Science. Получите одну из самых перспективных профессий за 24 месяца.
• «Fullstack-разработчик на Python». Станьте незаменимым специалистом-универсалом за 15 месяцев.
• «Уверенный старт в IT». Пройдите лучший курс для новичков: попробуйте 9 профессий и освойте подходящую именно вам.

Построим модель для прогноза количественных показателей по последовательности ДНК в PyTorch

В этом руководстве рассказывается, как PyTorch использует необработанные последовательности ДНК в качестве входных данных и вводит их в модель нейросети, чтобы прогнозировать количественный балл по самой последовательности.

План руководства

1. Синтез данных ДНК.
2. Подготовка данных для обучения в PyTorch.
3. Определение модели в PyTorch.
4. Определение функций обучающего цикла.
5. Запуск моделей.
6. Проверка прогнозов модели на тестовом наборе данных.
7. Визуализация свёрточных фильтров.
8. Итоги.

Полагаю, что читатель знаком со следующими аспектами машинного обучения:

• Нейросеть и свёрточная нейросеть (CNN).
• Эпохи обучения моделей.
• Разделение данных на наборы train/val/test.
• Функции потерь и сравнение кривых потерь train и val.

Также полагаю, что читатель знаком с некоторыми понятиями биологии. Это:

• Нуклеотиды в ДНК.
• Регуляторный мотив.
• Визуализация мотивов ДНК.

Я не утверждаю, что эти методы оптимальны. Я всё ещё учусь, и это лишь моё представление о решении задачи. Уверена, что есть решения элегантнее.

1. Синтез данных ДНК

Обычно учёные хотят прогнозировать балл связывания, силу экспрессии генов, классифицировать события связывания фактора транскрипции и тому подобное. У нас всё проще. В этом руководстве мы пытаемся научить модель обнаруживать в последовательности ДНК очень маленький и простой паттерн и присваивать ему соответствующий балл. Для этого мы используем глубокое обучение. Но практическая задача даже не в этом. Хочется просто убедиться, что правильно настроили элементы PyTorch и модель может учиться на входных данных, похожих на последовательность ДНК).

Допустим, у нас есть октамерная последовательность ДНК. Начислим ей определённое число баллов за каждый нуклеотид:

• A (аденин) = +20 баллов.
• C (цитозин) = +17 баллов.
• G (гуанин) = +14 баллов.
• T (тимин) = +11 баллов.

Для каждого октамера суммируем балл и усредняем его. Например:

AAAAAAAA даёт нам 20,0

mean(20 + 20 + 20 + 20 + 20 + 20 + 20 + 20) = 20,0

ACAAAAAA даёт 19,625

mean(20 + 17 + 20 + 20 + 20 + 20 + 20 + 20) = 19,625

Баллы выбраны для нуклеотидов произвольно — никакой настоящей биологии в этом нет! Мы просто назначаем последовательностям баллы, чтобы попрактиковаться в работе с PyTorch.

Во многих работах свёрточные нейросети (CNN) определяют «мотивы» или короткие паттерны ДНК. Они могут активировать или подавлять биологические реакции. Поэтому добавим в систему начисления баллов ещё один элемент: иногда нам нужно моделировать мотивы, влияющие на экспрессию генов. Допустим, последовательность получает +10 баллов при наличии мотива TAT в любой части октамера или -10 баллов при наличии мотива GCG. Повторюсь, к реальности это никакого отношения не имеет. Я просто показала вам на примере, как работает простой механизм моделирования активации и подавления.

Простая система подсчёта баллов для октамерных последовательностей ДНК

Вот реализация этой простой системы:

def kmers(k):
    '''Generate a list of all k-mers for a given k'''
    
    return [''.join(x) for x in product(['A','C','G','T'], repeat=k)]
    
seqs8 = kmers(8)
print('Total 8mers:',len(seqs8))
# prints: Total 8mers: 65536

score_dict = {
    'A':20,
    'C':17,
    'G':14,
    'T':11
}
                  
def score_seqs_motif(seqs):
    '''
    Calculate the scores for a list of sequences based on 
    the above score_dict
    '''
    data = []
    for seq in seqs:
        # get the average score by nucleotide
        score = np.mean([score_dict[base] for base in seq])
        
        # give a + or - bump if this k-mer has a specific motif
        if 'TAT' in seq:
            score += 10
        if 'GCG' in seq:
            score -= 10
        data.append([seq,score])
        
    df = pd.DataFrame(data, columns=['seq','score'])
    return df
    
mer8 = score_seqs_motif(seqs8)

Строим распределение баллов для октамерных последовательностей и видим, что они разделились на 3 группы:

• последовательности с мотивом GCG (~5 баллов);
• последовательности без мотивов (~15 баллов);
• последовательности с мотивом TAT (~25 баллов).

Распределение баллов для октамерных последовательностей

Теперь нам предстоит научить модель прогнозировать балл по последовательности ДНК.

2. Подготовка данных для обучения в PyTorch

Чтобы научить нейросеть прогнозировать балл, нужно дать ей входные данные в виде матрицы чисел. Поясню на примере. Пусть мы хотим разделить фотографии с кошками и без кошек на два набора. Сеть распознаёт изображение как матрицу значений пикселей и изучает паттерны, связанные с относительным расположением пикселей (например, паттерны, соответствующие кошачьим ушам или носу с усами).

Аналогичным образом в матрицу чисел нужно превратить последовательности ДНК (строки ACGT). Но как же нам сделать из нашей ДНК «кошку»?

Одна из распространённых стратегий — прямое кодирование [one-hot encoding] ДНК с одним активным состоянием. Для такого кодирования рассмотрим каждый нуклеотид как вектор длиной 4, где в трёх позициях стоят нули, а в одной — единица. Место единицы зависит от типа нуклеотида.

Благодаря прямому кодированию компьютер определит ДНК так же, как кошку на картинке!

def one_hot_encode(seq):
    """
    Given a DNA sequence, return its one-hot encoding
    """
    # Make sure seq has only allowed bases
    allowed = set("ACTGN")
    if not set(seq).issubset(allowed):
        invalid = set(seq) - allowed
        raise ValueError(f"Sequence contains chars not in allowed DNA alphabet (ACGTN): {invalid}")
        
    # Dictionary returning one-hot encoding for each nucleotide 
    nuc_d = {'A':[1.0,0.0,0.0,0.0],
             'C':[0.0,1.0,0.0,0.0],
             'G':[0.0,0.0,1.0,0.0],
             'T':[0.0,0.0,0.0,1.0],
             'N':[0.0,0.0,0.0,0.0]}
    
    # Create array from nucleotide sequence
    vec=np.array([nuc_d[x] for x in seq])
        
    return vec
  
# look at DNA seq of 8 As
a8 = one_hot_encode("AAAAAAAA")
print("AAAAAAAA:\n",a8)

# prints:
# AAAAAAAA:
# [[1. 0. 0. 0.]
# [1. 0. 0. 0.]
# [1. 0. 0. 0.]
# [1. 0. 0. 0.]
# [1. 0. 0. 0.]
# [1. 0. 0. 0.]
# [1. 0. 0. 0.]
# [1. 0. 0. 0.]]

# look at DNA seq of random nucleotides
s = one_hot_encode("AGGTACCT")
print("AGGTACCT:\n",s)
print("shape:",s.shape)

# prints:
# AGGTACCT:
# [[1. 0. 0. 0.]
# [0. 0. 1. 0.]
# [0. 0. 1. 0.]
# [0. 0. 0. 1.]
# [1. 0. 0. 0.]
# [0. 1. 0. 0.]
# [0. 1. 0. 0.]
# [0. 0. 0. 1.]]
# shape: (8, 4)

При помощи прямого кодирования можно подготовить наборы train, val и test. Соответствующая функция quick_split случайным образом выбирает индексы в кадре данных pandas для разделения (функция для этого есть и в sklearn).

В реальных задачах, не связанных с синтезом, может потребоваться более интеллектуальная стратегия разделения. Всё зависит от вашей прогностической задачи. В научных статьях разделение train/test часто производится по хромосоме или по другим особенностям расположения генома.

def quick_split(df, split_frac=0.8):
    '''
    Given a df of samples, randomly split indices between
    train and test at the desired fraction
    '''
    cols = df.columns # original columns, use to clean up reindexed cols
    df = df.reset_index()

    # shuffle indices
    idxs = list(range(df.shape[0]))
    random.shuffle(idxs)

    # split shuffled index list by split_frac
    split = int(len(idxs)*split_frac)
    train_idxs = idxs[:split]
    test_idxs = idxs[split:]
    
    # split dfs and return
    train_df = df[df.index.isin(train_idxs)]
    test_df = df[df.index.isin(test_idxs)]
        
    return train_df[cols], test_df[cols]
    
    
full_train_df, test_df = quick_split(mer8)
train_df, val_df = quick_split(full_train_df)

print("Train:", train_df.shape)
print("Val:", val_df.shape)
print("Test:", test_df.shape)

# prints: 
# Train: (41942, 2)
# Val: (10486, 2)
# Test: (13108, 2)

При подготовке данных в PyTorch важным шагом является применение объектов DataLoader и Dataset. Чтобы разобраться в этом, мне пришлось изрядно погуглить, а чтобы найти решение, пришлось прочесать кипу документов и постов со StackOverflow!

Dataset заключает данные в объект, которому можно легко задать отформатированные примеры X и маркировку Y для обучаемой модели. DataLoader получает Dataset и ряд других деталей о том, как объединять данные в серии и упрощает прохождение итераций обучения.

from torch.utils.data import Dataset, DataLoader

## Here is a custom defined Dataset object specialized for one-hot encoded DNA:

class SeqDatasetOHE(Dataset):
    '''
    Dataset for one-hot-encoded sequences
    '''
    def __init__(self,
                 df,
                 seq_col='seq',
                 target_col='score'
                ):
        # +--------------------+
        # | Get the X examples |
        # +--------------------+
        # extract the DNA from the appropriate column in the df
        self.seqs = list(df[seq_col].values)
        self.seq_len = len(self.seqs[0])
        
        # one-hot encode sequences, then stack in a torch tensor
        self.ohe_seqs = torch.stack([torch.tensor(one_hot_encode(x)) for x in self.seqs])
    
        # +------------------+
        # | Get the Y labels |
        # +------------------+
        self.labels = torch.tensor(list(df[target_col].values)).unsqueeze(1)
        
    def __len__(self): return len(self.seqs)
    
    def __getitem__(self,idx):
        # Given an index, return a tuple of an X with it's associated Y
        # This is called inside DataLoader
        seq = self.ohe_seqs[idx]
        label = self.labels[idx]
        
        return seq, label


## Here is how I constructed DataLoaders from Datasets.

def build_dataloaders(train_df,
                      test_df,
                      seq_col='seq',
                      target_col='score',
                      batch_size=128,
                      shuffle=True
                     ):
    '''
    Given a train and test df with some batch construction
    details, put them into custom SeqDatasetOHE() objects. 
    Give the Datasets to the DataLoaders and return.
    '''
    
    # create Datasets    
    train_ds = SeqDatasetOHE(train_df,seq_col=seq_col,target_col=target_col)
    test_ds = SeqDatasetOHE(test_df,seq_col=seq_col,target_col=target_col)

    # Put DataSets into DataLoaders
    train_dl = DataLoader(train_ds, batch_size=batch_size, shuffle=shuffle)
    test_dl = DataLoader(test_ds, batch_size=batch_size)
    
    return train_dl,test_dl
    
    
train_dl, val_dl = build_dataloaders(train_df, val_df)

Теперь наши объекты DataLoader можно использовать в обучающем цикле!

3. Определение модели PyTorch

Сначала мне хотелось попробовать метод свёрточных нейросетей, потому что в изучении данных генома он показал себя отлично. Для сравнения я взяла простую линейную модель. Вот несколько определений модели:

# very simple linear model
class DNA_Linear(nn.Module):
    def __init__(self, seq_len):
        super().__init__()
        self.seq_len = seq_len
        # the 4 is for our one-hot encoded vector length 4!
        self.lin = nn.Linear(4*seq_len, 1)

    def forward(self, xb):
        # reshape to flatten sequence dimension
        xb = xb.view(xb.shape[0],self.seq_len*4)
        # Linear wraps up the weights/bias dot product operations
        out = self.lin(xb)
        return out

# basic CNN model
class DNA_CNN(nn.Module):
    def __init__(self,
                 seq_len,
                 num_filters=32,
                 kernel_size=3):
        super().__init__()
        self.seq_len = seq_len
        
        self.conv_net = nn.Sequential(
            # 4 is for the 4 nucleotides
            nn.Conv1d(4, num_filters, kernel_size=kernel_size),
            nn.ReLU(inplace=True),
            nn.Flatten(),
            nn.Linear(num_filters*(seq_len-kernel_size+1), 1)
        ) 

    def forward(self, xb):
        # permute to put channel in correct order
        # (batch_size x 4channel x seq_len)
        xb = xb.permute(0,2,1) 
        
        #print(xb.shape)
        out = self.conv_net(xb)
        return out

Это не оптимальные модели, но надо же с чего-то начинать (ещё раз, мы просто учимся применять PyTorch в контексте ДНК).

• Линейная модель пытается прогнозировать балл по весу нуклеотидов на каждой позиции.
• Модель свёрточных нейросетей использует 32 фильтра с длиной (kernel_size) 3 для сканирования всей октамерной последовательности и поиска в ней паттернов-тримеров.

4. Определение функций обучающего цикла

Теперь нужно определить цикл обучения. Скажу честно, я не слишком-то уверена в этом решении. Я долго продиралась через ошибки несоответствия размерности матрицы. Наверное, есть красивые подходы! Но, может быть, и так сойдёт? Напишите мне пару слов, если у вас есть свои мысли на этот счёт.

Так или иначе, я пришла к такому определению стека функций:

# adds default optimizer and loss function
run_model()
    # loops through epochs
    fit()
        # loop through batches
        train_step()
            # calc train loss for batch
            loss_batch()
        val_step()
            # calc val loss for batch
            loss_batch()

# +--------------------------------+
# | Training and fitting functions |
# +--------------------------------+

def loss_batch(model, loss_func, xb, yb, opt=None,verbose=False):
    '''
    Apply loss function to a batch of inputs. If no optimizer
    is provided, skip the back prop step.
    '''
    if verbose:
        print('loss batch ****')
        print("xb shape:",xb.shape)
        print("yb shape:",yb.shape)
        print("yb shape:",yb.squeeze(1).shape)
        #print("yb",yb)

    # get the batch output from the model given your input batch 
    # ** This is the model's prediction for the y labels! **
    xb_out = model(xb.float())
    
    if verbose:
        print("model out pre loss", xb_out.shape)
        #print('xb_out', xb_out)
        print("xb_out:",xb_out.shape)
        print("yb:",yb.shape)
        print("yb.long:",yb.long().shape)
    
    loss = loss_func(xb_out, yb.float()) # for MSE/regression
    # __FOOTNOTE 2__
    
    if opt is not None: # if opt
        loss.backward()
        opt.step()
        opt.zero_grad()

    return loss.item(), len(xb)
    
def train_step(model, train_dl, loss_func, device, opt):
    '''
    Execute 1 set of batched training within an epoch
    '''
    # Set model to Training mode
    model.train()
    tl = [] # train losses
    ns = [] # batch sizes, n
    
    # loop through train DataLoader
    for xb, yb in train_dl:
        # put on GPU
        xb, yb = xb.to(device),yb.to(device)
        
        # provide opt so backprop happens
        t, n = loss_batch(model, loss_func, xb, yb, opt=opt)
        
        # collect train loss and batch sizes
        tl.append(t)
        ns.append(n)
    
    # average the losses over all batches    
    train_loss = np.sum(np.multiply(tl, ns)) / np.sum(ns)
    
    return train_loss

def val_step(model, val_dl, loss_func, device):
    '''
    Execute 1 set of batched validation within an epoch
    '''
    # Set model to Evaluation mode
    model.eval()
    with torch.no_grad():
        vl = [] # val losses
        ns = [] # batch sizes, n
        
        # loop through validation DataLoader
        for xb, yb in val_dl:
            # put on GPU
            xb, yb = xb.to(device),yb.to(device)

            # Do NOT provide opt here, so backprop does not happen
            v, n = loss_batch(model, loss_func, xb, yb)

            # collect val loss and batch sizes
            vl.append(v)
            ns.append(n)

    # average the losses over all batches
    val_loss = np.sum(np.multiply(vl, ns)) / np.sum(ns)
    
    return val_loss
    
def fit(epochs, model, loss_func, opt, train_dl, val_dl,device,patience=1000):
    '''
    Fit the model params to the training data, eval on unseen data.
    Loop for a number of epochs and keep train of train and val losses 
    along the way
    '''
    # keep track of losses
    train_losses = []    
    val_losses = []
    
    # loop through epochs
    for epoch in range(epochs):
        # take a training step
        train_loss = train_step(model,train_dl,loss_func,device,opt)
        train_losses.append(train_loss)

        # take a validation step
        val_loss = val_step(model,val_dl,loss_func,device)
        val_losses.append(val_loss)
        
        print(f"E{epoch} | train loss: {train_loss:.3f} | val loss: {val_loss:.3f}")

    return train_losses, val_losses


def run_model(train_dl,val_dl,model,device,
              lr=0.01, epochs=50, 
              lossf=None,opt=None
             ):
    '''
    Given train and val DataLoaders and a NN model, fit the mode to the training
    data. By default, use MSE loss and an SGD optimizer
    '''
    # define optimizer
    if opt:
        optimizer = opt
    else: # if no opt provided, just use SGD
        optimizer = torch.optim.SGD(model.parameters(), lr=lr)
    
    # define loss function
    if lossf:
        loss_func = lossf
    else: # if no loss function provided, just use MSE
        loss_func = torch.nn.MSELoss()
    
    # run the training loop
    train_losses, val_losses = fit(
                                epochs, 
                                model, 
                                loss_func, 
                                optimizer, 
                                train_dl, 
                                val_dl, 
                                device)

    return train_losses, val_losses

5. Запуск моделей

Для начала испытаем на наших октамерных последовательностях линейный подход.

# use GPU if available
DEVICE = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

# get the sequence length from the first seq in the df
seq_len = len(train_df['seq'].values[0])

# create Linear model object
model_lin = DNA_Linear(seq_len)
model_lin.to(DEVICE) # put on GPU

# run the model with default settings!
lin_train_losses, lin_val_losses = run_model(
    train_dl, 
    val_dl, 
    model_lin,
    DEVICE
)

Соберём данные о потерях train и val и набросаем график:

def quick_loss_plot(data_label_list,loss_type="MSE Loss",sparse_n=0):
    '''
    For each train/test loss trajectory, plot loss by epoch
    '''
    for i,(train_data,test_data,label) in enumerate(data_label_list):    
        plt.plot(train_data,linestyle='--',color=f"C{i}", label=f"{label} Train")
        plt.plot(test_data,color=f"C{i}", label=f"{label} Val",linewidth=3.0)

    plt.legend()
    plt.ylabel(loss_type)
    plt.xlabel("Epoch")
    plt.legend(bbox_to_anchor=(1,1),loc='upper left')
    plt.show()
    

lin_data_label = (lin_train_losses,lin_val_losses,"Lin")
quick_loss_plot([lin_data_label])

Обучение линейной модели и оценка кривых потерь

На первый взгляд, такой подход привнёс в обучение немного.

Теперь попробуем метод свёрточных нейросетей и построим график для него.

seq_len = len(train_df['seq'].values[0])

# create Linear model object
model_cnn = DNA_CNN(seq_len)
model_cnn.to(DEVICE) # put on GPU

# run the model with default settings!
cnn_train_losses, cnn_val_losses = run_model(
    train_dl, 
    val_dl, 
    model_cnn,
    DEVICE
)

cnn_data_label = (cnn_train_losses,cnn_val_losses,"CNN")
quick_loss_plot([lin_data_label,cnn_data_label])

Кривые потерь для линейного метода и метода свёрточных нейросетей (CNN)

Очевидно, что CNN может уловить паттерн, а линейный метод не может. Прогоним несколько последовательностей и посмотрим, что получится.

# oracle dict of true score for each seq
oracle = dict(mer8[['seq','score']].values)

def quick_seq_pred(model, desc, seqs, oracle):
    '''
    Given a model and some sequences, get the model's predictions
    for those sequences and compare to the oracle (true) output
    '''
    print(f"__{desc}__")
    for dna in seqs:
        s = torch.tensor(one_hot_encode(dna)).unsqueeze(0).to(DEVICE)
        pred = model(s.float())
        actual = oracle[dna]
        diff = pred.item() - actual
        print(f"{dna}: pred:{pred.item():.3f} actual:{actual:.3f} ({diff:.3f})")

def quick_8mer_pred(model, oracle):
    seqs1 = ("poly-X seqs",['AAAAAAAA', 'CCCCCCCC','GGGGGGGG','TTTTTTTT'])
    seqs2 = ("other seqs", ['AACCAACA','CCGGTGAG','GGGTAAGG', 'TTTCGTTT'])
    seqsTAT = ("with TAT motif", ['TATAAAAA','CCTATCCC','GTATGGGG','TTTATTTT'])
    seqsGCG = ("with GCG motif", ['AAGCGAAA','CGCGCCCC','GGGCGGGG','TTGCGTTT'])
    TATGCG =  ("both TAT and GCG",['ATATGCGA','TGCGTATT'])

    for desc,seqs in [seqs1, seqs2, seqsTAT, seqsGCG, TATGCG]:
        quick_seq_pred(model, desc, seqs, oracle)
        print()
        
# Ask the trained Linear model to make 
# predictions for some 8-mers
quick_8mer_pred(model_lin, oracle)

# prints:
# __poly-X seqs__
# AAAAAAAA: pred:23.230 actual:20.000 (3.230)
# CCCCCCCC: pred:13.582 actual:17.000 (-3.418)
# GGGGGGGG: pred:7.006 actual:14.000 (-6.994)
# TTTTTTTT: pred:17.767 actual:11.000 (6.767)

# __other seqs__
# AACCAACA: pred:18.818 actual:18.875 (-0.057)
# CCGGTGAG: pred:12.205 actual:15.125 (-2.920)
# GGGTAAGG: pred:13.826 actual:15.125 (-1.299)
# TTTCGTTT: pred:14.815 actual:12.125 (2.690)

# __with TAT motif__
# TATAAAAA: pred:22.146 actual:27.750 (-5.604)
# CCTATCCC: pred:16.931 actual:25.875 (-8.944)
# GTATGGGG: pred:12.141 actual:24.000 (-11.859)
# TTTATTTT: pred:18.266 actual:22.125 (-3.859)

# __with GCG motif__
# AAGCGAAA: pred:16.736 actual:8.125 (8.611)
# CGCGCCCC: pred:12.346 actual:6.250 (6.096)
# GGGCGGGG: pred:7.907 actual:4.375 (3.532)
# TTGCGTTT: pred:12.839 actual:2.500 (10.339)

# __both TAT and GCG__
# ATATGCGA: pred:15.664 actual:15.875 (-0.211)
# TGCGTATT: pred:14.771 actual:13.625 (1.146)

Из приведённых выше примеров видно, что линейная модель занижает прогноз для последовательности с большим количеством G и завышает его для последовательности с большим количеством T. По-видимому, она просто не замечает, что последовательности с GCG имеют необычно низкие баллы, а последовательности с TAT — необычно высокие. Но, поскольку в линейной модели невозможно учесть различия между GCG и GAG, она просто понижает балл для последовательности с G. Из нашей схемы подсчёта баллов известно, что это не так. Дело тут не в том, что G вредны или, если хотите, субвитальны, а в том, что вредны именно мотивы GCG.

# Ask the trained CNN model to make 
# predictions for some 8-mers
quick_8mer_pred(model_cnn, oracle)

# prints: 
# __poly-X seqs__
# AAAAAAAA: pred:19.722 actual:20.000 (-0.278)
# CCCCCCCC: pred:16.620 actual:17.000 (-0.380)
# GGGGGGGG: pred:13.771 actual:14.000 (-0.229)
# TTTTTTTT: pred:10.767 actual:11.000 (-0.233)

# __other seqs__
# AACCAACA: pred:18.530 actual:18.875 (-0.345)
# CCGGTGAG: pred:14.925 actual:15.125 (-0.200)
# GGGTAAGG: pred:14.900 actual:15.125 (-0.225)
# TTTCGTTT: pred:11.789 actual:12.125 (-0.336)

# __with TAT motif__
# TATAAAAA: pred:26.154 actual:27.750 (-1.596)
# CCTATCCC: pred:24.321 actual:25.875 (-1.554)
# GTATGGGG: pred:22.870 actual:24.000 (-1.130)
# TTTATTTT: pred:20.581 actual:22.125 (-1.544)

# __with GCG motif__
# AAGCGAAA: pred:8.689 actual:8.125 (0.564)
# CGCGCCCC: pred:6.531 actual:6.250 (0.281)
# GGGCGGGG: pred:5.029 actual:4.375 (0.654)
# TTGCGTTT: pred:3.055 actual:2.500 (0.555)

# __both TAT and GCG__
# ATATGCGA: pred:15.178 actual:15.875 (-0.697)
# TGCGTATT: pred:12.853 actual:13.625 (-0.772)

CNN лучше улавливает разницу между тримерными мотивами! При этом он хорошо переваривает последовательности и с мотивами и без них.

6. Проверка прогнозов модели на тестовом наборе данных

В любой задаче машинного обучения при оценке результатов важно убедиться, что модель может делать хорошие прогнозы по тестовым наборам, с которыми она никогда не сталкивалась в ходе обучения. Разница между реальными тестовыми последовательностями и прогнозами модели лучше всего видна на сравнительных графиках.

import altair as alt
from sklearn.metrics import r2_score

def parity_plot(model_name,df,r2):
    '''
    Given a dataframe of samples with their true and predicted values,
    make a scatterplot.
    '''
    plt.scatter(df['truth'].values, df['pred'].values, alpha=0.2)
    
    # y=x line
    xpoints = ypoints = plt.xlim()
    plt.plot(xpoints, ypoints, linestyle='--', color='k', lw=2, scalex=False, scaley=False)

    plt.ylim(xpoints)
    plt.ylabel("Predicted Score",fontsize=14)
    plt.xlabel("Actual Score",fontsize=14)
    plt.title(f"{model_name} (r2:{r2:.3f})",fontsize=20)
    plt.show()
    
def alt_parity_plot(model,df, r2):
    '''
    Make an interactive parity plot with altair
    '''
    chart = alt.Chart(df).mark_circle(size=100,opacity=0.4).encode(
        alt.X('truth:Q'),
        alt.Y('pred:Q'),
        tooltip=['seq:N']
    ).properties(
        title=f'{model} (r2:{r2:.3f})'
    ).interactive()
    
    chart.save(f'alt_out/parity_plot_{model}.html')
    display(chart)

def parity_pred(models, seqs, oracle,alt=False):
    '''Given some sequences, get the model's predictions '''
    dfs = {} # key: model name, value: parity_df
    
    for model_name,model in models:
        print(f"Running {model_name}")
        data = []
        for dna in seqs:
            s = torch.tensor(one_hot_encode(dna)).unsqueeze(0).to(DEVICE)
            actual = oracle[dna]
            pred = model(s.float())
            data.append([dna,actual,pred.item()])
        df = pd.DataFrame(data, columns=['seq','truth','pred'])
        r2 = r2_score(df['truth'],df['pred'])
        dfs[model_name] = (r2,df)
        
        #plot parity plot
        if alt: # make an altair plot
            alt_parity_plot(model_name, df, r2)
        else:
            parity_plot(model_name, df, r2)
            
# generate plots
seqs = test_df['seq'].values
models = [
    ("Linear", model_lin),
    ("CNN", model_cnn)
]
parity_pred(models, seqs, oracle)

Разница между реальными тестовыми последовательностями и прогнозами модели.

Сравнение графиков показывает, насколько точно модель прогнозирует отдельные последовательности. У идеальной модели все графики точно лягут на диагональ y=x. Это покажет, что модель точно прогнозирует последовательности. Отклонения от диагонали y=x означают, что прогнозы модели завышены или занижены.

Как мы видим, линейная модель может в какой-то мере прогнозировать линию тренда тестовых последовательностей, но обилие последовательностей в нижней и верхней областях распределения (при наличии мотивов) легко сбивают эту модель с толку.

CNN работает гораздо лучше и ожидаемо предсказывает значения, близкие к реальным, ведь архитектура нашего CNN использует ядра тримеров для поиска ключевых мотивов вдоль последовательности.

Впрочем, CNN — тоже не идеал. Можно было бы обучать его дольше или корректировать гиперпараметры, но цель здесь не в совершенстве. Мы сильно упростили задачу по сравнению с реальными регулятивными грамматиками. И я подумала, что для интерактивного выявления последовательностей, с которыми модели не справляются, было бы интересно использовать библиотеку визуализации Altair:

Заметьте, что вне диагонали наши последовательности, как правило, имеют множественные экземпляры мотивов! При подсчёте баллов мы просто давали последовательности ± поправки за один мотив. В действительности же разумно было бы несколько таких «бонусов» за несколько мотивов. В своём примере я просто даю «бонус» за 1 экземпляр мотива, но функцию расчёта баллов можно было бы составить иначе.

В любом случае здорово, когда модель замечает множественные экземпляры и прогнозирует их важность. Кажется, здесь я немного обманула её… Впрочем, R2 0,95 — это уже не так плохо.

7. Визуализация свёрточных фильтров

При обучении моделей CNN полезно визуализировать свёрточные фильтры первого слоя, чтобы лучше понять процесс обучения. При работе с изображениями свёрточные фильтры первого слоя часто изучают паттерны границ, цвета или текстуры — основные элементы изображения, при объединении которых мы получаем более сложные признаки.

В случае ДНК мы исходим из того, что свёрточные фильтры можно рассматривать как сканеры мотивов. По аналогии с матрицей весовых коэффициентов позиций для визуализации логотипов последовательностей, свёрточный фильтр похож на матрицу, показывающую определённый паттерн ДНК. Однако вместо того, чтобы быть точной последовательностью, он может содержать некоторую неопределённость относительно того, какие нуклеотиды показываются в той или иной части паттерна. Некоторые позиции могут быть строго определены (например, в позиции 2 всегда есть A; информативность высокая), в других же с почти равной вероятностью может быть множество нуклеотидов (энтропия высокая, информативность низкая).

Вычисления на скрытых слоях нейросетей могут быть очень сложными, и не каждый свёрточный фильтр будет представлять собой очевидный паттерн. Однако фильтры всё же выявляют информативные для объяснения прогноза паттерны.

Ниже приведены некоторые функции для визуализации свёрточных фильтров первого слоя как в виде необработанной тепловой карты, так и в виде логотипа.

import logomaker

def get_conv_layers_from_model(model):
    '''
    Given a trained model, extract its convolutional layers
    '''
    model_children = list(model.children())
    
    # counter to keep count of the conv layers
    model_weights = [] # we will save the conv layer weights in this list
    conv_layers = [] # we will save the actual conv layers in this list
    bias_weights = []
    counter = 0 

    # append all the conv layers and their respective weights to the list
    for i in range(len(model_children)):
        # get model type of Conv1d
        if type(model_children[i]) == nn.Conv1d:
            counter += 1
            model_weights.append(model_children[i].weight)
            conv_layers.append(model_children[i])
            bias_weights.append(model_children[i].bias)

        # also check sequential objects' children for conv1d
        elif type(model_children[i]) == nn.Sequential:
            for child in model_children[i]:
                if type(child) == nn.Conv1d:
                    counter += 1
                    model_weights.append(child.weight)
                    conv_layers.append(child)
                    bias_weights.append(child.bias)

    print(f"Total convolutional layers: {counter}")
    return conv_layers, model_weights, bias_weights

def view_filters(model_weights, num_cols=8):
    model_weights = model_weights[0]
    num_filt = model_weights.shape[0]
    filt_width = model_weights[0].shape[1]
    num_rows = int(np.ceil(num_filt/num_cols))
    
    # visualize the first conv layer filters
    plt.figure(figsize=(20, 17))

    for i, filter in enumerate(model_weights):
        ax = plt.subplot(num_rows, num_cols, i+1)
        ax.imshow(filter.cpu().detach(), cmap='gray')
        ax.set_yticks(np.arange(4))
        ax.set_yticklabels(['A', 'C', 'G','T'])
        ax.set_xticks(np.arange(filt_width))
        ax.set_title(f"Filter {i}")

    plt.tight_layout()
    plt.show()
    
conv_layers, model_weights, bias_weights = get_conv_layers_from_model(model_cnn)
view_filters(model_weights)

Возможно, это полезная информация, но обычно люди любят визуализировать последовательности с некоторой неопределённостью в виде логотипов мотивов, откладывая позиции мотива по X, а вероятности появления нуклеотидов в каждой из них — по Y. Часто для удобства визуализации эти вероятности преобразуются в биты (информацию).

Для преобразования необработанных свёрточных фильтров в визуальные матрицы весовых коэффициентов позиций обычно собирают активации фильтров: применяют веса фильтра вдоль кодируемой прямо последовательности и измеряют активацию фильтра (а также то, насколько хорошо веса соответствуют последовательности).

Матрицы весовых коэффициентов фильтра, соответствующие близкому совпадению с последовательностью, будут высокоактивированными (то есть дадут более высокие оценки совпадения). Собирая подпоследовательности ДНК, которые дают самые высокие оценки активации, мы можем создать весовую матрицу «высокоактивированных последовательностей» для каждого фильтра и, таким образом, визуализировать свёрточный фильтр в виде логотипа мотива.

Схема, показывающая степень сбора активированных участков последовательности и их преобразования в логотипы мотивов для фильтра свёртки

def get_conv_output_for_seq(seq, conv_layer):
    '''
    Given an input sequeunce and a convolutional layer, 
    get the output tensor containing the conv filter 
    activations along each position in the sequence
    '''
    # format seq for input to conv layer (OHE, reshape)
    seq = torch.tensor(one_hot_encode(seq)).unsqueeze(0).permute(0,2,1).to(DEVICE)

    # run seq through conv layer
    with torch.no_grad(): # don't want as part of gradient graph
        # apply learned filters to input seq
        res = conv_layer(seq.float())
        return res[0]
    

def get_filter_activations(seqs, conv_layer,act_thresh=0):
    '''
    Given a set of input sequences and a trained convolutional layer, 
    determine the subsequences for which each filter in the conv layer 
    activate most strongly. 
    
    1.) Run seq inputs through conv layer. 
    2.) Loop through filter activations of the resulting tensor, saving the
            position where filter activations were > act_thresh. 
    3.) Compile a count matrix for each filter by accumulating subsequences which
            activate the filter above the threshold act_thresh
    '''
    # initialize dict of pwms for each filter in the conv layer
    # pwm shape: 4 nucleotides X filter width, initialize to 0.0s
    num_filters = conv_layer.out_channels
    filt_width = conv_layer.kernel_size[0]
    filter_pwms = dict((i,torch.zeros(4,filt_width)) for i in range(num_filters))
    
    print("Num filters", num_filters)
    print("filt_width", filt_width)
    
    # loop through a set of sequences and collect subseqs where each filter activated
    for seq in seqs:
        # get a tensor of each conv filter activation along the input seq
        res = get_conv_output_for_seq(seq, conv_layer)

        # for each filter and it's activation vector
        for filt_id,act_vec in enumerate(res):
            # collect the indices where the activation level 
            # was above the threshold
            act_idxs = torch.where(act_vec>act_thresh)[0]
            activated_positions = [x.item() for x in act_idxs]

            # use activated indicies to extract the actual DNA
            # subsequences that caused filter to activate
            for pos in activated_positions:
                subseq = seq[pos:pos+filt_width]
                #print("subseq",pos, subseq)
                # transpose OHE to match PWM orientation
                subseq_tensor = torch.tensor(one_hot_encode(subseq)).T

                # add this subseq to the pwm count for this filter
                filter_pwms[filt_id] += subseq_tensor            
            
    return filter_pwms

def view_filters_and_logos(model_weights,filter_activations, num_cols=8):
    '''
    Given some convolutional model weights and filter activation PWMs, 
    visualize the heatmap and motif logo pairs in a simple grid
    '''
    model_weights = model_weights[0].squeeze(1)
    print(model_weights.shape)

    # make sure the model weights agree with the number of filters
    assert(model_weights.shape[0] == len(filter_activations))
    
    num_filts = len(filter_activations)
    num_rows = int(np.ceil(num_filts/num_cols))*2+1 
    # ^ not sure why +1 is needed... complained otherwise
    
    plt.figure(figsize=(20, 17))

    j=0 # use to make sure a filter and it's logo end up vertically paired
    for i, filter in enumerate(model_weights):
        if (i)%num_cols == 0:
            j += num_cols

        # display raw filter
        ax1 = plt.subplot(num_rows, num_cols, i+j+1)
        ax1.imshow(filter.cpu().detach(), cmap='gray')
        ax1.set_yticks(np.arange(4))
        ax1.set_yticklabels(['A', 'C', 'G','T'])
        ax1.set_xticks(np.arange(model_weights.shape[2]))
        ax1.set_title(f"Filter {i}")

        # display sequence logo
        ax2 = plt.subplot(num_rows, num_cols, i+j+1+num_cols)
        filt_df = pd.DataFrame(filter_activations[i].T.numpy(),columns=['A','C','G','T'])
        filt_df_info = logomaker.transform_matrix(filt_df,from_type='counts',to_type='information')
        logo = logomaker.Logo(filt_df_info,ax=ax2)
        ax2.set_ylim(0,2)
        ax2.set_title(f"Filter {i}")

    plt.tight_layout()
    

# just use some seqs from test_df to activate filters
some_seqs = random.choices(seqs, k=3000)

filter_activations = get_filter_activations(some_seqs, conv_layers[0],act_thresh=1)
view_filters_and_logos(model_weights,filter_activations)

# prints:
# Num filters 32
# filt_width 3
# torch.Size([32, 4, 3])

Из конкретно этого обучения CNN видно, что несколько фильтров уловили сильные мотивы TAT и GCG, в то время как другие фильтры сфокусировались и на других паттернах.

О важности визуализации свёрточных фильтров для интерпретации модели нет единого мнения. В глубоких моделях с несколькими свёрточными слоями свёрточные фильтры могут комбинироваться в скрытых слоях более сложными способами, поэтому фильтры первого слоя сами по себе могут быть не такими информативными (Ку и Эдди, 2019]). Большая часть этой области с тех пор перешла к механизмам внимания и другим методам объяснения результатов. Однако, если вам интересно визуализировать фильтры как потенциальные мотивы, эти функции могут быть полезны!

8. Итоги

В руководстве даются основы построения моделей CNN в PyTorch для обработки последовательностей ДНК. Практическая часть работы не отражает реалии биологических сигналов. Наш метод подсчёта баллов, скорее, моделирует наличие регуляторных мотивов в очень коротких последовательностях, которые легко понять человеку. Поведение PyTorch соответствовало ожиданиям. В этом маленьком примере мы проверили способность базовой CNN со сдвиговыми фильтрами к прогнозу расчёта баллов. CNN показала себя лучше, чем линейный метод, выявивший только абсолютные положения нуклеотидов (без учёта местных особенностей).

Подробнее о применении CNN к ДНК в естественных условиях можно узнать из следующих работ:

• DeepBind: Alipanahi et al 2015.
• DeepSea: Zhou and Troyanskaya 2015.
• Basset: Kelley et al 2016.

Надеюсь, эта статья поможет с PyTorch всем, кто интересуется биологией и ДНК, но мало знаком с машинным обучением.

9. Примечания

Первое

В нашем примере свёрточной нейросети используется одномерный свёрточный слой. Поскольку ДНК — не двухмерная картинка, достаточно прогнать Conv1D по всей длине последовательности, не сканируя «вниз» и «вверх». (Сдвиг фильтра «вверх» и «вниз» неприменим при кодировании матриц ДНК с одним активным состоянием. Нет смысла отделять ряды A и C от G и T — все четыре ряда должны быть строго представлены в последовательности).

Однако мне однажды понадобился инструмент анализа от keras, и тогда я нашла скрипт преобразования pytorch2keras. Скрипт преобразования умеет работать только со слоями Conv2d и выдаёт ошибки по слоям Conv1d в модели.

Если он понадобится и вам, вот пример переформатирования CNN в Conv2D, который обеспечивает беспроблемное сканирование, как для Conv1D:

class DNA_CNN_2D(nn.Module):
    def __init__(self,
                 seq_len,
                 num_filters=31,
                 kernel_size=3,
                ):
        super().__init__()
        self.seq_len = seq_len
        
        self.conv_net = nn.Sequential(
            nn.Conv2d(1, num_filters, kernel_size=(4,kernel_size)),
            # ^^ changed from 4 to 1 channel, moved 4 to kernel_size
            nn.ReLU(),
            nn.Flatten(),
            nn.Linear(num_filters*(seq_len-kernel_size+1), 1)
        ) 

    def forward(self, xb):
        # reshape view to batch_ssize x 4channel x seq_len
        # permute to put channel in correct order
        xb = xb.permute(0,2,1).unsqueeze(1)
        # ^^ Conv2D input fix
        
        #print(xb.shape)
        out = self.conv_net(xb)
        return out

Второе

Если ваша задача связана с классификацией, а не с регрессией, возможно, вы предпочтёте использовать CrossEntropyLoss. Учтите, что CrossEntropyLoss требуется формат, отличный от MSELoss. Попробуйте такой вариант:

loss = loss_func(xb_out, yb.long().squeeze (1))

А пока учёные моделируют ДНК, мы научим вас разрабатывать системы ИИ, которые решают широкий круг проблем:

• Профессия Data Scientist (24 месяца)
• Профессия Fullstack-разработчик на Python (15 месяцев)

Чтобы увидеть все курсы, кликните по баннеру: