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

Для поиска мотивов была выбрана бактерия Rhodococcus ruber.

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

from Bio import SeqIO
import pandas as pd
import numpy as np
from IPython.display import Image

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, длиной более 300 нуклеотидов, продуктами которых, не являются гипотетические белки, их координаты были сохранены в файл.

coordinates_table = pd.read_csv("coordinates.tsv", sep="\t", names=["start", "end", "strand", "name"])
with open("Rhodococcus_ruber.fasta", "r") as fasta_file:
    fasta = SeqIO.read(fasta_file, "fasta")
clear_sequences = []
with open("pr7/clear_sequences.fasta", "w") as file_write:
    
    for i in range(300):
        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)

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

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

[745882:747404](+)
[971143:972666](+)
[2287488:2289010](+)
[3211379:3212902](-)

print(str(fasta.seq[747384:747404].complement()))

CGCCGACCTAGTGGAGGAAA

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

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

Программа MEME была использована для построения PWM: meme meme_input.fasta -dna -minw 5 -maxw 8 -nmotifs 1

С выдачей в HTML формате можно ознакомится по ссылке. В результате был найден мотиф длиной 8 нуклеотидов, похожий на консенсусную последовательность Шайна — Дальгарно: AGGAGG. Он был обнаружен в 60 из 101 последовательности выборки. Ниже представлено LOGO мотива.

Image("logo1.png")

Проверка PWM

Для проверки была испольована программа fimo: fimo ./meme_out/meme.txt fimo_input.fasta

С выдачей в tsv фрмате можно ознакомится по ссылке. Было найдено 19 последовательностей из 200, что может говорить о не очень хорошем качестве полученной PWM.