Skip to content
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
146 changes: 146 additions & 0 deletions RHNO1.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction
from abc import ABC, abstractmethod


class FastQFilter:
def __init__(
self, input_file, output_file, min_length=0, min_quality=0, min_gc=0, max_gc=100
):
self.input_file = input_file
self.output_file = output_file
self.min_length = min_length
self.min_quality = min_quality
self.min_gc = min_gc
self.max_gc = max_gc

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

все логично и понятно, что приходит на вход, хорошо, что прописаны автоматические минимальные значения


def filter_fastq(self):
with open(self.output_file, "w") as output_handle:
for record in SeqIO.parse(self.input_file, "fastq"):

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

тут можно было бы добавить тест на существование файла с последовательностью, который берется на вход, если его нет, то можно выводить ошибку о его отсутствии

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
for record in SeqIO.parse(self.input_file, "fastq"):
for record in SeqIO.parse(self.input_file, "fastq"):
if not self.input_file():
raise FileNotFoundError(f"Файл '{self.input_file}' не найден!")

if self._passes_filter(record):
SeqIO.write(record, output_handle, "fastq")

def _passes_filter(self, record):
if len(record.seq) < self.min_length:
return False

if min(record.letter_annotations["phred_quality"]) < self.min_quality:
return False

gc_content = gc_fraction(record.seq)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

gc_content = GC(record.seq)
Для подсчета гц-состава можно использовать такой подход. Импортируется через from Bio.SeqUtils import GC

if gc_content < self.min_gc or gc_content > self.max_gc:
return False

return True


class BiologicalSequence(ABC):
@abstractmethod
def __len__(self):
pass

@abstractmethod
def __getitem__(self, index):
pass

@abstractmethod
def __str__(self):
pass

@abstractmethod
def check_alphabet(self):
pass


class NucleicAcidSequence(BiologicalSequence):
complement_dict = {"A": "", "T": "", "G": "", "C": ""}

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

не совсем понятно для чего этот словарь, если я правильно понимаю, он дальше меняется, но вероятно, тогда стоит задать пустой словарь

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Не очень поняла зачем здесь эта строка. Кажется, что она определяет словарь, который фактически не используется в классе NucleicAcidSequence, так как метод complement() использует словарь, определенный в дочерних классах DNASequence и RNASequence. Но лучше, конечно, чтобы в этом классе все определялось, а не в наследуемых.


def __init__(self, sequence):
self.sequence = sequence

def __len__(self):
return len(self.sequence)

def __getitem__(self, index):
return self.sequence[index]

def __str__(self):
return self.sequence

def check_alphabet(self):
valid_alphabet = set("ATGC")

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Вижу проверку на ДНК и на аминокислоты ниже, круто! Но не вижу на РНК, лучше ее тоже добавить

return set(self.sequence) <= valid_alphabet

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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


def complement(self):
complemented_sequence = [
self.complement_dict.get(base, base) for base in self.sequence
]
return self._create_instance("".join(complemented_sequence))

def _create_instance(self, sequence):
return self.__class__(sequence)

def gc_content(self):
gc_count = sum(base in "GCgc" for base in self.sequence)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Круто, что проверяется регистр! Можно для проверки алфавита тоже такое добавить

return gc_count / len(self.sequence)


class DNASequence(NucleicAcidSequence):
complement_dict = {"A": "T", "T": "A", "G": "C", "C": "G"}

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Вероятно этот комментарий супер вкусовщина. Более приятно читается словарь, когда он построчный. Но это сугубо вкусовщина

Suggested change
complement_dict = {"A": "T", "T": "A", "G": "C", "C": "G"}
complement_dict = {"A": "T",
"T": "A",
"G": "C",
"C": "G"}

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

следующие словарики я бы также подпарвил


def transcribe(self):
transcription_rules = {"A": "U", "T": "A", "G": "C", "C": "G"}
transcribed_sequence = "".join(
transcription_rules.get(base, base) for base in self.sequence
)
return RNASequence(transcribed_sequence)


class RNASequence(NucleicAcidSequence):
complement_dict = {"A": "U", "U": "A", "G": "C", "C": "G"}

def codons(self):
return [self.sequence[i : i + 3] for i in range(0, len(self.sequence), 3)]

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Может оказаться так, что последовательность не кратна 3, поэтому на такой случай стоит что-то придумать, например показывать пользователю ошибку

Suggested change
return [self.sequence[i : i + 3] for i in range(0, len(self.sequence), 3)]
if len(self.sequence) % 3 != 0:
raise ValueError("Sequence length is not a multiple of 3")
else:
return [self.sequence[i : i + 3] for i in range(0, len(self.sequence), 3)]



class AminoAcidSequence(BiologicalSequence):
def __init__(self, sequence):
self.sequence = sequence

def __len__(self):
return len(self.sequence)

def __getitem__(self, index):
return self.sequence[index]

def __str__(self):
return self.sequence

def check_alphabet(self):

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Можно еще прописать вывод ошибок, если подаются символы, отличные от ДНК/РНК/аминокислот, без них у меня работает вообще с любой последовательностью.

valid_alphabet = set("ACDEFGHIKLMNPQRSTVWY")
return set(self.sequence) <= valid_alphabet

def get_molecular_weight(self):
molecular_weights = {
"A": 89.09,
"C": 121.16,
"D": 133.10,
"E": 147.13,
"F": 165.19,
"G": 75.07,
"H": 155.16,
"I": 131.18,
"K": 146.19,
"L": 131.18,
"M": 149.21,
"N": 132.12,
"P": 115.13,
"Q": 146.15,
"R": 174.20,
"S": 105.09,
"T": 119.12,
"V": 117.15,
"W": 204.23,
"Y": 181.19,
}

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

очень красивый словарик, легко читаемый

return sum(molecular_weights[aa] for aa in self.sequence)

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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