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
121 changes: 121 additions & 0 deletions TCF4.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,121 @@
from Bio import SeqIO
from Bio.SeqUtils import GC
Comment on lines +1 to +2
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Работа с SeqIO.parse реализована хорошо, единственное, он избыточен в паре мест.


def filter_fastq(input_path: str, quality_threshold: int, output_filename="final_filtered.fastq", gc_bounds=(40, 60), length_bounds=(50, 350)):
filename = input_path
records = SeqIO.parse(filename, "fastq")
Comment on lines +5 to +6
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Зачем вводить переменную filename? Так ты просто заставляешь два имени ссылаться на один и тот же объект. Можно было оставить input_path:

Suggested change
filename = input_path
records = SeqIO.parse(filename, "fastq")
records = SeqIO.parse(input_path, "fastq")

###quality filter
good_reads = (rec for rec in records if min(rec.letter_annotations["phred_quality"]) >= quality_threshold)
result_quality = SeqIO.write(good_reads, "good_quality.fastq", "fastq")
result_quality_GC = SeqIO.parse("good_quality.fastq", "fastq")
###GC content filter
min_gc_content = gc_bounds[0]
max_gc_content = gc_bounds[1]
Comment on lines +12 to +13
Copy link
Copy Markdown

@angrygeese angrygeese Mar 4, 2024

Choose a reason for hiding this comment

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

gc_bounds , по условию задания, может, в том числе, принимать только верхний порог в виде int. В этом случае у тебя здесь упадёт код.

Suggested change
min_gc_content = gc_bounds[0]
max_gc_content = gc_bounds[1]
if isinstance(gc_bounds, int):
gc_bounds = (0, gc_bounds)
min_gc_content, max_gc_content = gc_bounds

GC_quality_filt = []

for sequence in result_quality_GC:
if min_gc_content <= GC(sequence.seq) <= max_gc_content:
GC_quality_filt.append(sequence)

result_quality = SeqIO.write(GC_quality_filt, "good_quality_GC.fastq", "fastq")
result_quality_GC_length = SeqIO.parse("good_quality_GC.fastq", "fastq")

##length filter
filtered_GC_quality_length = []

for sequence in result_quality_GC_length:
if len(sequence.seq) >= length_bounds[0] and len(sequence.seq) <= length_bounds[1]:
Comment on lines +26 to +27
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Suggested change
for sequence in result_quality_GC_length:
if len(sequence.seq) >= length_bounds[0] and len(sequence.seq) <= length_bounds[1]:
length_min, length_max = length_bounds
for sequence in result_quality_GC_length:
seq_len = len(sequence.seq)
if seq_len >= length_min and seq_len <= length_max:

filtered_GC_quality_length.append(sequence)

result_quality = SeqIO.write(filtered_GC_quality_length, output_filename, "fastq")

print(result_quality)
Comment on lines +5 to +32
Copy link
Copy Markdown

@angrygeese angrygeese Mar 4, 2024

Choose a reason for hiding this comment

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

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

Suggested change
filename = input_path
records = SeqIO.parse(filename, "fastq")
###quality filter
good_reads = (rec for rec in records if min(rec.letter_annotations["phred_quality"]) >= quality_threshold)
result_quality = SeqIO.write(good_reads, "good_quality.fastq", "fastq")
result_quality_GC = SeqIO.parse("good_quality.fastq", "fastq")
###GC content filter
min_gc_content = gc_bounds[0]
max_gc_content = gc_bounds[1]
GC_quality_filt = []
for sequence in result_quality_GC:
if min_gc_content <= GC(sequence.seq) <= max_gc_content:
GC_quality_filt.append(sequence)
result_quality = SeqIO.write(GC_quality_filt, "good_quality_GC.fastq", "fastq")
result_quality_GC_length = SeqIO.parse("good_quality_GC.fastq", "fastq")
##length filter
filtered_GC_quality_length = []
for sequence in result_quality_GC_length:
if len(sequence.seq) >= length_bounds[0] and len(sequence.seq) <= length_bounds[1]:
filtered_GC_quality_length.append(sequence)
result_quality = SeqIO.write(filtered_GC_quality_length, output_filename, "fastq")
print(result_quality)
def filter_fastq_corr(input_path: str, quality_threshold: int, output_filename="final_filtered_corr.fastq", gc_bounds=(40, 60), length_bounds=(50, 350)):
if isinstance(gc_bounds, int):
gc_bounds = (0, gc_bounds)
min_gc_content, max_gc_content = gc_bounds
length_min, length_max = length_bounds
records = SeqIO.parse(input_path, "fastq")
good_reads = []
for rec in records:
rec_len, rec_gc_content = len(rec.seq), GC(rec.seq)
if (
min(rec.letter_annotations["phred_quality"]) >= quality_threshold
and max_gc_content >= rec_gc_content >= min_gc_content
and length_max >= rec_len >= length_min
):
good_reads.append(rec)
result_quality = SeqIO.write(good_reads, output_filename, "fastq")
print(result_quality)


#filter_fastq("example_fastq.fastq", 15)


from abc import ABC, abstractmethod

class InvalidInputError(ValueError):
pass

class BiologicalSequence(ABC, str):
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

При таком наследовании BiologicalSequence и его наследники переймут все методы класса str, которые ты не переопределишь. Например, у экземпляров твоего DNASequence есть методы count и transtale.

Suggested change
class BiologicalSequence(ABC, str):
class BiologicalSequence(ABC):

@abstractmethod
def __init__(self, seq):
self.seq = seq
Comment on lines +43 to +45
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Насколько я успел посмотреть, в абстрактных классах не используют связку __init__ и @abstractmethod.

Обычно делают так: в родительском, абстрактном классе объявляют переменную класса __slots__. Это такой способ зарезервировать память под набор атрибутов. Затем в __init__ класса-наследника их определяют.

Suggested change
@abstractmethod
def __init__(self, seq):
self.seq = seq
__slots__ = ('seq',)

Если вдруг станет интересно, что можно посмотреть:

  1. Исходный код абстрактного класса для последовательностей в Python.
  2. Особенности использования __slots__ детально разобраны тут.

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

upd: Я всё-таки встретил связку __init__ и @abstractmethod.


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

def __getitem__(self, index):
return self.seq[int(index)]
Comment on lines +47 to +51
Copy link
Copy Markdown

@angrygeese angrygeese Mar 6, 2024

Choose a reason for hiding this comment

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

Я бы выполнил это задание в таком ключе, определив каждый метод как абстрактный, без конкретной реализации:

Suggested change
def __len__(self):
return len(self.seq)
def __getitem__(self, index):
return self.seq[int(index)]
@abstractmethod
def __len__(self):
raise NotImplementedError
@abstractmethod
def __getitem__(self, index):
raise NotImplementedError

Питон допускает реализацию в абстрактном методе, но это нужно для отдельных случаев с множественным наследованием.


def __repr__(self):
return __str__(self.seq)
Comment on lines +53 to +54
Copy link
Copy Markdown

@angrygeese angrygeese Mar 3, 2024

Choose a reason for hiding this comment

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

Если ты хочешь вызвать метод __str__ у self.seq, то синтаксис должен быть таким:

    def __repr__(self):
        return self.seq.__str__()

Или можно сделать это же через встроенную функцию str:

    def __repr__(self):
        return str(self.seq)

У тебя получилась смесь того и другого, поэтому код упадёт с ошибкой.

Suggested change
def __repr__(self):
return __str__(self.seq)
def __repr__(self):
return self.seq.__str__()

Но, вообще, лучше оставить метод в абстрактном классе без определения:

Suggested change
def __repr__(self):
return __str__(self.seq)
@abstractmethod
def __repr__(self):
raise NotImplementedError

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

Когда пишите код в комментах PR, добавляйте python после ``` - так синтаксис будет подсвечен как надо


def check_nucleic_acid(self):
unique_chars = set(self.seq)
nucleotides_dna = set('ATGCatgc')
nucleotides_rna = set('AUGCaugc')
if unique_chars <= nucleotides_dna:
seq = 'dna'
elif unique_chars <= nucleotides_rna:
seq = 'rna'
else:
raise InvalidInputError()
return seq_type
Comment on lines +56 to +66
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Идея правильная - метод сравнивает уникальные элементы с "алфавитом" конкретного полимера. Но:

  1. unique_chars = set(self.seq) предполагает, что self.seq у тебя итерируемый, а ты этого явно нигде не проверяешь. Т.е. если создать экземпляр, передав DNASequence число 1000, эта проверка упадёт не по предусмотренному тобой сценарию.
  2. получается, ты внутри абстрактного класса реализуешь проверку только для двух конкретных наследников. Хотя, например, такая логика применима и к аминокислотным, но для них этот метод работать не будет. Логично реализовать проверку, запрашивая у класса его алфавит. При этом, в complement у тебя всё отлично, и метод опирается на атрибуты класса, для которого его вызывают.
    Здесь можно сделать похожим образом - добавить классу атрибут self.alphabet и сделать как-то так:
Suggested change
def check_nucleic_acid(self):
unique_chars = set(self.seq)
nucleotides_dna = set('ATGCatgc')
nucleotides_rna = set('AUGCaugc')
if unique_chars <= nucleotides_dna:
seq = 'dna'
elif unique_chars <= nucleotides_rna:
seq = 'rna'
else:
raise InvalidInputError()
return seq_type
def check_alphabet(self) -> bool:
return set(self.seq).issubset(self.alphabet)

Теперь все потомки BiologicalSequence будут при проверке сравнивать свой собственный алфавит и свою собственную последовательность.


class NucleicAcidSequence(BiologicalSequence):
def __init__(self, seq):
super().__init__(seq)
Copy link
Copy Markdown

@angrygeese angrygeese Mar 6, 2024

Choose a reason for hiding this comment

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

Благодаря этому моменту я разобрался в деталях работы super(). И что им, оказывается, можно задействовать смысловой код из абстрактного метода класса-родителя.

self.check_nucleic_acid()
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Неплохо было бы реализовать на основе check_nucleic_acid проверку, что пользователь пытается создать правильный класс для последовательности. В текущей реализации я могу вызвать DNASequence('AUGC'), получить экземпляр класса, и никто ни на что ругаться не будет.

Для этого можно объявить отдельный класс-прослойку между BiologicalSequence и NucleotideSequence, назовём его BiologicalSequenceInitializer. В этом классе последовательно реализуем геттер и сеттер для алфавита и для самой последовательности. Заодно мы сможем переделать метод check_nucleic_acid в check_alphabet и задействовать его внутри сеттера последовательности:

class BiologicalSequence(ABC):
    
    __slots__ = ('seq', 'alphabet')

    ...

class BiologicalSequenceInitializer(BiologicalSequence):
    def __init__(self):
        self._alphabet = None  # защищённый атрибут
        self._seq = None    # защищённый атрибут

    # инициализируем алфавит:
    
    ## геттер для алфавита
    @property
    def alphabet(self) -> set:
        return self._alphabet

    ## сеттер для алфавита
    @alphabet.setter
    def alphabet(self, alphabet: set):
        if not isinstance(alphabet, set):
            raise TypeError('Must be a set.')
        self._alphabet = alphabet

    # инициализируем последовательность:
        
    ## геттер для последовательности
    @property
    def seq(self) -> str:
        return self._seq

    ## сеттер для последовательности
    @seq.setter
    def seq(self, seq: str):
        if not (isinstance(seq, str)):
            raise TypeError('Must be a string.')
        if not self.check_alphabet(seq):  # задействуем проверку алфавита внутри сеттера
            raise TypeError(f'Sequence does not match alphabet for {repr(self)}. Please, use only: {self._alphabet}.')
        self._seq = seq

    ...

    def check_alphabet(self, seq) -> bool:
        return set(seq).issubset(self.alphabet)

Теперь достаточно просто задать в наследниках атрибуты self.seq и selt.alphabet присваиванием:

class NucleicAcidSequence(BiologicalSequenceInitializer):
    
    ...
    
class RNASequence(NucleicAcidSequence):
    def __init__(self, seq: str):
        self.alphabet = set('AUGCaugc')
        self.seq = seq

и будут задействоваться сеттеры из BiologicalSequenceInitializer. Кроме того, сеттер сработает и при попытке определить атрибут вне класса, т.е.:

dna_sequence1 = DNASequence('ATGC')
dna_sequence1.seq = 'AUGC'

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

Вот это конечно объект комментов, ахах..... Мощь

self.length = len(self.seq)

def complement(self):
list_input = list(self.seq)
for i in range(len(self.seq)):
if list_input[i] in self.complement_dict:
list_input[i] = self.complement_dict[list_input[i]]
Comment on lines +77 to +78
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Огонь - метод опирается на атрибуты класса.

return "".join(list_input)

class DNASequence(NucleicAcidSequence):
complement_dict = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G', 'a': 't', 't': 'a', 'g': 'c', 'c': 'g'}
def __init__(self, seq):
super().__init__(seq)
#self.complement()

def transcribe(self):
list_input = list(self.seq)
for i in range(len(self.seq)):
if (list_input[i] == 'T'):
list_input[i] = 'U'
elif (list_input[i] == 't'):
list_input[i]='u'
return "".join(list_input)
Comment on lines +87 to +94
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Есть способ сократить код:

Suggested change
def transcribe(self):
list_input = list(self.seq)
for i in range(len(self.seq)):
if (list_input[i] == 'T'):
list_input[i] = 'U'
elif (list_input[i] == 't'):
list_input[i]='u'
return "".join(list_input)
def transcribe(self) -> RNASequence:
return RNASequence(str(self).translate(str.maketrans(('Tt'), ('Uu'))))


class RNASequence(NucleicAcidSequence):
complement_dict = {'A': 'U', 'U': 'A', 'G': 'C', 'C': 'G', 'a': 'u', 'u': 'a', 'g': 'c', 'c': 'g'}
def __init__(self, seq):
super().__init__(seq)
#self.complement()

class AminoAcidSequence(BiologicalSequence):
def __init__(self, seq):
self.seq = seq
Comment on lines +102 to +104
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Здесь есть маленькая деталь: ни в классе-родителе, ни в самом AminoAcidSequence не реализована проверка алфавита последовательности на корректность. Поэтому я могу, например, создать экземпляр AminoAcidSequence, передав конструктору класса число 1000. И потом запросить для него частоту аминокислот.


def amino_acid_frequency(self):
"""Calculates molecular weight of a protein
Comment on lines +106 to +107
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Не из той функции взята докстринга.

Arguments:
- seq (str) 1-letter coded protein sequence
Return:
- int, molecular weight (g/mol) rounded to integer"""
freq_dict = {}
for letter in self.seq:
if letter in freq_dict:
freq_dict[letter] += 1
else:
freq_dict[letter] = 1
Comment on lines +114 to +117
Copy link
Copy Markdown

Choose a reason for hiding this comment

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

Логика абсолютно правильная, есть способ записать её несколько короче.

У словарей в Python есть метод get, он принимает два аргумента - ключ и любой аргумент. Если ключ в словаре есть, возвращает значение из словаря по этому ключу, если нет - возвращает второй аргумент.

Suggested change
if letter in freq_dict:
freq_dict[letter] += 1
else:
freq_dict[letter] = 1
freq_dict[letter] = freq_dict.get(letter, 0) + 1

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

Его используют сильно реже, но тут кажется отличный кейс для get

for letter in freq_dict:
freq_dict[letter] = round(freq_dict[letter] / len(self.seq) * 100, 2)
return freq_dict