E-box

Расскажу немного о сигнале E-box. E-box - это энхансерная (enhancer box) нуклеотидная последовательность, встречающаяся в некоторых эукариотических промоторах. Обычно этот сигнал связывется белками с характерным мотивом "спираль-петля-спираль". Пример последовательности E-box'а: CAXXTG.

Источник.

import numpy as np
import pandas as pd
import sympy as sy
import seaborn as sns
from sklearn.model_selection import train_test_split

from Bio import Entrez
from Bio import SeqIO

import os

Entrez.email = "sample@gmail.com"
stream = Entrez.efetch(db='nucleotide', id='Y16535', rettype='gb', retmode='text')

record = SeqIO.read(stream, "genbank")
stream.close()
seq = str(record.seq)
seq[323:329]

'CACTTG'

PWM

В рамках этого задания будем анализировать окрестность ATG в геноме человека.

Теперь соберём данные для обучения и тестирования PWM-модели. Для этого возьмём некоторое количество генов с первой хромосомы. Данные для обучения и теста - контекст первого ATG в гене, негативный контроль - окружение не первого ATG в гене.

genes_set = pd.read_csv('../../../../Lab/DBs/china/fromrawgenome/bed_diced.bed', sep='\t',
                       names=['chr',
                              'start',
                              'end',
                              'strand',
                              'annot',
                              'name',
                              'smth',
                              'source',
                              'seq'])

genes_set.iloc[:1000, :]

	chr	start	end	strand	annot	name	smth	source	seq
0	chr1	11869.0	14409.0	+	transcribed_unprocessed_pseudogene	DDX11L1	DDX11L1	gencode	TTAACTTGCCGTCAGCCTTTTCTTTGACCTCTTCTTTCTGTTCATG...
1	chr1	14404.0	29570.0	-	unprocessed_pseudogene	WASH7P	WASH7P	gencode	TTCTGCTCAGTTCTTTATTGATTGGTGTGCCGTTTTCTCTGGAAGC...
2	chr1	17368.0	17436.0	-	miRNA	hsa-mir-6859-1	hsa-mir-6859-1	GB_snomirna	CTACAGAGGCGACATGGGGGTCAGGCAAGCTGACACCCGCTGTCCT...
3	chr1	17369.0	17436.0	-	miRNA	MIR6859-1	MIR6859-1	gencode	TACAGAGGCGACATGGGGGTCAGGCAAGCTGACACCCGCTGTCCTG...
4	chr1	29554.0	31109.0	+	lncRNA	MIR1302-2HG	MIR1302-2HG	gencode	TGCACACGGCTCCCATGCGTTGTCTTCCGAGCGTCAGGCCGCCCCT...
...	...	...	...	...	...	...	...	...	...
995	chr1	24919345.0	24919416.0	-	miRNA	MIR6731	MIR6731	gencode	TGGGGAGAGTGGGGAATAGAGGCAGGTGGTTGGCACCTGGAGCTTC...
996	chr1	24961345.0	24963097.0	+	lncRNA	RUNX3-AS1	RUNX3-AS1	gencode	AAAAGGGAAGTCATCTTCCGTTGGTCTATGGGGACCCTGAGGTAGT...
997	chr1	24968423.0	24970865.0	-	lncRNA	AL445471.1	AL445471.1	gencode	GGCATCATCTAGATTTTACTATTTCTTTCCATCCTTTACAATGTCC...
998	chr1	25023502.0	25023586.0	+	miRNA	hsa-mir-4425	hsa-mir-4425	GB_snomirna	GTGCTTTACATGAATGGTCCCATTGAATCCCAACAGCTTTGCGAAG...
999	chr1	25023503.0	25023586.0	+	miRNA	MIR4425	MIR4425	gencode	TGCTTTACATGAATGGTCCCATTGAATCCCAACAGCTTTGCGAAGT...

1000 rows × 9 columns

start_atg_seqs = []
downstream_atg_seqs = []

for i, row in genes_set.iloc[:, :].iterrows():
    if row['strand'] == '+':
        start_coord = row['seq'].find('ATG')
        
        if start_coord > 6 and len(row['seq'][start_coord-7:start_coord+6]) == 13:
            start_atg_seqs.append(row['seq'][start_coord-7:start_coord+6])
        
            down_coord = row['seq'][start_coord+1:].find('ATG')
            if row['seq'][down_coord-7:down_coord+6] != '':
                downstream_atg_seqs.append(row['seq'][down_coord-7:down_coord+6])

train_starts, test_starts = train_test_split(start_atg_seqs, test_size=0.4)

with open('train_seq.fa', 'w') as fa:
    print(*train_starts, sep='\n', file=fa)

Мы получили наши выборки, теперь настало время расчитать частоты букв в каждой из позиций. GC-состав будем считать равным $0.444$.

Для удобства оформим позиционно-весовую матрицу в класс PWM.

class PWM():
    def __init__(ma, alignment, gc_cont, pseudocount=0.01):
        
        cont_dee = {'A': (1 - gc_cont)/2,
                    'T': (1 - gc_cont)/2,
                    'G': gc_cont/2,
                    'C': gc_cont/2}
    
        matrix = []
    
        for i in range(len(alignment[0])):
            freqs = {'A': 0, 'T': 0, 'G': 0, 'C': 0}
        
            for seq in alignment:
                freqs[seq[i]] += 1
            
            for key in freqs:
                freqs[key] = np.log((freqs[key] + pseudocount)/cont_dee[key]/len(alignment))
            
            matrix.append(list(freqs.values()))
        
        matrix = np.array(matrix).T
    
        df = pd.DataFrame({key: matrix[chi] for chi, key in enumerate(['A', 'T', 'G', 'C'])}).T
        df.columns = np.arange(1, len(alignment[0]) + 1)
        
        ma.df = df
        
        ma.min_result = sum(min(col) for col in matrix)
        ma.max_result = sum(max(col) for col in matrix)
    
    def weights(ma):
        return ma.df
    
    def predict(ma, seq):
        outp = 0
        for j, liter in enumerate(seq):
            outp += ma.df.iloc[:, j][liter]
            proba = (outp - ma.min_result)/(ma.max_result - ma.min_result)
        return proba

pwm = PWM(train_starts, 0.444)
pwm.weights()

	1	2	3	4	5	6	7	8	9	10	11	12	13
A	-0.173972	-0.029164	-0.267759	-0.201197	0.052022	0.004591	-0.106610	1.280135	-13.105904	-13.105904	-0.121881	-0.146057	-0.148177
T	-0.238946	-0.252464	-0.100297	-0.174941	-0.348533	-0.379845	-0.643636	-13.105904	1.280135	-13.105904	-0.334231	-0.167461	-0.268025
G	0.264145	0.232365	0.175052	0.204534	0.318181	0.222222	0.250365	-12.880960	-12.880960	1.505078	0.426722	0.115162	0.311118
C	0.151700	0.052422	0.200792	0.184200	-0.075774	0.132037	0.362206	-12.880960	-12.880960	-12.880960	-0.033294	0.216087	0.097073

Итак, мы расчитали нашу pwm. Попробуем визуализировать её чуть понятнее - построим тепловую карту:

sns.heatmap(pwm.weights(), cmap='coolwarm')

<Axes: >

Теперь "прогоним" нашу матрицу на тренировочной и тестовой выборках:

positive_scores = [pwm.predict(pos_seq) for pos_seq in test_starts]
negative_scores = [pwm.predict(neg_seq) for neg_seq in downstream_atg_seqs]

np.mean(positive_scores), np.mean(negative_scores)

(0.9982553614468179, 0.4187664680451599)

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

Построенное лого: