Поиск последовательности Шайна-Дальгарно

Поиск мотивов осуществлялся в геноме бактерии Defluviitoga tunisiensis.

Подготовка последовательностей

Приведенный ниже скрипт принадлежит Дмитрию Звездину (c изменениями).

In [3]:
from Bio import SeqIO
import pandas as pd
import numpy as np
from IPython.display import Image
In [8]:
with open("features.gb", "r") as file:
    table = SeqIO.read(file, "gb")
coordinates = []
with open("coordinates.tsv", "w") as file_write:
    for feature in table.features[1:]:
        start, end, strand = int(feature.location.start), int(feature.location.end), feature.location.strand
        try:
            locus_tag = feature.qualifiers["locus_tag"][0]
        except KeyError:
            locus_tag = f"unknown_locus_tag_{start}_{end}"
        if (feature.type == "CDS") and ("pseudo" not in feature.qualifiers.keys()) and (end - start > 300) and (feature.qualifiers["product"][0] != "hypothetical protein"):
            coordinates.append("\t".join([str(start), str(end), str(strand), locus_tag]))
    print(*coordinates, sep="\n", file=file_write)

Файл в формате genbank был скачан из записи по ссылке. Далее из него были получены координаты тех CDS, которые не имеют статус "unreviewed" и при этом имеют длину более 300 нуклеотидов. Координаты (start, end) и расположение на цепи (1 - на прямой, -1 - на обратной) были сохранены в файл coordinates.tsv.

In [23]:
coordinates_table = pd.read_csv("coordinates.tsv", sep="\t", names=["start", "end", "strand", "name"])
with open("Defluviitoga_tunisiensis.fasta", "r") as fasta_file:
    fasta = SeqIO.read(fasta_file, "fasta")
clear_sequences = []
with open("clear_sequences.fasta", "w") as file_write:
    
    for i in range(301):
        name = coordinates_table.loc[i, "name"]
        if str(coordinates_table.loc[i, "strand"]) == "1":
            start = int(coordinates_table.loc[i, "start"])
            sequence = str(fasta.seq[start - 20: start])
        else:
            start = int(coordinates_table.loc[i, "end"])
            sequence = str(fasta.seq[start: start + 20].reverse_complement())
        clear_sequences.append(f">{name}")
        clear_sequences.append(sequence)
    print(*clear_sequences, sep="\n", file=file_write)

По координатам из файла из генома Defluviitoga tunisiensis были выбраны 300 последовательностей длины 20 нуклеотидов в upstream генов. Последовательности были сохранены в файл clear_sequences.fasta и разбиты на 2 набора: 100 последовательностей для MEME и 200 для FIMO.

In [29]:
with open('clear_sequences.fasta', 'r') as f, open('meme_sequences.fasta', 'w') as f1, open('fimo_sequences.fasta', 'w') as f2:
  i = 0
  line = f.readlines()
  while i < 202:
    f1.write(line[i])
    i += 1
  if i >= 202:
    for i in range(202, 602):
      f2.write(line[i])
In [24]:
for feature in table.features[1:]:
    try:
        if ("RNA" in feature.type) and ("16S" in feature.qualifiers["product"][0]):
            print(feature.location)
    except KeyError:
        pass

[938625:940137](+)
[1334080:1335592](+)
[1873559:1875071](+)

Далее был взят участок в 20 нуклеотидов 3' конца 16 рРНК и добавлен в файл для построения PWM.

In [25]:
print(str(fasta.seq[940117:940137].complement()))

CGCGGACCTAGTGGAGGAAA

Построение позиционной весовой матрицы

Программа MEME была использована для построения PWM:

meme meme_sequences.fasta -dna -minw 5 -maxw 8 -nmotifs 1

Выдача meme в формате HTML доступна по ссылке.

Был найден мотив длиной 6 нуклеотидов в 84 из 101 последовательностях из выборки.

Ниже представлено LOGO мотива.

In [30]:
Image("logo1.png")
Out[30]:

Проверка PWM

Для проверки была испольована программа fimo:

fimo --thresh 0.0005 ./meme_out/meme.txt fimo_sequences.fasta

По умолчанию программа fimo использует порог e-value, равный 10е-4. Однако при таком пороге не обнаруживалось ни одного появления мотива в последовательностях, в связи с чем порог был увеличен до 5*10е-4.

С выдачей в tsv формате можно ознакомиться по ссылке. Было найдено 79 последовательностей из 200 (39.5%). Вывод: PWM обладает средним качеством.