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
190 changes: 190 additions & 0 deletions PRICKLE4.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,190 @@
import os
from abc import ABC, abstractmethod
from Bio import SeqIO
from Bio import SeqUtils


class BiologicalSequence(ABC, str):

Choose a reason for hiding this comment

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

Здесь не стоило наследоваться от строки, все же это абстрактный класс, ну и вообще, как оказалось, наследование от базового класса здесь не предполагалось.


def __init__(self, seq):
self.seq = seq
Comment on lines +7 to +10
Copy link

Choose a reason for hiding this comment

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

Абстрактный класс должен содержать только абстрактные методы, которые определяют, что должно быть в дочерних классах

Comment on lines +8 to +10

Choose a reason for hiding this comment

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

Как сама автор кода указывает ниже, здесь стоило оставить только абстрактные методы. Тот же комментарий к методам ниже.


#да, Никита уже сказал, что так лучше не делать, а оставить сдесь только абстрактные методы,
# но я зацепилась за идею не повторять проверки, раз суть для всех одна, только алфавиты разные,
# а переделать уже не успеваю
def check_alphabet(self):
if set(self.seq) <= self.alphabet:
raise UnexpectedSymbolInSeqError

# не поняла, не доделала
@abstractmethod
def to_print_seq(self):
return self.seq
pass
Comment on lines +19 to +23

Choose a reason for hiding this comment

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

Здесь, если я правильно понимаю, предполагалась реализация "вывода на печать в удобном виде"? В таком случае, стоило это сделать через str. Вообще, как я понимаю, предполагалось, что наш объект имеет атрибут, в котором лежит строка с последовательностью. Тогда можно было бы сделать так:

Suggested change
# не поняла, не доделала
@abstractmethod
def to_print_seq(self):
return self.seq
pass
def __str__(self):
return self.seq

Как я понимаю, здесь логика была иной, учитывая наследование от строки. Такой объект тоже можно было бы перевести в тип строка, например, итерируясь по нему, добавляя элементы в list, а затем переведя лист в строку. Это не очень осмысленно, но если нужен именно тип строка, так можно было бы сделать.
Повторюсь только, что в абстрактном классе делать этого не стоило.


Comment on lines +20 to +24

Choose a reason for hiding this comment

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

Абстрактные методы в абстрактном классе определяют некий интерфейс, который будет унаследован дочерними классами. Тут подразумевалось, что в родительском классе Biological sequence вы скажите, что дочерние классы должны уметь считать длину, получать индексы и прочее. Например, про длину, это можно записать как

Suggested change
@abstractmethod
def to_print_seq(self):
return self.seq
pass
@abstractmethod
def __len__(self):
pass

И потом переопределить нашу заглушку в дочернем, написав что-то осмысленное вместо pass
Как-то так)


class UnexpectedSymbolInSeqError(ValueError):
pass
Comment on lines +26 to +27

Choose a reason for hiding this comment

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

Приятно, что сделана специальная ошибка.



class AminoAcidSequence(BiologicalSequence):
alphabet = {'M', 'O', 'v', 'D', 'f', 'N', 'c', 'A', 'R', 'W', 'I', 'm', 'L', 's', 'H', 'q', 'w', 'V', 'n', 'i',
'g', 'F', 'S', 'e', 'l', 'U', 'P', 'Q', 'K', 'Y', 'u', 'y', 'd', 'h', 'k', 'r', 't', 'G', 'o', 'E',
'p', 'T', 'C', 'a'}
masses = {'A': 71.08, 'R': 156.2, 'N': 114.1, 'D': 115.1, 'C': 103.1, 'E': 129.1, 'Q': 128.1, 'G': 57.05, 'H': 137.1,
'I': 113.2, 'L': 113.2, 'K': 128.2, 'M': 131.2, 'F': 147.2, 'P': 97.12, 'S': 87.08, 'T': 101.1,
'W': 186.2, 'Y': 163.2, 'V': 99.13, 'U': 168.05, 'O': 255.3, 'a': 71.08, 'r': 156.2, 'n': 114.1, 'd': 115.1,
'c': 103.1, 'e': 129.1, 'q': 128.1, 'g': 57.05, 'h': 137.1, 'i': 113.2, 'l': 113.2, 'k': 128.2, 'm': 131.2,
'f': 147.2, 'p': 97.12, 's': 87.08, 't': 101.1, 'w': 186.2, 'y': 163.2, 'v': 99.13, 'u': 168.05, 'o': 255.3}
Comment on lines +31 to +38

Choose a reason for hiding this comment

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

Такое количество алфавитов кажется избыточным. Если вы хотите проверить последовательность и определить массу, вы можете использовать всего один словарь)

Comment on lines +31 to +38

Choose a reason for hiding this comment

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

Тут стоит отметить, что хотя решение не является неправильным, alphabet здесь не обязательная переменная (кстати, переменные стоит писать капсом). Ведь можно использовать set(masses.keys()). Это, конечно, не ошибка, но на будущее, это как будто несколько более экономно, не хранить избыточную информацию (особенно учитывая, что в реальных примерах она может весить гигабайты).


def __init__(self, seq):
super().__init__(seq)
self.seq = seq
self.check_alphabet()
def molecular_weight(self) -> float:
"""
Function calculates molecular weight of the amino acid chain
Input Parameters:
each letter refers to one-letter coded proteinogenic amino acids
Returns:
(float) Molecular weight of tge given amino acid chain in Da
"""
m = 0
for acid in str(self.seq):
m += self.masses[acid]
return m

class NucleicAcidSequence(BiologicalSequence):

Choose a reason for hiding this comment

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

Здесь есть некоторая проблема - gc_content и complement определяются здесь, но DNASequence и RNASequence наследуются от абстрактного класса. Получается, что этих методов у них нет :( . Так делать не стоило, следовало наследоваться от NucleicAcidSequence и тогда все было бы хорошо.

def __init__(self, seq):
super().__init__(seq)
self.check_alphabet()

def gc_content(self) -> float:
"""
Function counts GC-content in sequence, and returns result in %
"""
n = 0
for nucl in self.seq:
if nucl == 'c' or nucl == 'g' or nucl == 'C' or nucl == 'G':
n += 1
Comment on lines +67 to +69

Choose a reason for hiding this comment

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

Можно для упрощения кода использовать .count().

Copy link
Member Author

Choose a reason for hiding this comment

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

Так зато все считается за один прогон по строке, а не два. Думаю это не так существенно, один вариант чуть эффективнее, другой - чуть читаемее. Кому что больше нравится.

return 100 * n / len(self.seq)
Comment on lines +67 to +70

Choose a reason for hiding this comment

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

Вместо цикла можно использовать методы строк, типа

Suggested change
for nucl in self.seq:
if nucl == 'c' or nucl == 'g' or nucl == 'C' or nucl == 'G':
n += 1
return 100 * n / len(self.seq)
n = self.seq.upper().count('G') + self.seq.upper().count('C')
return 100 * n / len(self.seq)

либо использовать regex :)


def complement(self):
"""
Function return complement sequence
"""
complementary_dna = self.seq.maketrans(self.dict_comp)

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.

Очень интересный метод, я раньше не встречал, это здорово!

res = self.seq.translate(complementary_dna)
Comment on lines +76 to +77
Copy link

Choose a reason for hiding this comment

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

Интересно реализовано

return NucleicAcidSequence(res)

class DNASequence(BiologicalSequence):

Choose a reason for hiding this comment

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

Класс DNASequence должен наследоваться от класса NucleicAcidSequence, а не от абстрактного класса

Copy link

Choose a reason for hiding this comment

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

По идее не должно работать, если бы абстрактный класс содержал абстрактные методы. И получается, что проверка алфавита отсутствует

alphabet = {'A', 'g', 't', 'G', 'T', 'a', 'c', 'C'}
dict_trans = {'A': 'A', 'C': 'C', 'T': 'U', 'G': 'G', 'a': 'a', 'c': 'c', 't': 'u', 'g': 'g'}
dict_comp = {'A': 'T', 'C': 'G', 'T': 'A', 'G': 'C', 'a': 't', 'c': 'g', 't': 'a', 'g': 'c'}
Comment on lines +81 to +83

Choose a reason for hiding this comment

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

Снова заметка о том, что переменная alphabet, строго говоря, не нужна. В RNASequence та же ситуация будет. Ну и на самом деле вещи вроде 'A': 'A' можно было не делать, просто заменить T на U и t на u.

def __init__(self, seq):
super().__init__(seq)


def transcribe(self):
"""
Function return transcript of DNA sequence
"""
transcribe = self.seq.maketrans(self.dict_trans)
res = self.seq.translate(transcribe)
Comment on lines +88 to +93

Choose a reason for hiding this comment

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

Усложнили немного)
Тут можно обойтись просто заменой T на U, вместо использования отдельного алфавита под транскрипцию

Copy link
Member Author

Choose a reason for hiding this comment

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

Но в целом так тоже окей)


return RNASequence(res)

Choose a reason for hiding this comment

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

Круто, что предусмотрели, что возвращается объект другого класса!



class RNASequence(BiologicalSequence):

Choose a reason for hiding this comment

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

То же самое наследование от NucleicAcidSequence

Copy link

Choose a reason for hiding this comment

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

То же, что для DNASequence

alphabet = {'U', 'A', 'g', 'G', 'a', 'c', 'C', 'u'}
dict_comp = {'A': 'U', 'C': 'G', 'U': 'A', 'G': 'C', 'a': 'u', 'c': 'g', 'u': 'a', 'g': 'c'}
Comment on lines +99 to +100
Copy link

Choose a reason for hiding this comment

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

Хорошо, что предусмотрены регистры во всех алфавитах

def __init__(self, seq):
super().__init__(seq)


#фильтратор
def filter_gc(records, gc_bounds_both_side=(0, 100)) -> list:
"""
This function selects sequences with the GC content of your interest
:parameters:
records: records from fastq parced by SeqIO
gc_bound: interval for the of acceptable GC content, in %
:return:(dict) new dictionary consists of selected sequences

Choose a reason for hiding this comment

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

Следует поправить тип возвращаемого объекта.

"""
Comment on lines +106 to +113
Copy link

Choose a reason for hiding this comment

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

Видимо, забыли исправить докстрингу. Функция возвращает список, а не словарь.

new_records = []
for record in records:
if (gc_bounds_both_side[1]/100 >= SeqUtils.gc_fraction(record.seq) >= gc_bounds_both_side[0]/100):
new_records.append(record)

return new_records


def filter_length(records, length_bounds_both_side=(0, 2 ** 32)) -> list:
"""
This function selects sequences with the length of your interest
:parameters:
records: records from fastq parced by SeqIO
length_bound: interval for the of acceptable sequense length in number of nucleotide
:return:(dict) new dictionary consists of selected sequences

Choose a reason for hiding this comment

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

Аналогично.

"""
new_records = []
for record in records:
if (length_bounds_both_side[1] >= len(record.seq) >= length_bounds_both_side[0]):
new_records.append(record)
print(new_records)
Copy link

Choose a reason for hiding this comment

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

Suggested change
print(new_records)

Choose a reason for hiding this comment

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

Не очень понятно, почему здесь принт... Видимо, забыли убрать?

return new_records



def filter_quality(records, quality_threshold=0) -> list:
"""
This function selects FASTQ sequences with appropriate average nucleotide read quality
parameters:
seqs: dictionary of FASTQ sequences {name: (sequence, quality)}
quality_treshold: threshold value for average quality per nucleotide (phred33 scale)
:return:(dict) recordes for selected sequences

Choose a reason for hiding this comment

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

И еще раз.

"""
new_records = []
for record in records:
if (sum(record.letter_annotations["phred_quality"])/len(record.seq) >= quality_threshold):
new_records.append(record)
print(new_records)
Copy link

Choose a reason for hiding this comment

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

Suggested change
print(new_records)

Choose a reason for hiding this comment

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

И снова принт.

return new_records


def fastq_filtration(input_fastq, gc_bounds=(0, 100), length_bounds=(0, 2 ** 32), quality_treshold=0, output_fastq=''):
"""
This function provides you the opportunity to filter the FASTQ file to select sequences
that fit requirements on 5 parameters: input and output(optional) files, length, GC composition,
and quality of the reed
:parameters
input_fastq: path to fastq file
gc_bounds: (tuple) interval for the of acceptable GC content, in %
length_bounds: (tuple) interval for the of acceptable sequense length in number of nucleotide
quality_treshold: (float) threshold value for average quality per nucleotide (phred33 scale)
output_fastq = name of output file, ./fastq_filtrator_resuls/output_fastq, if it is not defined,
it will be the same of the input file

"""
if not os.path.isdir('fastq_filtrator_resuls'):
os.mkdir('fastq_filtrator_resuls')
Comment on lines +169 to +170

Choose a reason for hiding this comment

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

Папка для результатов - отличное решение)

if output_fastq == '':
output_fastq = os.path.join('fastq_filtrator_resuls', os.path.basename(input_fastq))
else:
output_fastq = os.path.join('fastq_filtrator_resuls', output_fastq + ".fasta")
if type(gc_bounds) == float or type(gc_bounds) == int:
gc_bounds_both_side = (0, gc_bounds)
else:
gc_bounds_both_side = gc_bounds
if type(length_bounds) == int:
length_bounds_both_side = (0, length_bounds)
else:
length_bounds_both_side = length_bounds
records = list(SeqIO.parse(input_fastq, "fastq"))

Choose a reason for hiding this comment

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

Думаю, так тоже рабочий вариант) Или же еще можно итерироваться по объектам SeqIO, и например, скормить только последовательность в виде seq_record.seq функцию filter_lenghth, кроме того метод .letter_annotations сам по себе спокойно работает с классом SeqRecord)

filtered_records_l = filter_length(records, length_bounds_both_side)
filtered_records_gc = filter_gc(filtered_records_l, gc_bounds_both_side)
filtered_records_q = filter_quality(filtered_records_gc, quality_treshold)
SeqIO.write((record for record in filtered_records_q), output_fastq, "fastq")
return
Comment on lines +187 to +188

Choose a reason for hiding this comment

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

Работа хорошая, мне понравилась) что то взяла себе на заметку. Удачи!


#fastq_filtration('example_fastq.fastq', length_bounds= (10, 30), gc_bounds=(20, 30), output_fastq='')