import pandas as pd
import numpy as np

Введение

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

Выбор генома

Для исследования был выбран геном Methanocaldococcus jannaschii (ссылка на REBASE, ссылка на геном).
Данные по рестриктазам второго типа данной археи доступны в файле M_jannaschii.txt.

Mja_df = pd.read_csv('M_jannaschii.txt', sep='\t', header=0)

Mja_df[Mja_df['Specificity'] > ''].sort_values(by=['Specificity'])

	ORF	Gene	Most similar	Specificity	Name
17	1208	R	226 aa conserved hypothetical protein	CCGG	MjaVIP
18	???	M	DNA methylase N-4/N-6 domain protein	CCGG	M.MjaVI
9	984	R	restriction enzyme, type II R/M system 1	CTAG	MjaI
10	985	M	modification methylase, type II R/M system 1	CTAG	M.MjaI
4	598	M	modification methylase, type II R/M system 2	GATC	M.MjaIII
6	600	R	restriction enzyme, type II R/M system 2	GATC	MjaIII
25	1448	M	modification methylase, type II R/M system	GGNCC	M.MjaII
26	1449	R	370 aa hypothetical protein	GGNCC	MjaII
29	1498	M	modification methylase, type II R/M system	GTAC	M.MjaV
31	1500	R	230 aa hypothetical protein	GTAC	MjaV
21	1327	R	245 aa hypothetical protein	GTNNAC	MjaIV
22	???	M	modification methylase, type II R/M system	GTNNAC	M.MjaIV

Mja_sites = (Mja_df.dropna())['Specificity'].unique()

print(Mja_sites)

['GATC' 'CTAG' 'CCGG' 'GTNNAC' 'GGNCC' 'GTAC']

Получение списка сайтов

Список потенциальных сайтов рестрикции был получен из таблицы.

RE_df = pd.read_csv('TypeII_REs.tsv', sep='\t', header=0)

with open('sites.txt', 'w') as f:
    f.write(' '.join( (RE_df[ RE_df['Recognition site'] != '-' ])['Recognition site'].unique() ))

Отбор недопредставленных сайтов

Для полученных потенциальных сайтов была найдена их представленность в геноме Methanocaldococcus jannaschii при помощи контраста, вычисленного методом С.Карлина и соавторов.
(команда: cbcalc M_jannaschii.fasta --burge -o out.tsv)

cbcalc_df = pd.read_csv('out.tsv', sep='\t', header=0)
del cbcalc_df['Sequence ID']
cbcalc_df = cbcalc_df.sort_values(by=['O/E ratio (BCK)'])

Представленность сайтов, для которых для данной археии в REBASE присутствуют рестриктазы, приведена в таблице ниже. Можно видеть, что все они, как и ожидалось, недопредставлены (O/E ratio < 1), но степень несоответствия ожидаемому в значительной степени разнится.

cbcalc_df[cbcalc_df['Site'].isin(Mja_sites)]

	Site	Observed	Expected (BCK)	O/E ratio (BCK)	Total
58	CTAG	90	1555.32	0.058	1664915
23	GATC	253	2238.50	0.113	1664915
78	GTNNAC	168	1299.95	0.129	1664887
187	GTAC	334	1073.90	0.311	1664915
100	GGNCC	817	1254.99	0.651	1664901
74	CCGG	573	739.95	0.774	1664915

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

cbcalc_df.head(n=6)

	Site	Observed	Expected (BCK)	O/E ratio (BCK)	Total
96	GGCCNNNNNGGCC	0	0.42	0.000	1664790
32	GCGCGC	0	3.57	0.000	1664887
6	CGATCG	0	1.60	0.000	1664887
193	GTCGAC	0	1.50	0.000	1664887
194	GGCCGGCC	0	0.32	0.000	1664859
58	CTAG	90	1555.32	0.058	1664915

Для дальнейшего рассмотрения были отобраны сайты длиной более одного нуклеотида с представленностью менее 0.8 (поскольку это было названо "разумным значением" в задании). Таких сайтов нашлось 24.

chosen_sites = cbcalc_df[
    (0 < cbcalc_df['O/E ratio (BCK)']) &
    (cbcalc_df['O/E ratio (BCK)'] < 0.8) &
    (cbcalc_df['Site'].apply(len) > 1)
]

chosen_sites

	Site	Observed	Expected (BCK)	O/E ratio (BCK)	Total
58	CTAG	90	1555.32	0.058	1664915
23	GATC	253	2238.50	0.113	1664915
78	GTNNAC	168	1299.95	0.129	1664887
187	GTAC	334	1073.90	0.311	1664915
156	GCGC	119	371.63	0.320	1664915
206	GGCC	689	1610.10	0.428	1664915
170	CATATG	103	225.90	0.456	1664887
117	GGGCCC	30	52.09	0.576	1664887
115	CTCGAG	6	9.22	0.651	1664887
100	GGNCC	817	1254.99	0.651	1664901
67	GGNNCC	656	1001.45	0.655	1664887
38	CCCGGG	52	76.82	0.677	1664887
11	CGRYCG	20	28.53	0.701	1664887
179	ACGCGT	4	5.52	0.724	1664887
49	CCTAGG	3	4.07	0.738	1664887
147	YGGCCR	102	137.21	0.743	1664887
183	CCCGC	110	145.97	0.754	1664901
197	CTTAAG	228	297.50	0.766	1664887
153	GGWCC	405	523.20	0.774	1664901
74	CCGG	573	739.95	0.774	1664915
207	TGGCCA	55	70.49	0.780	1664887
164	GCNGC	1393	1775.71	0.784	1664901
2	GASTC	529	672.05	0.787	1664901
166	AGGCCT	42	52.96	0.793	1664887

Отбор экспериментально проверенных эндонуклеаз

Из эндонуклеаз, распознающих отобранные сайты, были взяты те, активность которых была проверена экспериментально. Список их названий в REBASE (в формате list-файла для seqret) доступен по ссылке, всего 125 белков.

chosen_RE = RE_df[
    RE_df['Recognition site'].isin(chosen_sites['Site']) &
    (RE_df['Putative'] == 'no')
]

chosen_RE

	#:REBASE name	Sequence AC	System type	Protein type	Recognition site	Length	Locus	Protein ID	Uniprot AC	GI	Putative	Complex name
40	AatI	NEBC13653	Type II	R	AGGCCT	274	-	-	-	-	no	-
152	AbrI	X62690	Type II	R	CTCGAG	222	-	CAA44567.1	P26919	38661	no	-
505	ApaI	HQ327698	Type II	R	GGGCCC	355	-	ADO24182.1	E3VXA1	308229533	no	-
997	AvaII	BA000019	Type II	R	GGWCC	230	-	BAB72890.1	Q8YYB7	17130279	no	-
1068	BalI	D82028	Type II	R	TGGCCA	260	-	BAA11515.1	P71102	1620877	no	-
...	...	...	...	...	...	...	...	...	...	...	...	...
23583	TliI	CP006670	Type II	R	CTCGAG	281	OCC_13281	EHR77447.1	H3ZRP1	374741014	no	-
23734	TspMI	EF426476	Type II	R	CCCGGG	379	-	ABO38150.1	A4L311	133753294	no	-
23777	Tvu2HI	CP039710	Type II	R	GGCC	316	FA954_10130	QCV55937.1	-	-	no	-
24482	XcyI	M98768	Type II	R	CCCGGG	333	-	AAA27609.1	P30773	155372	no	-
24573	XhoI	NEBC13645	Type II	R	CTCGAG	246	-	-	-	-	no	-

125 rows × 12 columns

RE_name = np.array(chosen_RE['#:REBASE name'])

with open('RE_list.txt', 'w') as f:
    for name in RE_name:
        f.write(f'TypeII_REs.fasta:{name}\n')

Поиск гомологов отобранных эндонуклеаз

Для поиска гомологов отобранных эндонуклеаз был проведён tblastn с геномом археи в качестве базы и последовательностями рестриктаз к качастве query. Выдача blast доступна в файле

blast_df = pd.read_csv('blast.out', sep='\t', header=0)
blast_df['Recognition site'] = np.array((RE_df[RE_df['#:REBASE name'].isin(blast_df['qseqid'])])['Recognition site'])
blast_df['in Mja (REBASE)'] = np.array( blast_df['qseqid'].isin(Mja_df['Name']) )
blast_df = blast_df.sort_values(by=['sstart'])

blast_df

	qseqid	qstart	qend	qcovs	sstart	send	evalue	Recognition site	in Mja (REBASE)
5	HpyAIII	5	273	97	531176	532030	3.000000e-44	GATC	False
10	MjaIII	1	290	100	531176	532045	2.000000e-161	GATC	True
3	DpnII	6	287	98	531179	532036	5.000000e-57	GATC	False
14	R1.Ssu11318I	6	303	97	531182	532039	5.000000e-47	GATC	False
2	BstXII	20	308	93	531182	532045	3.000000e-81	GATC	False
19	R1.SsuDAT1I	6	303	97	531182	532039	5.000000e-47	GATC	False
18	R1.Ssu8074I	6	303	97	531182	532039	5.000000e-47	GATC	False
17	R1.Ssu4961I	6	303	97	531182	532039	5.000000e-47	GATC	False
16	R1.Ssu4109I	6	303	97	531182	532039	1.000000e-39	GATC	False
15	R1.Ssu2479I	6	303	97	531182	532039	5.000000e-47	GATC	False
20	R2.Ssu11318I	19	299	93	531194	532036	4.000000e-72	GATC	False
23	R2.Ssu4961I	19	299	93	531194	532036	4.000000e-72	GATC	False
25	R2.SsuDAT1I	19	299	93	531194	532036	4.000000e-72	GATC	False
21	R2.Ssu2479I	19	299	93	531194	532036	4.000000e-72	GATC	False
7	LlaDCHI	15	299	94	531194	532036	9.000000e-71	GATC	False
6	LlaAI	15	299	94	531194	532036	1.000000e-70	GATC	False
24	R2.Ssu8074I	19	299	93	531194	532036	4.000000e-72	GATC	False
22	R2.Ssu4109I	19	299	93	531194	532036	6.000000e-72	GATC	False
28	Ssu220I	4	206	98	531419	532039	2.000000e-56	GATC	False
26	Ssu211I	2	204	99	531419	532039	2.000000e-53	GATC	False
27	Ssu212I	2	203	99	531422	532039	9.000000e-53	GATC	False
8	MjaI	1	222	100	915786	916451	5.000000e-120	CTAG	True
1	BfaIB	10	192	93	915855	916421	4.000000e-46	CTAG	False
0	BfaIA	19	200	89	915876	916421	1.000000e-44	CTAG	False
13	MthZI	25	190	82	915882	916400	8.000000e-45	CTAG	False
4	Hpy299IX	31	225	85	1153230	1152631	2.000000e-51	CCGG	False
11	MjaIV	1	245	100	1277076	1277810	3.000000e-154	GTNNAC	True
9	MjaII	1	370	100	1418572	1419681	0.000000e+00	GGNCC	True
29	TdeIII	33	228	76	1418677	1419321	9.000000e-17	GGNCC	False
12	MjaV	1	230	100	1472915	1472226	8.000000e-135	GTAC	True

Находок достаточно много, с хорошим e-value и достаточно большим покрытием. Можно видеть, что часто на приблизительно один локус выравнивается множество находок, из которых одна - с покрытием ровно 100% - сооветствует белку Mja согласно REBASE, а остальные имеют покрытие немного меньше, но всё ещё хорошее e-value, распознают ту же специфичную последовательность и, вероятно, просто являются близкими гомологами этого белка.

При этом из отмеченных в REBASE blast не нашёл семь белков:

Mja_df[(Mja_df['Specificity'] > '') & (~Mja_df['Name'].isin(blast_df['qseqid']))].sort_values(by=['Specificity'])

	ORF	Gene	Most similar	Specificity	Name
17	1208	R	226 aa conserved hypothetical protein	CCGG	MjaVIP
18	???	M	DNA methylase N-4/N-6 domain protein	CCGG	M.MjaVI
10	985	M	modification methylase, type II R/M system 1	CTAG	M.MjaI
4	598	M	modification methylase, type II R/M system 2	GATC	M.MjaIII
25	1448	M	modification methylase, type II R/M system	GGNCC	M.MjaII
29	1498	M	modification methylase, type II R/M system	GTAC	M.MjaV
22	???	M	modification methylase, type II R/M system	GTNNAC	M.MjaIV

Шесть из них являются метилазами, а не рестриктазами, соответственно, находились вне области рассмотрения. А вот вместо рестриктазы MjaVIP нашёлся другой белок, распознающий ту же последовательность (CCGG): Hpy299IX. Покрытие всего 85 процентов, то есть это лишь гомолог, а не сам белок. Сам же белок не нашёлся, поскольку не ещё не проходил экспериментальную проверку:

RE_df[RE_df['#:REBASE name'] == 'MjaVIP']

	#:REBASE name	Sequence AC	System type	Protein type	Recognition site	Length	Locus	Protein ID	Uniprot AC	GI	Putative	Complex name
15540	MjaVIP	L77117	Type II	R	CCGG	226	MJ_1208.1	AAB99221.1	P81326	2826378	yes	-

Вывод

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