From 0ebdc003cfe46fbb048f2672c33c5567b792d2b1 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Sat, 7 Oct 2023 17:57:12 +0300 Subject: [PATCH 01/27] Create bioinf_tools.py --- bioinf_tools.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 bioinf_tools.py diff --git a/bioinf_tools.py b/bioinf_tools.py new file mode 100644 index 0000000..e69de29 From bb282ed8e787907019d1015e4717f56265c4448f Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Sat, 7 Oct 2023 18:08:17 +0300 Subject: [PATCH 02/27] Create bioinf_modules with dna_rna_tools.py, protein_tools.py, dictionaries.py --- bioinf_modules/dictionaries.py | 106 +++++++++++ bioinf_modules/dna_rna_tools.py | 71 ++++++++ bioinf_modules/protein_tools.py | 308 ++++++++++++++++++++++++++++++++ 3 files changed, 485 insertions(+) create mode 100644 bioinf_modules/dictionaries.py create mode 100644 bioinf_modules/dna_rna_tools.py create mode 100644 bioinf_modules/protein_tools.py diff --git a/bioinf_modules/dictionaries.py b/bioinf_modules/dictionaries.py new file mode 100644 index 0000000..f4a1ada --- /dev/null +++ b/bioinf_modules/dictionaries.py @@ -0,0 +1,106 @@ +amino_acids = { + "A": "Ala", + "C": "Cys", + "D": "Asp", + "E": "Glu", + "F": "Phe", + "G": "Gly", + "H": "His", + "I": "Ile", + "K": "Lys", + "L": "Leu", + "M": "Met", + "N": "Asn", + "P": "Pro", + "Q": "Gln", + "R": "Arg", + "S": "Ser", + "T": "Thr", + "V": "Val", + "W": "Trp", + "Y": "Tyr", + "a": "ala", + "c": "cys", + "d": "asp", + "e": "glu", + "f": "phe", + "g": "gly", + "h": "his", + "i": "ile", + "k": "lys", + "l": "leu", + "m": "met", + "n": "asn", + "p": "pro", + "q": "gln", + "r": "arg", + "s": "ser", + "t": "thr", + "v": "val", + "w": "trp", + "y": "tyr", +} +translation_rule = { + "F": "UUU", + "f": "uuu", + "L": "CUG", + "l": "cug", + "I": "AUU", + "i": "auu", + "M": "AUG", + "m": "aug", + "V": "GUG", + "v": "gug", + "P": "CCG", + "p": "ccg", + "T": "ACC", + "t": "acc", + "A": "GCG", + "a": "gcg", + "Y": "UAU", + "y": "uau", + "H": "CAU", + "h": "cau", + "Q": "CAG", + "q": "cag", + "N": "AAC", + "n": "aac", + "K": "AAA", + "k": "aaa", + "D": "GAU", + "d": "gau", + "E": "GAA", + "e": "gaa", + "C": "UGC", + "c": "ugc", + "W": "UGG", + "w": "ugg", + "R": "CGU", + "r": "cgu", + "S": "AGC", + "s": "agc", + "G": "GGC", + "g": "ggc", +} +amino_acid_weights = { + "A": 89.09, + "C": 121.16, + "D": 133.10, + "E": 147.13, + "F": 165.19, + "G": 75.07, + "H": 155.16, + "I": 131.17, + "K": 146.19, + "L": 131.17, + "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, +} diff --git a/bioinf_modules/dna_rna_tools.py b/bioinf_modules/dna_rna_tools.py new file mode 100644 index 0000000..14e148a --- /dev/null +++ b/bioinf_modules/dna_rna_tools.py @@ -0,0 +1,71 @@ +complements_rule = {"A": "T", "T": "A", "G": "C", "C": "G", "a": "t", "t": "a", + "g": "c", "c": "g", "U": "A", "u": "a"} + + +def check_user_input(sequence): + if not all(i in "".join(complements_rule.keys()) for i in sequence): + raise ValueError("Invalid sequence given") + if "T" in sequence.upper() and "U" in sequence.upper(): + raise ValueError("New nucleic acid discovered, go get Nobel") + + +def isrna(sequence): + return ("U" in sequence or "u" in sequence) + + +def transcribe(sequence): + if isrna(sequence): + transcript = sequence.replace("U", "T").replace("u", "t") + else: + transcript = sequence.replace("T", "U").replace("t", "u") + return transcript + + +def reverse(sequence): + return sequence[::-1] + + +def complement(sequence): + if isrna(sequence): + complements_rule["A"] = "U" + complements_rule["a"] = "u" + complement_strand = "" + for letter in sequence: + complement_strand += complements_rule[letter] + complements_rule["A"] = "T" + complements_rule["a"] = "t" + return complement_strand + + +def reverse_complement(sequence): + if isrna(sequence): + complements_rule["A"] = "U" + complements_rule["a"] = "u" + reverse_complement_strand = "" + complement_strand = sequence[::-1] + for letter in complement_strand: + reverse_complement_strand += complements_rule[letter] + complements_rule["A"] = "T" + complements_rule["a"] = "t" + return reverse_complement_strand + + +procedures_functions = {"transcribe": transcribe, + "reverse": reverse, + "complement": complement, + "reverse_complement": reverse_complement} + + +def run_dna_rna_tools(*args): + procedure = args[-1] + if procedure not in procedures_functions.keys(): + raise ValueError("Wrong procedure") + sequences = list(args[:-1]) + processed_sequences = [] + for sequence in sequences: + check_user_input(sequence) + processed_sequences.append(procedures_functions[procedure](sequence)) + if len(processed_sequences) == 1: + return processed_sequences[0] + else: + return processed_sequences diff --git a/bioinf_modules/protein_tools.py b/bioinf_modules/protein_tools.py new file mode 100644 index 0000000..0de1815 --- /dev/null +++ b/bioinf_modules/protein_tools.py @@ -0,0 +1,308 @@ +import dictionaries + + +def three_one_letter_code(sequences: str) -> list: + """ + Reverse the protein sequences from one-letter to three-letter format and vice-versa + + Case 1: get three-letter sequence\n + Use one-letter amino-acids sequences of any letter case + + Case 2: get one-letter sequence\n + Use three-letter amino-acid separated by "-" sequences. + Please note that sequences without "-" are parsed as one-letter code sequences\n + Example: for sequence "Ala" function will return "Ala-leu-ala" + + Arguments: + - sequences (tuple[str] or list[str]): protein sequences to convert\n + Example: ["WAG", "MkqRe", "msrlk", "Met-Ala-Gly", "Met-arg-asn-Trp-Ala-Gly", "arg-asn-trp"] + + Return: + - list: one-letter/three-letter protein sequences\n + Example: ["Met-Ala-Gly", "Met-arg-asn-Trp-Ala-Gly", "arg-asn-trp", "WAG", "MkqRe", "rlk"] + """ + inversed_sequences = [] + for sequence in sequences: + inversed_sequence = "" + if "-" not in sequence: + for letter in sequence: + inversed_sequence += dictionaries.amino_acids[letter] + "-" + inversed_sequence = inversed_sequence[:-1] + inversed_sequences.append(inversed_sequence) + else: + aa_splitted = sequence.split("-") + for aa in aa_splitted: + inversed_sequence += list(dictionaries.amino_acids.keys())[ + list(dictionaries.amino_acids.values()).index(aa) + ] + inversed_sequences.append(inversed_sequence) + return inversed_sequences + + +def define_molecular_weight(sequences: str) -> dict: + """ + Define molecular weight of the protein sequences + + Use one-letter amino-acids sequences of any letter case + The molecular weight is: + - a sum of masses of each atom constituting a molecule + - expressed in units called daltons (Da) + - rounded to hundredths + + Arguments: + - sequences (tuple[str] or list[str]): protein sequences to convert + + Return: + - dictionary: protein sequences as keys and molecular masses as values\n + Example: {"WAG": 332.39, "MkqRe": 690.88, "msrlk": 633.86} + """ + sequences_weights = {} + for sequence in sequences: + sequence_weight = 0 + for letter in sequence: + sequence_weight += dictionaries.amino_acid_weights[letter.upper()] + sequence_weight -= (len(sequence) - 1) * 18 # deduct water from peptide bond + sequences_weights[sequence] = round(sequence_weight, 2) + return sequences_weights + + +def search_for_motifs( + sequences: (tuple[str] or list[str]), motif: str, overlapping: bool +) -> dict: + """ + Search for motifs - conserved amino acids residues in protein sequence + + Search for one motif at a time\n + Search is letter case sensitive\n + Use one-letter aminoacids code for desired sequences and motifs\n + Positions of AA in sequences are counted from 0\n + By default, overlapping matches are counted + + Arguments: + - sequences (tuple[str] or list[str]): sequences to check for given motif within\n + Example: sequences = ["AMGAGW", "GAWSGRAGA"] + - motif (str]: desired motif to check presense in every given sequence\n + Example: motif = "GA" + - overlapping (bool): count (True) or skip (False) overlapping matches. (Optional)\n + Example: overlapping = False + Return: + - dictionary: sequences (str] as keys , starting positions for presented motif (list) as values\n + Example: {"AMGAGW": [2], "GAWSGRAGA": [0, 7]} + """ + new_line = "\n" + all_positions = {} + for sequence in sequences: + start = 0 + positions = [] + print(f"Sequence: {sequence}") + print(f"Motif: {motif}") + if motif in sequence: + while True: + start = sequence.find(motif, start) + if start == -1: + break + positions.append(start) + if overlapping: + start += 1 + else: + start += len(motif) + print_pos = ", ".join(str(x) for x in positions) + print_pos = f"{print_pos}{new_line}" + print( + f"Motif is present in protein sequence starting at positions: {print_pos}" + ) + else: + print(f"Motif is not present in protein sequence{new_line}") + all_positions[sequence] = positions + return all_positions + + +def search_for_alt_frames(sequences: str, alt_start_aa: str) -> dict: + """ + Search for alternative frames in a protein sequences + + Search is not letter case sensitive\n + Without an alt_start_aa argument search for frames that start with methionine ("M") + To search frames with alternative start codon add alt_start_aa argument\n + In alt_start_aa argument use one-letter code + + The function ignores the last three amino acids in sequences + + Arguments: + - sequences (tuple[str] or list[str]): sequences to check + - alt_start_aa (str]: the name of an amino acid that is encoded by alternative start AA (Optional)\n + Example: alt_start_aa = "I" + + Return: + - dictionary: the number of a sequence and a collection of alternative frames + """ + alternative_frames = {} + num_position = 0 + for sequence in sequences: + alternative_frames[sequence] = [] + for amino_acid in sequence[1:-3]: + alt_frame = "" + num_position += 1 + if amino_acid == alt_start_aa or amino_acid == alt_start_aa.swapcase(): + alt_frame += sequence[num_position:] + alternative_frames[sequence].append(alt_frame) + num_position = 0 + return alternative_frames + + +def convert_to_nucl_acids(sequences: list, nucl_acids: str) -> dict: + """ + Convert protein sequences to RNA or DNA sequences. + + Use the most frequent codons in human. The source - https://www.genscript.com/tools/codon-frequency-table\n + All nucleic acids (DNA and RNA) are showed in 5"-3" direction + + Arguments: + - sequences (tuple[str] or list[str]): sequences to convert + - nucl_acids (str]: the nucleic acid that is prefered\n + Example: nucl_acids = "RNA" - convert to RNA\n + nucl_acids = "DNA" - convert to DNA\n + nucl_acids = "both" - convert to RNA and DNA + Return: + - dictionary: nucleic acids (str) as keys, collection of sequences (list) as values + """ + rule_of_translation = sequences[0].maketrans(dictionaries.translation_rule) + rule_of_transcription = sequences[0].maketrans("AaUuCcGg", "TtAaGgCc") + nucl_acid_seqs = {"RNA": [], "DNA": []} + for sequence in sequences: + rna_seq = sequence.translate(rule_of_translation) + dna_seq = rna_seq.translate(rule_of_transcription) + if nucl_acids == "RNA": + nucl_acid_seqs["RNA"].append(rna_seq) + if sequence == sequences[-1]: + del nucl_acid_seqs["DNA"] + if nucl_acids == "DNA": + nucl_acid_seqs["DNA"].append(dna_seq) + if sequence == sequences[-1]: + del nucl_acid_seqs["RNA"] + if nucl_acids == "both": + nucl_acid_seqs["RNA"].append(rna_seq) + nucl_acid_seqs["DNA"].append(dna_seq) + return nucl_acid_seqs + + +procedures_to_functions = { + "search_for_motifs": search_for_motifs, + "search_for_alt_frames": search_for_alt_frames, + "convert_to_nucl_acids": convert_to_nucl_acids, + "three_one_letter_code": three_one_letter_code, + "define_molecular_weight": define_molecular_weight, +} + + +def check_and_parse_user_input( + sequences: list[str] or tuple[str], **kwargs +) -> dict and str: + """ + Check if user input can be correctly processed\n + Parse sequences and arguments for desired procedure + + Arguments: + - sequences (list[str] or tuple[str]): sequences to process + - **kwargs - needed arguments for completion of desired procedure + + Return: + - string: procedure name + - dictionary: a collection of procedure arguments and their values + """ + if len(sequences) == 0: + raise ValueError("No sequences provided") + procedure = kwargs["procedure"] + if procedure not in procedures_to_functions.keys(): + raise ValueError("Wrong procedure") + allowed_inputs = set(dictionaries.amino_acids.keys()).union( + set(dictionaries.amino_acids.values()) + ) + allowed_inputs.add("-") + if procedure != "three_one_letter_code": + allowed_inputs -= set(dictionaries.amino_acids.values()) + for sequence in sequences: + allowed_inputs_seq = allowed_inputs.copy() + if procedure == "three_one_letter_code" and "-" in sequence: + allowed_inputs_seq -= set(dictionaries.amino_acids.keys()) + if not all( + aminoacids in allowed_inputs_seq for aminoacids in sequence.split("-") + ): + raise ValueError("Invalid sequence given") + else: + allowed_inputs_seq.remove("-") + allowed_inputs_seq -= set(dictionaries.amino_acids.values()) + if not all(aminoacids in allowed_inputs_seq for aminoacids in sequence): + raise ValueError("Invalid sequence given") + procedure_arguments = {} + if procedure == "search_for_motifs": + if "motif" not in kwargs.keys(): + raise ValueError("Please provide desired motif") + procedure_arguments["motif"] = kwargs["motif"] + if "overlapping" not in kwargs.keys(): + procedure_arguments["overlapping"] = True + else: + procedure_arguments["overlapping"] = kwargs["overlapping"] + elif procedure == "search_for_alt_frames": + if "alt_start_aa" not in kwargs.keys(): + procedure_arguments["alt_start_aa"] = "M" + else: + if len(kwargs["alt_start_aa"]) > 1: + raise ValueError("Invalid alternative start AA") + procedure_arguments["alt_start_aa"] = kwargs["alt_start_aa"] + elif procedure == "convert_to_nucl_acids": + if "nucl_acids" not in kwargs.keys(): + raise ValueError("Please provide desired type of nucl_acids") + if kwargs["nucl_acids"] not in {"DNA", "RNA", "both"}: + raise ValueError("Invalid nucl_acids argument") + procedure_arguments["nucl_acids"] = kwargs["nucl_acids"] + procedure_arguments["sequences"] = sequences + return procedure_arguments, procedure + + +def run_protein_tools(sequences: list[str] or tuple[str], **kwargs: str): + """ + Main function to process protein sequence by one of the developed tools.\n + Run one procedure at a time: + - Search for conserved amino acids residues in protein sequence + - Search for alternative frames in a protein sequences + - Convert protein sequences to RNA or DNA sequences + - Reverse the protein sequences from one-letter to three-letter format and vice-versa + - Define molecular weight of the protein sequences + + All functions except *search_for_alt_frames* are letter case sensitive\n + Provide protein sequence in one letter code.\n + You can obtain one letter code from three letter code with *three_one_letter_code*\n + If more information needed please see README or desired docstring + + Arguments: + - sequences (list[str] or tuple[str]): sequences to process + - procedure (str]: desired procedure: + - "search_for_motifs" + - "search_for_alt_frames" + - "convert_to_nucl_acids" + - "three_one_letter_code" + - "define_molecular_weight" + + For "search_for_motif" procedure provide: + - motif (str]: desired motif to check presense in every given sequence\n + Example: motif = "GA" + - overlapping (bool): count (True) or skip (False) overlapping matches. (Optional)\n + Example: overlapping = False + + For "search_for_alt_frames" procedure provide: + - alt_start_aa (str]: the name of an amino acid that is encoded by alternative start codon (Optional)\n + Example: alt_start_aa = "I" + + For "convert_to_nucl_acids" procedure provide: + - nucl_acids (str]: the nucleic acid to convert to\n + Example: nucl_acids = "RNA"\n + nucl_acids = "DNA"\n + nucl_acids = "both" + + Return: + - dict: Dictionary with processed sequences. Depends on desired tool\n + Please see Readme or desired docstring + """ + procedure_arguments, procedure = check_and_parse_user_input(sequences, **kwargs) + return procedures_to_functions[procedure](**procedure_arguments) From a9106af12c6bcff1c26205623d02ae5be16679be Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Sat, 7 Oct 2023 20:00:29 +0300 Subject: [PATCH 03/27] Increase performance,add dockstrings to dna_rna_tools.py, modify dictionaries.py --- bioinf_modules/dictionaries.py | 8 +- bioinf_modules/dna_rna_tools.py | 169 ++++++++++++++++++++++++-------- 2 files changed, 132 insertions(+), 45 deletions(-) diff --git a/bioinf_modules/dictionaries.py b/bioinf_modules/dictionaries.py index f4a1ada..4ffb682 100644 --- a/bioinf_modules/dictionaries.py +++ b/bioinf_modules/dictionaries.py @@ -1,4 +1,4 @@ -amino_acids = { +AMINO_ACIDS = { "A": "Ala", "C": "Cys", "D": "Asp", @@ -40,7 +40,7 @@ "w": "trp", "y": "tyr", } -translation_rule = { +TRANSLATION_RULE = { "F": "UUU", "f": "uuu", "L": "CUG", @@ -82,7 +82,7 @@ "G": "GGC", "g": "ggc", } -amino_acid_weights = { +AMINO_ACID_WEIGHTS = { "A": 89.09, "C": 121.16, "D": 133.10, @@ -104,3 +104,5 @@ "W": 204.23, "Y": 181.19, } +NUCL_ACIDS_COMPLEMENT_RULE = {"A": "T", "T": "A", "G": "C", "C": "G", "a": "t", "t": "a", + "g": "c", "c": "g", "U": "A", "u": "a"} diff --git a/bioinf_modules/dna_rna_tools.py b/bioinf_modules/dna_rna_tools.py index 14e148a..5b17213 100644 --- a/bioinf_modules/dna_rna_tools.py +++ b/bioinf_modules/dna_rna_tools.py @@ -1,19 +1,55 @@ -complements_rule = {"A": "T", "T": "A", "G": "C", "C": "G", "a": "t", "t": "a", - "g": "c", "c": "g", "U": "A", "u": "a"} +import dictionaries -def check_user_input(sequence): - if not all(i in "".join(complements_rule.keys()) for i in sequence): - raise ValueError("Invalid sequence given") - if "T" in sequence.upper() and "U" in sequence.upper(): - raise ValueError("New nucleic acid discovered, go get Nobel") +def check_user_input(sequences: str or list[str] or tuple[str], procedure): + """ + Check if user input can be correctly processed.\n + Arguments: + - sequences (str list[str] or tuple[str]): sequences to process. -def isrna(sequence): - return ("U" in sequence or "u" in sequence) + Return: + - sequences (list): sequences to process. + """ + if procedure not in DNA_RNA_PROCEDURES_FUNCTIONS.keys(): + raise ValueError("Wrong procedure") + if isinstance(sequences, str): + sequences = sequences.split() + for sequence in sequences: + if not set(sequence.upper()).issubset( + dictionaries.NUCL_ACIDS_COMPLEMENT_RULE.keys() + ): + raise ValueError("Invalid sequence given") + if "T" in sequence.upper() and "U" in sequence.upper(): + raise ValueError("Neither DNA or RNA sequence given") + return sequences + + +def isrna(sequence: str) -> bool: + """ + Check if given sequence is RNA sequence.\n + + Arguments: + - sequence (str): sequence to check. + + Return: + - boolean True/False. + """ + return "U" in sequence or "u" in sequence def transcribe(sequence): + """ + Transcribe nucleic acid sequence.\n + Also works for reverse transcription.\n + Letter case sensitive.\n + + Arguments: + - sequence (str): sequence to transcribe + + Return: + - trancript (sequence): transcribed sequence + """ if isrna(sequence): transcript = sequence.replace("U", "T").replace("u", "t") else: @@ -21,51 +57,100 @@ def transcribe(sequence): return transcript -def reverse(sequence): +def reverse(sequence: str) -> str: + """ + Convert nucleic acid sequence into its reverse counterpart.\n + Letter case sensitive.\n + + Arguments: + - sequence (str): sequence to convert + + Return: + - str: reversed sequence + """ return sequence[::-1] -def complement(sequence): +def complement(sequence: str) -> str: + """ + Convert nucleic acid sequence into its complement counterpart.\n + Letter case sensitive.\n + + Arguments: + - sequence (str): sequence to convert + + Return: + - str: comlepement sequence + """ if isrna(sequence): - complements_rule["A"] = "U" - complements_rule["a"] = "u" + dictionaries.NUCL_ACIDS_COMPLEMENT_RULE["A"] = "U" complement_strand = "" for letter in sequence: - complement_strand += complements_rule[letter] - complements_rule["A"] = "T" - complements_rule["a"] = "t" + if letter.islower(): + complement_strand += dictionaries.NUCL_ACIDS_COMPLEMENT_RULE[letter].lower() + else: + complement_strand += dictionaries.NUCL_ACIDS_COMPLEMENT_RULE[letter] + dictionaries.NUCL_ACIDS_COMPLEMENT_RULE["A"] = "T" return complement_strand -def reverse_complement(sequence): - if isrna(sequence): - complements_rule["A"] = "U" - complements_rule["a"] = "u" - reverse_complement_strand = "" - complement_strand = sequence[::-1] - for letter in complement_strand: - reverse_complement_strand += complements_rule[letter] - complements_rule["A"] = "T" - complements_rule["a"] = "t" - return reverse_complement_strand +def reverse_complement(sequence: str) -> str: + """ + Convert nucleic acid sequence into its reverse-complement counterpart.\n + Letter case sensitive.\n + Arguments: + - sequence (str): sequence to convert -procedures_functions = {"transcribe": transcribe, - "reverse": reverse, - "complement": complement, - "reverse_complement": reverse_complement} + Return: + - reverse_complement_strand (str): reverse-comlepement sequence + """ + reverse_complement_strand = complement(sequence) + reverse_complement_strand = reverse(reverse_complement_strand) + return reverse_complement_strand -def run_dna_rna_tools(*args): - procedure = args[-1] - if procedure not in procedures_functions.keys(): - raise ValueError("Wrong procedure") - sequences = list(args[:-1]) - processed_sequences = [] +DNA_RNA_PROCEDURES_FUNCTIONS = { + "transcribe": transcribe, + "reverse": reverse, + "complement": complement, + "reverse_complement": reverse_complement, +} + + +def run_dna_rna_tools( + sequences: str or list[str] or tuple[str], procedure: str +) -> str or dict: + """ + Process nucleic acid sequence by one of the developed tools.\n + Run one procedure at a time: + - Transcribe nucleic acid sequence. + - Convert nucleic acid sequence into its reverse counterpart. + - Convert nucleic acid sequence into its complement counterpart. + - Convert nucleic acid sequence into its reverse-complement counterpart. + + All functions are letter case sensitive.\n + If only one sequence provided - *sequences* can be string.\n + If more - please provide *sequences* as list or tuple.\n + If more information needed please see README or desired docstring. + + Arguments: + - sequences (str or list[str] or tuple[str]): sequences to process. + - procedure (str]: desired procedure: + - "transcribe" + - "reverse" + - "complement" + - "reverse_complement" + + Return: + - processed_sequence (dict or str): Sequences, processed with desired tool. + """ + sequences = check_user_input(sequences, procedure) + processed_sequences = {} for sequence in sequences: - check_user_input(sequence) - processed_sequences.append(procedures_functions[procedure](sequence)) + processed_sequences[sequence] = DNA_RNA_PROCEDURES_FUNCTIONS[procedure]( + sequence + ) if len(processed_sequences) == 1: - return processed_sequences[0] - else: - return processed_sequences + return list(processed_sequences.values())[0] + return processed_sequences From eb189cb90b129211b587fc149f84c5e96c177fba Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Sun, 8 Oct 2023 11:30:31 +0300 Subject: [PATCH 04/27] Increase performance, delete lower case in all dictionaries, update ptorein_tools.py --- bioinf_modules/dictionaries.py | 45 +----------- bioinf_modules/protein_tools.py | 119 ++++++++++++++++++++------------ 2 files changed, 77 insertions(+), 87 deletions(-) diff --git a/bioinf_modules/dictionaries.py b/bioinf_modules/dictionaries.py index 4ffb682..c95f54a 100644 --- a/bioinf_modules/dictionaries.py +++ b/bioinf_modules/dictionaries.py @@ -18,69 +18,29 @@ "T": "Thr", "V": "Val", "W": "Trp", - "Y": "Tyr", - "a": "ala", - "c": "cys", - "d": "asp", - "e": "glu", - "f": "phe", - "g": "gly", - "h": "his", - "i": "ile", - "k": "lys", - "l": "leu", - "m": "met", - "n": "asn", - "p": "pro", - "q": "gln", - "r": "arg", - "s": "ser", - "t": "thr", - "v": "val", - "w": "trp", - "y": "tyr", + "Y": "Tyr" } TRANSLATION_RULE = { "F": "UUU", - "f": "uuu", "L": "CUG", - "l": "cug", "I": "AUU", - "i": "auu", "M": "AUG", - "m": "aug", "V": "GUG", - "v": "gug", "P": "CCG", - "p": "ccg", "T": "ACC", - "t": "acc", "A": "GCG", - "a": "gcg", "Y": "UAU", - "y": "uau", "H": "CAU", - "h": "cau", "Q": "CAG", - "q": "cag", "N": "AAC", - "n": "aac", "K": "AAA", - "k": "aaa", "D": "GAU", - "d": "gau", "E": "GAA", - "e": "gaa", "C": "UGC", - "c": "ugc", "W": "UGG", - "w": "ugg", "R": "CGU", - "r": "cgu", "S": "AGC", - "s": "agc", "G": "GGC", - "g": "ggc", } AMINO_ACID_WEIGHTS = { "A": 89.09, @@ -104,5 +64,4 @@ "W": 204.23, "Y": 181.19, } -NUCL_ACIDS_COMPLEMENT_RULE = {"A": "T", "T": "A", "G": "C", "C": "G", "a": "t", "t": "a", - "g": "c", "c": "g", "U": "A", "u": "a"} +NUCL_ACIDS_COMPLEMENT_RULE = {"A": "T", "T": "A", "G": "C", "C": "G", "U": "A"} diff --git a/bioinf_modules/protein_tools.py b/bioinf_modules/protein_tools.py index 0de1815..995dc2d 100644 --- a/bioinf_modules/protein_tools.py +++ b/bioinf_modules/protein_tools.py @@ -1,7 +1,7 @@ import dictionaries -def three_one_letter_code(sequences: str) -> list: +def three_one_letter_code(sequences: (tuple[str] or list[str])) -> list: """ Reverse the protein sequences from one-letter to three-letter format and vice-versa @@ -23,23 +23,35 @@ def three_one_letter_code(sequences: str) -> list: """ inversed_sequences = [] for sequence in sequences: - inversed_sequence = "" + inversed_sequence = [] if "-" not in sequence: for letter in sequence: - inversed_sequence += dictionaries.amino_acids[letter] + "-" - inversed_sequence = inversed_sequence[:-1] - inversed_sequences.append(inversed_sequence) + if letter.islower(): + inversed_sequence.append( + dictionaries.AMINO_ACIDS[letter.capitalize()].lower() + ) + else: + inversed_sequence.append(dictionaries.AMINO_ACIDS[letter]) + inversed_sequences.append("-".join(inversed_sequence)) else: aa_splitted = sequence.split("-") for aa in aa_splitted: - inversed_sequence += list(dictionaries.amino_acids.keys())[ - list(dictionaries.amino_acids.values()).index(aa) - ] - inversed_sequences.append(inversed_sequence) + aa_index = list(dictionaries.AMINO_ACIDS.values()).index( + aa.capitalize() + ) + if aa[0].islower(): + inversed_sequence.append( + list(dictionaries.AMINO_ACIDS.keys())[aa_index].lower() + ) + else: + inversed_sequence.append( + list(dictionaries.AMINO_ACIDS.keys())[aa_index] + ) + inversed_sequences.append("".join(inversed_sequence)) return inversed_sequences -def define_molecular_weight(sequences: str) -> dict: +def define_molecular_weight(sequences: (tuple[str] or list[str])) -> dict: """ Define molecular weight of the protein sequences @@ -60,7 +72,7 @@ def define_molecular_weight(sequences: str) -> dict: for sequence in sequences: sequence_weight = 0 for letter in sequence: - sequence_weight += dictionaries.amino_acid_weights[letter.upper()] + sequence_weight += dictionaries.AMINO_ACID_WEIGHTS[letter.upper()] sequence_weight -= (len(sequence) - 1) * 18 # deduct water from peptide bond sequences_weights[sequence] = round(sequence_weight, 2) return sequences_weights @@ -117,7 +129,9 @@ def search_for_motifs( return all_positions -def search_for_alt_frames(sequences: str, alt_start_aa: str) -> dict: +def search_for_alt_frames( + sequences: (tuple[str] or list[str]), alt_start_aa: str +) -> dict: """ Search for alternative frames in a protein sequences @@ -150,7 +164,9 @@ def search_for_alt_frames(sequences: str, alt_start_aa: str) -> dict: return alternative_frames -def convert_to_nucl_acids(sequences: list, nucl_acids: str) -> dict: +def convert_to_nucl_acids( + sequences: (tuple[str] or list[str]), nucl_acids: str +) -> dict: """ Convert protein sequences to RNA or DNA sequences. @@ -166,27 +182,35 @@ def convert_to_nucl_acids(sequences: list, nucl_acids: str) -> dict: Return: - dictionary: nucleic acids (str) as keys, collection of sequences (list) as values """ - rule_of_translation = sequences[0].maketrans(dictionaries.translation_rule) - rule_of_transcription = sequences[0].maketrans("AaUuCcGg", "TtAaGgCc") + rule_of_translation = str.maketrans(dictionaries.TRANSLATION_RULE) + # add lower case pairs, because only upper case pairs are stored in dictionaries + rule_of_translation.update( + str.maketrans( + dict( + (k.lower(), v.lower()) for k, v in dictionaries.TRANSLATION_RULE.items() + ) + ) + ) nucl_acid_seqs = {"RNA": [], "DNA": []} for sequence in sequences: rna_seq = sequence.translate(rule_of_translation) - dna_seq = rna_seq.translate(rule_of_transcription) if nucl_acids == "RNA": nucl_acid_seqs["RNA"].append(rna_seq) - if sequence == sequences[-1]: - del nucl_acid_seqs["DNA"] - if nucl_acids == "DNA": + elif nucl_acids == "DNA": + dna_seq = rna_seq.replace("U", "T").replace("u", "t") nucl_acid_seqs["DNA"].append(dna_seq) - if sequence == sequences[-1]: - del nucl_acid_seqs["RNA"] - if nucl_acids == "both": + elif nucl_acids == "both": + dna_seq = rna_seq.replace("U", "T").replace("u", "t") nucl_acid_seqs["RNA"].append(rna_seq) nucl_acid_seqs["DNA"].append(dna_seq) + if nucl_acids == "RNA": + del nucl_acid_seqs["DNA"] + if nucl_acids == "DNA": + del nucl_acid_seqs["RNA"] return nucl_acid_seqs -procedures_to_functions = { +PROTEINS_PROCEDURES_TO_FUNCTIONS = { "search_for_motifs": search_for_motifs, "search_for_alt_frames": search_for_alt_frames, "convert_to_nucl_acids": convert_to_nucl_acids, @@ -196,12 +220,12 @@ def convert_to_nucl_acids(sequences: list, nucl_acids: str) -> dict: def check_and_parse_user_input( - sequences: list[str] or tuple[str], **kwargs + sequences: (str, tuple[str] or list[str]), **kwargs ) -> dict and str: """ Check if user input can be correctly processed\n Parse sequences and arguments for desired procedure - + Arguments: - sequences (list[str] or tuple[str]): sequences to process - **kwargs - needed arguments for completion of desired procedure @@ -210,29 +234,34 @@ def check_and_parse_user_input( - string: procedure name - dictionary: a collection of procedure arguments and their values """ - if len(sequences) == 0: - raise ValueError("No sequences provided") + if isinstance(sequences, str): + sequences = sequences.split() + if "" in sequences or len(sequences) == 0: + raise ValueError("Empty sequence provided") procedure = kwargs["procedure"] - if procedure not in procedures_to_functions.keys(): + if procedure not in PROTEINS_PROCEDURES_TO_FUNCTIONS.keys(): raise ValueError("Wrong procedure") - allowed_inputs = set(dictionaries.amino_acids.keys()).union( - set(dictionaries.amino_acids.values()) + allowed_inputs = set(dictionaries.AMINO_ACIDS.keys()) + allowed_inputs = allowed_inputs.union( + set(k.lower() for k in dictionaries.AMINO_ACIDS.keys()) ) - allowed_inputs.add("-") - if procedure != "three_one_letter_code": - allowed_inputs -= set(dictionaries.amino_acids.values()) + if procedure == "three_one_letter_code": + allowed_inputs = allowed_inputs.union(set(dictionaries.AMINO_ACIDS.values())) + allowed_inputs = allowed_inputs.union( + set(v.lower() for v in dictionaries.AMINO_ACIDS.values()) + ) for sequence in sequences: allowed_inputs_seq = allowed_inputs.copy() if procedure == "three_one_letter_code" and "-" in sequence: - allowed_inputs_seq -= set(dictionaries.amino_acids.keys()) - if not all( - aminoacids in allowed_inputs_seq for aminoacids in sequence.split("-") - ): + allowed_inputs_seq -= set(dictionaries.AMINO_ACIDS.keys()) + allowed_inputs_seq -= set( + k.lower() for k in dictionaries.AMINO_ACIDS.keys() + ) + allowed_inputs_seq.union(set("-")) + if not set(sequence.split("-")).issubset(allowed_inputs_seq): raise ValueError("Invalid sequence given") else: - allowed_inputs_seq.remove("-") - allowed_inputs_seq -= set(dictionaries.amino_acids.values()) - if not all(aminoacids in allowed_inputs_seq for aminoacids in sequence): + if not set(sequence).issubset(allowed_inputs_seq): raise ValueError("Invalid sequence given") procedure_arguments = {} if procedure == "search_for_motifs": @@ -260,9 +289,9 @@ def check_and_parse_user_input( return procedure_arguments, procedure -def run_protein_tools(sequences: list[str] or tuple[str], **kwargs: str): +def run_protein_tools(sequences: (str, tuple[str] or list[str]), **kwargs: str) -> dict: """ - Main function to process protein sequence by one of the developed tools.\n + Process protein sequence by one of the developed tools.\n Run one procedure at a time: - Search for conserved amino acids residues in protein sequence - Search for alternative frames in a protein sequences @@ -271,12 +300,14 @@ def run_protein_tools(sequences: list[str] or tuple[str], **kwargs: str): - Define molecular weight of the protein sequences All functions except *search_for_alt_frames* are letter case sensitive\n + If only one sequence provided - *sequences* can be string.\n + If more - please provide *sequences* as list or tuple.\n Provide protein sequence in one letter code.\n You can obtain one letter code from three letter code with *three_one_letter_code*\n If more information needed please see README or desired docstring Arguments: - - sequences (list[str] or tuple[str]): sequences to process + - sequences (str, list[str] or tuple[str]): sequences to process - procedure (str]: desired procedure: - "search_for_motifs" - "search_for_alt_frames" @@ -305,4 +336,4 @@ def run_protein_tools(sequences: list[str] or tuple[str], **kwargs: str): Please see Readme or desired docstring """ procedure_arguments, procedure = check_and_parse_user_input(sequences, **kwargs) - return procedures_to_functions[procedure](**procedure_arguments) + return PROTEINS_PROCEDURES_TO_FUNCTIONS[procedure](**procedure_arguments) From 0acfa301fc20bf8b1a60103a058a51f7f475e013 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Sun, 8 Oct 2023 15:34:26 +0300 Subject: [PATCH 05/27] Create fastq_tools.py, functions and dockstrings within --- bioinf_modules/fastq_tools.py | 164 ++++++++++++++++++++++++++++++++++ 1 file changed, 164 insertions(+) create mode 100644 bioinf_modules/fastq_tools.py diff --git a/bioinf_modules/fastq_tools.py b/bioinf_modules/fastq_tools.py new file mode 100644 index 0000000..5adf85e --- /dev/null +++ b/bioinf_modules/fastq_tools.py @@ -0,0 +1,164 @@ +import dictionaries + + +def run_fastq_filter( + seqs: dict[str, tuple[str] | list[str]] = { + "@SRX079804:1:SRR292678:1:1101:21885:21885": ( + "ACAGCAACATAAACATGATGGGATGGCGTAAGCCCCCGAGATATCAGTTTACCCAGGATAAGAGATTAAATTATGAGCAACATTATTAA", + "FGGGFGGGFGGGFGDFGCEBB@CCDFDDFFFFBFFGFGEFDFFFF;D@DD>C@DDGGGDFGDGG?GFGFEGFGGEF@FDGGGFGFBGGD", + ), + "@SRX079804:1:SRR292678:1:1101:30161:30161": ( + "GAACGACAGCAGCTCCTGCATAACCGCGTCCTTCTTCTTTAGCGTTGTGCAAAGCATGTTTTGTATTACGGGCATCTCGAGCGAATC", + "DFFFEGDGGGGFGGEDCCDCEFFFFCCCCCB>CEBFGFBGGG?DE=:6@=>AD?D8DCEE:>EEABE5D@5:DDCA;EEE-DCD", + ), + }, + gc_bounds: (int | float | tuple[int | float] | list[int | float]) = (0, 100), + length_bounds: (tuple[int]) = (0, 2**32), + quality_threshold: (int) = 0, +) -> dict: + """ + Filter out fastq reads by several parameters: + - GC content + - length + - mean nucleotide quality + Please provide GC content bounds in percentages.\n + Please provide threshold quality in Phred-33 scale.\n + For GC content and length bounds: if one number is provided, bounds from 0 to number are considered.\n + + Arguments: + - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered + - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds + - length_bounds (int | tuple[int] | list[int]): read length thresholds + - quality_thresholds: read Phred-33 scaled quality thresholds + + Return: + - seqs_filtered (dict): similar dictionary as input, bad reads are filtered out + """ + seqs, gc_bounds, length_bounds, quality_threshold = check_user_input( + seqs, gc_bounds, length_bounds, quality_threshold + ) + return fastq_filter(seqs, gc_bounds, length_bounds, quality_threshold) + + +def check_user_input( + seqs: dict[str, tuple[str] | list[str]], + gc_bounds: (int | float | tuple[int | float] | list[int | float]), + length_bounds: (int | tuple[int] | list[int]), + quality_threshold: int, +): + """ + Check if user input can be correctly processed\n + + Arguments: + - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered + - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds + - length_bounds (int | tuple[int] | list[int]): read length thresholds + - quality_thresholds: read Phred-33 scaled quality thresholds + + Return: + - same arguments as input, checked for correctness + """ + if not isinstance(seqs, dict): + raise ValueError("Please provide sequences info with a dictionary") + for seq_name in seqs.keys(): + if not isinstance(seq_name, str): + raise ValueError("Invalid sequence name given") + if not set(seqs[seq_name][0].upper()).issubset({"A", "T", "G", "C"}): + raise ValueError("Invalid read sequence given") + if not set(seqs[seq_name][1]).issubset(dictionaries.QUALITY_SYMBOLS): + raise ValueError("Invalid quality sequence given") + return seqs, gc_bounds, length_bounds, quality_threshold + + +def fastq_filter( + seqs: dict[str, tuple[str] | list[str]], + gc_bounds: (int | float | tuple[int | float] | list[int | float]), + length_bounds: (int | tuple[int] | list[int]), + quality_threshold: (int | float), +) -> dict: + """ + Parse checked input and filter out bad reads.\n + + Arguments: + - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered + - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds + - length_bounds (int | tuple[int] | list[int]): read length thresholds + - quality_thresholds: read Phred-33 scaled quality thresholds + + Return: + - seqs_filtered (dict): similar dictionary as input, bad reads are filtered out + """ + seqs_filtered = {} + for seq_name in seqs.keys(): + seq = seqs[seq_name][0] + seq_qual = seqs[seq_name][1] + gc_result = is_gc_good(seq, gc_bounds) + len_result = is_len_good(seq, length_bounds) + qual_result = is_qual_good(seq_qual, quality_threshold) + if gc_result and len_result and qual_result: + seqs_filtered[seq_name] = seqs[seq_name] + return seqs_filtered + + +NEW_LINE = "\n" # needed for output in f-strings + + +def is_gc_good( + seq: str, + gc_bounds: (int | float | tuple[int | float] | list[int | float]), +) -> bool: + """ + Check GC content of a given read.\n + Please provide bounds in percentages.\n + If one number is provided, bounds from 0 to number are considered.\n + + Arguments: + - seq (str): read to check it's length + - gc_bounds (int | float | tuple[int] | list[int]): GC content thresholds, by which reads are filtered + + Return: + - condition (bool): True - GC content is within bounds, False - read is to be filtered + """ + if isinstance(gc_bounds, int) or isinstance(gc_bounds, float): + gc_bounds = (0, gc_bounds) + gc_content = seq.count("G") + seq.count("C") + gc_content = gc_content / len(seq) * 100 + print(f"Read: {seq}") + print(f"GC Content: {round(gc_content, 4)}") + return gc_bounds[0] <= gc_content <= gc_bounds[1] + + +def is_len_good(seq: str, length_bounds: (int | tuple[int] | list[int])) -> bool: + """ + Check length of a given read.\n + If one number is provided, bounds from 1 to number are considered.\n + + Arguments: + - seq (str): read to check it's length + - length_bounds (int | tuple[int] | list[int]): length thresholds, by which reads are filtered + + Return: + - condition (bool): True - length is within bounds, False - read is to be filtered + """ + if isinstance(length_bounds, int): + length_bounds = (1, length_bounds) + print(f"Read Length: {round(len(seq), 4)}") + return length_bounds[0] <= len(seq) <= length_bounds[1] + + +def is_qual_good(seq_qual: str, quality_threshold: int | float) -> bool: + """ + Check mean nucleotide sequencing quality for a given read.\n + Reads with mean quality less then threshold are filtered out.\n + Please provide Phred-33 Q-Scores as input.\n + + Arguments: + - seq_qual (str): quality sequence in Phred-33 scale for a read + - quality_threshold (int | float): threshold, by which reads are filtered + + Return: + - condition (bool): True - quality is above threshold, False - read is to be filtered + """ + mean_quality = sum(ord(symbol) - 33 for symbol in seq_qual) / len(seq_qual) + print(f"Mean Nucleotide Quality: {round(mean_quality, 4)}{NEW_LINE}") + return mean_quality > quality_threshold From 7671daa26b4990b99ea3c6808ad8747916f5b6bc Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Sun, 8 Oct 2023 15:40:00 +0300 Subject: [PATCH 06/27] Rename bioinf_modules to src --- {bioinf_modules => src}/dictionaries.py | 43 ++++++++++++++++++++++++ {bioinf_modules => src}/dna_rna_tools.py | 0 {bioinf_modules => src}/fastq_tools.py | 0 {bioinf_modules => src}/protein_tools.py | 0 4 files changed, 43 insertions(+) rename {bioinf_modules => src}/dictionaries.py (74%) rename {bioinf_modules => src}/dna_rna_tools.py (100%) rename {bioinf_modules => src}/fastq_tools.py (100%) rename {bioinf_modules => src}/protein_tools.py (100%) diff --git a/bioinf_modules/dictionaries.py b/src/dictionaries.py similarity index 74% rename from bioinf_modules/dictionaries.py rename to src/dictionaries.py index c95f54a..7efc1e3 100644 --- a/bioinf_modules/dictionaries.py +++ b/src/dictionaries.py @@ -65,3 +65,46 @@ "Y": 181.19, } NUCL_ACIDS_COMPLEMENT_RULE = {"A": "T", "T": "A", "G": "C", "C": "G", "U": "A"} +QUALITY_SYMBOLS = { + ".", + "2", + "5", + "+", + "-", + ",", + "8", + "4", + "B", + "&", + ")", + "#", + "7", + ">", + "0", + "!", + "=", + "C", + "D", + "1", + "%", + "A", + "9", + "?", + "F", + "@", + "G", + "$", + "H", + ":", + "/", + ";", + "I", + "<", + "3", + "(", + "E", + '"', + "*", + "6", + "'", +} diff --git a/bioinf_modules/dna_rna_tools.py b/src/dna_rna_tools.py similarity index 100% rename from bioinf_modules/dna_rna_tools.py rename to src/dna_rna_tools.py diff --git a/bioinf_modules/fastq_tools.py b/src/fastq_tools.py similarity index 100% rename from bioinf_modules/fastq_tools.py rename to src/fastq_tools.py diff --git a/bioinf_modules/protein_tools.py b/src/protein_tools.py similarity index 100% rename from bioinf_modules/protein_tools.py rename to src/protein_tools.py From cd6bdccad683e43235cbee51f813157be742c0f9 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Sun, 8 Oct 2023 18:52:04 +0300 Subject: [PATCH 07/27] Increase performance, update bioinf_tools.py --- bioinf_tools.py | 144 ++++++++++++++++++ src/__pycache__/dictionaries.cpython-311.pyc | Bin 0 -> 2027 bytes src/__pycache__/dna_rna_tools.cpython-311.pyc | Bin 0 -> 4454 bytes src/__pycache__/fastq_tools.cpython-311.pyc | Bin 0 -> 7341 bytes src/__pycache__/protein_tools.cpython-311.pyc | Bin 0 -> 15212 bytes src/dna_rna_tools.py | 38 ----- src/fastq_tools.py | 71 +++------ src/protein_tools.py | 50 ------ 8 files changed, 164 insertions(+), 139 deletions(-) create mode 100644 src/__pycache__/dictionaries.cpython-311.pyc create mode 100644 src/__pycache__/dna_rna_tools.cpython-311.pyc create mode 100644 src/__pycache__/fastq_tools.cpython-311.pyc create mode 100644 src/__pycache__/protein_tools.cpython-311.pyc diff --git a/bioinf_tools.py b/bioinf_tools.py index e69de29..1986edf 100644 --- a/bioinf_tools.py +++ b/bioinf_tools.py @@ -0,0 +1,144 @@ +from src import dna_rna_tools +from src import fastq_tools +from src import protein_tools + + +def run_dna_rna_tools( + sequences: str or list[str] or tuple[str], procedure: str +) -> str or dict: + """ + Process nucleic acid sequence by one of the developed tools.\n + Run one procedure at a time: + - Transcribe nucleic acid sequence. + - Convert nucleic acid sequence into its reverse counterpart. + - Convert nucleic acid sequence into its complement counterpart. + - Convert nucleic acid sequence into its reverse-complement counterpart. + + All functions are letter case sensitive.\n + If only one sequence provided - *sequences* can be string.\n + If more - please provide *sequences* as list or tuple.\n + If more information needed please see README or desired docstring. + + Arguments: + - sequences (str or list[str] or tuple[str]): sequences to process. + - procedure (str]: desired procedure: + - "transcribe" + - "reverse" + - "complement" + - "reverse_complement" + + Return: + - processed_sequence (dict or str): Sequences, processed with desired tool. + """ + sequences = dna_rna_tools.check_user_input(sequences, procedure) + processed_sequences = {} + for sequence in sequences: + processed_sequences[sequence] = dna_rna_tools.DNA_RNA_PROCEDURES_FUNCTIONS[ + procedure + ](sequence) + if len(processed_sequences) == 1: + return list(processed_sequences.values())[0] + return processed_sequences + + +def run_protein_tools(sequences: (str, tuple[str] or list[str]), **kwargs: str) -> dict: + """ + Process protein sequence by one of the developed tools.\n + Run one procedure at a time: + - Search for conserved amino acids residues in protein sequence + - Search for alternative frames in a protein sequences + - Convert protein sequences to RNA or DNA sequences + - Reverse the protein sequences from one-letter to three-letter format and vice-versa + - Define molecular weight of the protein sequences + + All functions except *search_for_alt_frames* are letter case sensitive\n + If only one sequence provided - *sequences* can be string.\n + If more - please provide *sequences* as list or tuple.\n + Provide protein sequence in one letter code.\n + You can obtain one letter code from three letter code with *three_one_letter_code*\n + If more information needed please see README or desired docstring + + Arguments: + - sequences (str, list[str] or tuple[str]): sequences to process + - procedure (str]: desired procedure: + - "search_for_motifs" + - "search_for_alt_frames" + - "convert_to_nucl_acids" + - "three_one_letter_code" + - "define_molecular_weight" + + For "search_for_motif" procedure provide: + - motif (str]: desired motif to check presense in every given sequence\n + Example: motif = "GA" + - overlapping (bool): count (True) or skip (False) overlapping matches. (Optional)\n + Example: overlapping = False + + For "search_for_alt_frames" procedure provide: + - alt_start_aa (str]: the name of an amino acid that is encoded by alternative start codon (Optional)\n + Example: alt_start_aa = "I" + + For "convert_to_nucl_acids" procedure provide: + - nucl_acids (str]: the nucleic acid to convert to\n + Example: nucl_acids = "RNA"\n + nucl_acids = "DNA"\n + nucl_acids = "both" + + Return: + - dict: Dictionary with processed sequences. Depends on desired tool\n + Please see Readme or desired docstring + """ + procedure_arguments, procedure = protein_tools.check_and_parse_user_input( + sequences, **kwargs + ) + return protein_tools.PROTEINS_PROCEDURES_TO_FUNCTIONS[procedure]( + **procedure_arguments + ) + + +def run_fastq_filter( + seqs: dict[str, tuple[str] | list[str]] = { + "@SRX079804:1:SRR292678:1:1101:21885:21885": ( + "ACAGCAACATAAACATGATGGGATGGCGTAAGCCCCCGAGATATCAGTTTACCCAGGATAAGAGATTAAATTATGAGCAACATTATTAA", + "FGGGFGGGFGGGFGDFGCEBB@CCDFDDFFFFBFFGFGEFDFFFF;D@DD>C@DDGGGDFGDGG?GFGFEGFGGEF@FDGGGFGFBGGD", + ), + "@SRX079804:1:SRR292678:1:1101:30161:30161": ( + "GAACGACAGCAGCTCCTGCATAACCGCGTCCTTCTTCTTTAGCGTTGTGCAAAGCATGTTTTGTATTACGGGCATCTCGAGCGAATC", + "DFFFEGDGGGGFGGEDCCDCEFFFFCCCCCB>CEBFGFBGGG?DE=:6@=>AD?D8DCEE:>EEABE5D@5:DDCA;EEE-DCD", + ), + }, + gc_bounds: (int | float | tuple[int | float] | list[int | float]) = (0, 100), + length_bounds: (tuple[int]) = (0, 2**32), + quality_threshold: (int) = 0, + verbose=True, +) -> dict: + """ + Filter out fastq reads by several parameters: + - GC content + - length + - mean nucleotide quality + Please provide GC content bounds in percentages.\n + Please provide threshold quality in Phred-33 scale.\n + For GC content and length bounds: if one number is provided, bounds from 0 to number are considered.\n + + Arguments: + - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered + - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds + - length_bounds (int | tuple[int] | list[int]): read length thresholds + - quality_thresholds: read Phred-33 scaled quality thresholds + Examples are provided as default values for each argument + + Return: + - seqs_filtered (dict): similar dictionary as input, bad reads are filtered out + """ + ( + seqs, + gc_bounds, + length_bounds, + quality_threshold, + verbose, + ) = fastq_tools.check_user_input( + seqs, gc_bounds, length_bounds, quality_threshold, verbose + ) + return fastq_tools.fastq_filter( + seqs, gc_bounds, length_bounds, quality_threshold, verbose + ) diff --git a/src/__pycache__/dictionaries.cpython-311.pyc b/src/__pycache__/dictionaries.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..4a24e0308a0c930fc1a68b13c469785f08ebe7e8 GIT binary patch literal 2027 zcmY+EOKcle6o%(`;?!{-b`xisrp-HT;x?_@rf-_mj^j?;I?ls(t40MAa$Lt!5=Y~- zQdm_&2*id33k2$}P+(KB>=GdaU4Yab)m?$ufsjV3U@tza9sgB=_J zqr4yNz5ID>uJjzEN>Yvk!oEqYz813XSu!G0J z7{|e5d>riI32=~4fEu3!6MPD6huBcm@n}8f@cPu!A#Tl&^uEuE*Gf8$1X0a2D+4>)>I&0rqna z9OQX0&kMZxDr*2^8txPZmmUUJ^mUDwAJIw^7YbUus3@iM)W)i!w9?bN4Mnr)H?|dJ zlx}X^Q*=$~i|Z9dru26e%_)7ZhIv_~XE!i@L+LqLm&5-G=FcmAzOthzFZ*mMTEu#* zu4qZtqrZ%CtE#A=^g_i_R8)F#-BS9^)eV@n4eVNRUC~WNw-nu0dZ|>x-AcGXN-R|} ztucIxu%C>j(^5#|`exd+;%=dok`*4LnnE>|LQ+%4OtF%gw#L25rtD(M&Zg{QrsKNR zkCHVgXalN$kKg(5+sw^mHT?6V?WyZ;B&(Cs9=(~YqVU%^zm=>$F@JqK_L!5^nAmsA z$?Cg*EIsY7tR<^cvaEv1?wo6wJRr5lkk7VJda#?n?WQFKt531JGd^W#OPNlQ61r;qWr}Bkd zs+e8KmzPVqjOzTiT)vgbnsdd1>dlwZId4Tdy)eI+%gkr;MQ?s&sg%lPi?_>#+w(IE zxx!c4OWad)7809VJBhVKsZp^SiTm~YtJTD8eeK?6WoxIA$j&4dcX!t7TjkS@%4TBj z*1582)v9%?wkqX@wU*%8+D@&$wQAKWjg#BEPUy;Jo$qZ_rmc1Cj~_Cv@h9GOTPGqr z&>ac9XxjIFx2?l|aV-e5o3&sFqx-}|p0|jHJ#Q6n^SoWW!}ExE)bst~ot}4z$2{*A zH#|Qe-sAZ}@m|jli68d7PrTpr0r5f4hs1|H9}yqb!jTs(9eCMs>1%s<&}vU~$NZo% zKWOYfLE(N-xE~bm2Zj4V;eJrK9~ABfh5JF_eo(j{6z&Iw`$6G;P`Doy{=cB;YY6A! zxfVrI!%u3sNew%QG4GKYc%*PQsbRN>yGaeZRoqVscas`+M8@5uhP_|hO={R(;%-vI z?iP2G8ukHkH>qJC6nB#v_91aMsbTktyGaduK-^7g*hAuOQo|k*N7Cqv7ChAm9x3|T q9-eER_JlS%_GRJg#y4HhCZG2$$V==XPSv$y;LGl$}V&*|DzW}AY4RKOCC72eU%1|h1t6TWAEXQ=8$Mzsx= zpzrAO2)dp$x>p$zRO<%-y#n(HW<6=l{E$&S4Q5nCs<+vC1e!c)XhXt9RI*93@jC2@je>=Q&u7g_ljXt+mQ-cLY7TP^aB+i*O1PY;s6T|tyPQlowUEzYs}&!y zF{*N1E=Z-S%wn}bq+=vfX(vvlGuMOf5LsHT$`c9%uXo>&Rp{NnR%00 zSdUF%QN08E9iK{54QHXK@1V*U_nW;in1PYZ%mcHdnwCgviA7iHayGY8G6~bi&?Emy*wkD%7!Rx(fL`CZ<a@|7s?= zBB@C}xd28~lB@D+b~!mN^Hm!qZ;U2q*3>()l=-fLSCVhM^GfCiL|m4Mm<1Ci6~ZS4 zDVrhisj{3?68SZ=gQvjAP*h|rRFDX)#lol1fbs#Vgps%NuhGl<(O37Puj6{mhVYW=Qi69e6so;-Advo|xf@RS@HPp&ahpo_;v3hX;)Cz-AjXE6D8{ zWvX=3dk_q5UE5AQOzB-0jIIlMu-^#w?+1tXg2M-as1dldAGo#`xTXh2jKIi0osJIQ zzYvfLa<*Ol!i~E6w7ysL&=n(eMGIUh`_F6*ZYLfl^w4=DbYAn*x1{Xe`JMDzNl&DV zL`u6z--A%E*7toqlr}f_!(>#kuD&nn^*wl8=(z4Y4EbbCyf~p`I#bG zMR}*-pe4AObNv^sZf3!-2{ut_0pHj9m*rkIU*N!tEK42hLrY*y49+}H3OMmu^*kRHQ zV=d=F0b-rUU!by$RNGkRSJYb$at8D=9n?Aw)_AK8Hn~%ks8RQ!_<6r^wtrncP*l+Yn<2 zUIpL8oLZZmr`}&Sjs7J}V$dIlPr*PwKwm_W zw|#T*$M4^N9~>9FzZ~h?j|}WZ21@+yB|VZfBFX*8)xF5or_QHydL(T`(wn|b-)E;` zEBhlF8+i^qqo4z1JK(SsY#%wB37X}!f2caGEP>6`azO$!=Cefo zzuT3ss&brB;i?&_>u|_`Q=k6_9?}DwCFkhT+uR&XXaSP+!@z5Fka|BjdN&G$x5nYO zv^D-9rLkw9bBrAz@N2w0C@^Xx#Y5=5MqO~wRRt$_7ft8pHXM-R(AaV0X-(nCAS9^` zDRzlPN7IZY8i|fPt&-LuecWIkq)sER#yc;t$FE)8O)XCR`J2FpmscP^qS+K=T;Cpo zK#d&ESM{g|hh`KonY4+lP!M-nA%Jwi;gU`urbkqA@;%6ITvh=Iixm?(n;XjjMzfM& zx^q^btEVm1`^MMgBQJrxH{nx$2?Sz9w=G66uIuHiBR>ro-sGoF-Fvwl?ExPzhk7cA zb-fMA1ErhpvcL1^sh^~D|2e~dj-~`}Gvz?vuX~Ny8-MIJVl(@(n|raFdTh~%E$V?q zEzk$=XLb%fFjMUGidgH$8->dAePM%9ZTK!8`lD34 zWtjCZzX^J3UOe19#BU2lHBB;2O=@PZ=PWi=k4GBZCKLcxM5^B`XHOK{a?ym2vXLan z%EG5y1Y&_HHne+9=Pn!EWsU7S;4W%|X`LH2xKU`9{rwu-Z`n7_T9=ijM^z&#Djuj} z-Vb$?rT=t{{kkBHFE#} literal 0 HcmV?d00001 diff --git a/src/__pycache__/fastq_tools.cpython-311.pyc b/src/__pycache__/fastq_tools.cpython-311.pyc new file mode 100644 index 0000000000000000000000000000000000000000..fe9d9859c3ce5fd9e359d90fd86b0684df40bc83 GIT binary patch literal 7341 zcmeHLTWs6b89tOmiMm*dE#I9uOkF3I?}?qXa2&gdlctU@O`5btGXyRzQnnVEa!4wT zU3v9_KKLOIPLU!C)&T03v?{s+-mr)4VOxu>K(Qw=AW(rouman|vL|7|fFe)(&!H}) z?6|?!0(%%m{p7j*=jS=+|NifskE^OY2(Fe-o)K1p2>lf&iWez?#G`L1gl-}p@l*_* zN4!m-1e;(NY)!?QNN|`_%Pxc>rc&N6(1H`^LaV}Cxh=_eY@N^Jd}y<|wl#;R1vUH_~086&{E1S;a;8L*V$-RIrBN3i`De*p1DZo zNveZ)cj4VSx1P#&Ft$|Emb6D4unfK2A7R=(EL>yrLii#Zjj#z>kl3iWm{8a-C$h5w z8;(no5LRMK?4lG83$iTm9uH}BRGLdH2%;hnkz@~xmowYJN5jgGWJT&`mBeC9z^@nB zBpZv$B>TdJ&=4EpWaSbo2^=r8N*s3+iN+K`0^Qi09=7BYY)4d7pm`(~=kU|eY04|< zHs?VXC&$2aQ33ndlKaR7?G`sH_E6Gn)`*S8X0k?;&8jxJs9b1p8hn=T85m&YFc%Ye zwouWmfmrj3pd=)*G=$uOz;eZbWjWcHjBYlJX9(0p;!<%+M3@YsV!7{hOyFbz&(P&4 zFW`wt2qH`tw2s8tx1!2C%dvQ3qH&RvmX=$_#mfLKe5ocY;*RH{mjzM6NHO97p;Je7 zYD}k2>eRV%6EYmd9*iarF?=ky7L^u1)mrJ?i_2`WM(1TgnvKgsXX#rFMIn#wiONw? zRyZ(UcVLrs&ud&PA)Jt;xTHHS3QMwXhcyL4P;@%6xF|?E6P4w}EY9$izA`Gvdf?2f zM<>V6y)kq4jj3ZVPo51qCG0a2{yGVxpoAescNNUh88CSU4~2~L@ls^j$J2caqS6=c zdlh^`?z0Y*XWa_a35Li zGfF%jlY1AJ^c`d_XRxcy5Lc6$p)+if+yDc)ivG11`RdpC+s9S!9?iQak0^SQdf;!& zHl0%a6PkY_%S_~I8*h%kH-3HM#zY=D=}{uqmr1xmSXTJ`1-Vb)F!EPKzw;Gt# z0+ZjZNDZ`Wf%dNhz4rsXpCz@viE_hHV!;u$W<;wQ$@)g}cFH~ajmy9O?A3Qx-$^}} z7VfunWm~$)2Od`a!mANIh9ANBAl)euciB1){u654L{YRt zLkz>!H(l0P0WoZ5#R-^bLU5icqxp)fqRd3VRv5(;C6lfsWsP_%ZgZ4{a$7~YM=`~w zB$KpPlzWqIDUkF?RY|Ypwc3%Slg?rIwrg)DX~1zi_m?7$pXMZhfg)BIIL?WH8U$V0 zIKcWWfEfYa-v7gRj4wun?;FM;63!Gunh_EqkjT-6XpECgKx_i|(!2m9mct>QGJ^30 zl_ej{ddMMR27!Tc2EycAJkIN0c!vmOBz%y--$PO%-&TkxHDC=@fOcgwrAC|wsCvnO zL9Bw1Q{bM77{D;@pz!V``i+E)jUq@{5EAwv-DjYlF$mH=$ZChb{62_Pl;4Z$8dj$) zAm*YE<`IZBCg-nNfBu@7I-~h_ui10Kh7w)@z}2><~uDKu7vFTuXL9H9q z>IT<559*&uKcDRxQR_#w`q4FS&RdJI>tL#xJThML$nEXSJby=An@V-4zP4|3-o~t> zQECCp9#bCa8F+?GzKyV-;`F1RfSXnDY-$^xU1wp(l!6hya9Mp_a`04=ir7}3$i4Xr z#{49X&fBewZyEYqu_x_HiDx2|L=|;_d-INEh=x$pCLJ&{T^^A1XXtiC?-a7saVk~C z?*i?1-evu|t=cCkGtMUMVcQkkMNA@4+zVw)8H0^n(H6E*2(Zr0eOAOegR`1wUOd8b z28|#92l$2;WCb-cc3QxTSl|s|V29kCAQ!5RN7%S1uwr6i7PDZv$fx*jw$Lmh#TVFq zOj!$+u&)Bml!0{t-P;aQ*$yBrMK%cv2uQP0WuIK&mc-hpVZN38)HDQTrVA!7v2V>s z!}A52ZAKMiWUR~+zJiZpe#myr0#gA?&q)aZsJKZz3<>6vFnGpzbh9sVF$g$421kbC zEy%oqT=?We!mIEh4`SK*3Xtj{cG(Bo#|)Z2#6H9xVqlz0!fQ&m1F6?(!u@qR44Y&{ zw*w66^rA#cB)o41K*T<4kU9cv9>B!{`IoUyybcH|Nv*gd7AN5^{}}|(tyaX;tcPx& zRb6{D7ol5ich7IQe|J!AKdQAK1Tbw%`SZv|AAG=6t%>XMnz$ya%r1@D zm1TB8Ue-Tw`_=SN=90Q=zqV^XB;i(>L5&&AGJ^nN>u0VFr-riLgK*y`GSjrN|Cc+B zcOyURI7kr8f?Rk%PLS(w@bdJ?Wx^wy)MMr`yv@v7wnG&_BG^ZiPaU3s0*OyS zA50vr0v~8jvf@Kw-?lY2Jf*|$|3Mwzu&{9E^PQW~!uR|Tj$TmxpE)@u0n1g`tJz6% zhy~l=F32O^AEJpnOWio5AB+81^nnQ34d*4BG3>lDUa%#9-Hq#E@h1G`pMn5?-HH9R z=WbASwQH_+?4Dh@re@eHGRpJs$B}jSVU5C;-NY6Vn#AUPPxeKKyPG`Dd>#k>SHTS^h}B8zBF0A?SsmOL zR?k@w-~dQWgkwTnF&!#I#kta120o0~BIM##U_=U*C@GL@h;uSThv=iKrn~h%<%zZpXv<<@~mM%QE1Mq(|1Y!bFur8D6Lj4M5Wv14Z$fIQl z)~}(P6osxekOOPm!#$v2Iz;LAxWtESx~;cgcgjl(vv8z=kBbW7^|5<+P7sAFi_*yQ zZu2grcQ_V@A(lsaOBJp{Gr0%E-{H54vh|~IXI^|?P93`b_KmmGuc!_Cj7+m6&bQn& zh1ceEd2OCC4$2Fd6QW@olUN}5QidCynLhFQ%;fm=2_OvEeaYF6u@my)bxM{C9Jg0c z1{PPvA_0Fn0sbb>(vI1t2>V8X{eWXqPoUZd)od%r&{}LTFH#vSyF5J^DbP* zpvIQG8^0l|I7aQHnzpk5il52?sP{w*a2bP|cBX_|-n<*9e8}m^SK&8g6(h*EW0lD} z>#0LIw=ZwU&*oie`c_wl&V)Bkd~#9i7*d;u?w$L*?JN3=F73#)+B{9_eM7tRG}O-{ zt2j_Y4HI4X*}5lvBGY!~l-k<6QM(cTO_SO>l($0(X${XjvWm_fRCj71kKnN}aIfxH Qhso<#$I0V~f~ zH%?m^3$(F3vyBznp!z1J4C-5-S>LACw?q9h*1}rZGS#n63|cncnkH6 z51?A~JpYs+v{VaAAZ0AnDO_o5{}lTFOyz+jZhr5Z@>J+=3#}g}$?eC-kSSAvqrz2E zuce{A1MYsE)TPXY`h_co)Ye@#CXFyd^Mu$9;O945_0+&#~l`pf?>ZnfheVLjX| z+i^|mpKndB>TVj7NSPSuYWP&kC^*isS6Pl{sl+Tx&2#Yt8;Mb>ilL^t_#73FvHmEV zNU$80h(kq=Rqvt=omv?jf7c0qD7Zs!Y65<+TBOZunBFOJaxI+d8k#?0X-Lq z#rnoC1mu8-}EXOT8jy)GEFqqQ8 zp`l&EqWqflb2?BvS;A+}nj?HPevRd_CVGCJjWJm}6T$OCa}hAWw&2+4_(Ujpa`e<> zwk%A~M-p^2vczVMax24I@klIdAgfqYKc~7>P0L^{n9vyU6=295S3) zR|L%jETWLQ%CeUHWwvt5+pJ0cAhI@^4&hE+^or;I$mDDwQqU-D;tm#;&n`h^`&RjSY8sQ={aV|o`LhItWaF;S8 zp+r0$xpYQ}x0|RcE@Y zbJMd|@*H?nv(bGb-F-srJ|)5F8Q$=`oc6pddd4Kr*b}|3y5=haid>-n3zVA#Hyl75-%QqgkQrxY!ny+Jd6}gk(P<;PdmkmIr~m zfpxc7=a=gID=%ki_T&hgqvi=JR&|^8%@6Gob!3AYN>f836_lu;SbsvQKe17NCS89f zQ&n?k?AF-(6SpUT#|9Vh*x*{VOtw0?qxXCfTUieN5( zGV{@l5Ev5Mj!SLFMSD=P2L(GwD44M|3WmlnVe$!$V~U(8?3n%)T-+fdieq}|jb$A} zF6*&7>H~LFpFmBO<1N+Q1BViQm0J$2Gf*)8OAZ2@-JAP?UXc@))M?Euqgw*!h|w>B ztFTmyS|p0xtE8SW{1kfqOzk-d&hr59aPv)UMvb%uUSYLyik*(cSZXdFWy1?mnxn3< zk(t>91(s`D%POm~=QV!lYn}pwy`{T|Y7l6OUjWlRP0i6f&w^2!rdT>WOVJ6i=%DlY zL?p40h{R?nI!}gGYzcdPo@1dALxC@_5Q`*uDojVCaL>@u1T>}EPtiW5G8bO}O~xS4 zU>2@8mPyR=^60)RBc-(=szzSi&cJ~=Mk*G2Cph1bUCE zMzD`MeDpxz$Pv{T!h_!4z>&ktSQgl9e1TT%7cNwwRpT>q$otN^3YH>tVR&@p%+zGo zxBxbV%j%;n*a$2}sErLlWnSs?A$b|`Ts5Lnk$H1OEIL5 zHeb0N-X2B`-h%>^atq7wzLK6Dt!LOYQxWDa!|A_r}(PVuSyYORVn|Exjuhn{6HI z6OUgN+lHjJAyD;e?V_bsu(W1uU=JF&M#QjjE}n=?XKisXeo?TquqMMu4EY0%4{MH} zfovg)0NMhyhXKLbOMkm;WOU0W0Ji3Ysw>bJoizT0xL@+MlZ2*@AxUmm2flIRRf7AW zl}H-5hlEhWq}2<{CICjJD_Bo(Wb8oFyuG%9cEIIlDyJlo)N!b-8S5)k%#}20^i>H| zlxy7_V=UW)H^%xdku;wwf!W$Vmn})lPl%)*jeXgokw_VOvsG#_m!_6YZmXpVBQYlp zNh8*df0r}>$gtCMFm3qLm}F@VbPsZP@nrLqAL}5NyNX<9R8Oh@ zCx4|yw7~6<>48?GXetqzQ*eP&C&Cv)2Jn$%d?XRMs#Mxe`N~XSuN^VV05_n~io9!t zbmSq(+?k(gGvM&7Wn%SMUtXfCOTS+z^O~#N^ZW@lx$m+yd)EaEHGVnWx%M(@`GX7c@#H-I{E_~T$(z_tenXVnb^GH#4v`(9)0%n?w4b$3>-7~#c=b<^Sy z%3oMJ0tK9+rA@N536?hexOvSelsCgES}4gv2^K12bFQ8dz--$ZaypBnJ>#JsIPN-9 z?D~1p(<^y;SB)9BXYItT{?-0Wl{e$`ZaAs5lS);K&V7<|pWxh=af45xuWX0cYp=fl z((3T)@aI+D`;)2a4=&!jC|2!T9nMt(>70x3Ha@W4waWAR$OpgTo&gC;&w%6^SQ*XK z?-e`;QnOINxiY$Gugo|v{`I2Zya*?N7e_e&Dn~hhKBf4rg>bk8*U017p@37gpONfm z1pApxOUH_1byl=A3zp`L%_A5*vRh&YUpItjCPdR=Y}!G<=un|d7H%+P-6mI76Uj8dH^37q$D9GC%X>32VJM2<<>(BbQD*kJ!t}? zt7nYGeaKD%;-LVCQmV^QOnIw2a`P;5uWm9)9T;#E{VN!DsR0MEuw0Brb`*>nJty00 zdW*GIKtRD1qi?Dxf>&yeFDUUeIz|;s9*m(J>Vo>J6QT4gN92&2#nCuwjs zMP_1gFohTrMrmFNjA|U9Vh2K|@NF4?#foZ1r)8Jc=J)PY7xR~ISiLN)Fb)O#2H_@6-mih92i8FN?Ncm-=}i_x1Spa4z1&nH;w z){cmly$}{lb*Ijy4&Dm@O0z;#2#S@LGM+lng?7&x1YR0fUddE=*WOH5_=F1I`o+fs z8~ta~{b$Agb5j4gwEvvoKL<*A?e$bNHJ6%`8V`!jgR9nzvn6#za&~Swd(+O|N28+i znB+VrIFJ1PVhiCY9y+g4h{#6i;qOZily`{Y2+sD3XZzHa7Evtg zAc1B8s#}g-&@W0HG32|%_Wf<2MA0hRr(-EC$J#p*lCLw6>%OOh=oT^!|Av#TMbW6j zz)Is%kP;*ry`tnBK*8&aR5=RiZaI?1(%x?#1sVp0#%aK^XqL7djFBn9vOCl>GnOrw z7vf{jBK0^I2iq~tVIMysTP}>4XBXz^SOA9x`S=1C zW&vl;B}69ApYFse(3+5#|9AKRn9DCh?j{dK(w%XOIqJxr#y z=neEht0_TS@1oYUPC<l?Vt_mj7c%*LwFly;^yv_Ezot-rL@_;DeF7BOi?3 z8x@>wIeorfp+(^e8$fU%^e|ajJXk6rvI}HiNg0ZnI`azTmeyE8EJJy96D5 z0)IS;yxt&wyO*eMdf>n7|G8ai9~SFQOLeCOL+$5=+6_Zf+R(IVvHjTbBgdxQS&-Xi z?~v4ccBA*Tbnk0o@9R?U>sZc^P^v@_t+E`h6tMP0Z?KkqML_XHZ!(&`BA|dYAQa_9 z6Sb}rsypDkF^=7&GGkAUfCAhkSKWqdZ`!q2c;SfXIx4x2u9!2n>b2^$oKdvbtsRo= zdp7LtX?y#6t!VF-?A-!BGmhG|Lum&kIH*ib~P|V#7S8Wg;X@(K3WUNTl$iLEAZWlfMSwJ6v1At_2;UR#P%_6tyWluczu$%O7Lz? zT$|8(D~wL@2P%s5Z%oM+s~?cavKDMK@(>d!y`mhu@#kRco3@V;@%h&c&g9c<2UhCQNEQiD>DD#Zg|^)uOdIYASvw%&`Lsk zw@bNLc9hF5`u~P<0Z$#a;x(g8bjn^BQ{hTGI`GsBk`x(h(g3-mzfT&A=XD3!nF=%u zS6b>?Y2getY!)o04ZBjmvo(|NZx<`qs^wIB2YI)mZ)&GZKJ!!>DdZw-IZ9Cvul=@d zT0VY2-uw?OC))SlP9LywZx7@rO|UNP$`|%J;{YFO7xTSEJC$oW#2V01(14IRbfDk` zQ?YCMCdLV^Dp0G~FF>DGXsF509|G@{3h!qb*EgW;2HGx6`wunxo@88wezap^s;*#^ zU!s?*wYBibi)g6*OaQMn3a_I~?XG#%t`yHUz10>tX1rTwrtax2fWKNKha>@4_4o=_ z#54-4Hqh6>$YM|bBu>gjrYQ&?a*(Q+hy8?M2*6!Nc|#5XspujCiLl!nwqwZX;w;){ zSU`qyR4l(85GVNZSxelIS>W<$E?+6X`>q5cs9OdL8wer6R;UEfNMS{ZB<|mTp@AN2*1FB;! z;sIeQK+6#>pIJk|2ZsXXBw6k7+zs)fi3sDLQD!0@~TLPI|Fn+)nOdOjly+t3Qa~qewJ6Y?@-@Ient4a(MM<)ctPU!1YUS)_#wvVZdKESRAssquh;FY1-N%>pzloFnh>jAkg8r-v1ZKHmEXJZ{TtuE-?|q0AaF0R9!%Hyg*yLc zxkvE!h~>Rfd9PsU%{XgT-dxmt76ULk9PhXYb;|3>SHbn6MR^_0|l3Lc9U z!xEJCVaYxWz1Bfg0*aM#=(QSB*mh(86I4L#p|NA5@j$xqfY{h0HTFQdCS$7-JY8v9 zmtgC9bn@|$M-z`GU_bH0)`wB4?bt@!$#mOEv29pt8x~tmOD(5G`)R@6gliN&)Xml| z%+&F8>v6I5gw%Qhf71{uX+uGiHWV~zMOUqI_yi~bYiM)s-n4!1!^w5`hp&J9x@Zqb z_JCjyY`SY!uits+);p=NSi4`U-7mU3BzMP#`(WCAP&hOqy3a`NGq5ro{hKwEu;-{) z(T??y3sM?n{ z?-R`XGB(#*+k=j~9Uu7b`32iPI5%AlLSw(^8jxHAf_)&f&%ZwN@%IHw{n`jhdD%QE z(2Mi2G^?CP;j~!Ash);c#vATpBf7ersv4xCaDLu3N)vbmfCJ=*RDQz`|jHxEZu#aYJm8+eS zty$PNDA>@6-}i~Ge#zA@*!we`J&y*Y&f{p{uZR|3X5Rsj-I`4_fDD>aR>9^I3_dxg zwgBAZ~d$H|Z1650kIoqFv4M`(*_k|LY7Qa}%n^yl>Oi(A5Pj?fmT$quq^YfJR-@X?$eesN3q%@Nw-_er9z zBh%QUoXv-GHq%}m#4Pvf1cJ)BPLkyQjMKe(ey#WRMWMPgb#Pt(@lm0|pVOmPK=#I) QE8(0G-Q str: "complement": complement, "reverse_complement": reverse_complement, } - - -def run_dna_rna_tools( - sequences: str or list[str] or tuple[str], procedure: str -) -> str or dict: - """ - Process nucleic acid sequence by one of the developed tools.\n - Run one procedure at a time: - - Transcribe nucleic acid sequence. - - Convert nucleic acid sequence into its reverse counterpart. - - Convert nucleic acid sequence into its complement counterpart. - - Convert nucleic acid sequence into its reverse-complement counterpart. - - All functions are letter case sensitive.\n - If only one sequence provided - *sequences* can be string.\n - If more - please provide *sequences* as list or tuple.\n - If more information needed please see README or desired docstring. - - Arguments: - - sequences (str or list[str] or tuple[str]): sequences to process. - - procedure (str]: desired procedure: - - "transcribe" - - "reverse" - - "complement" - - "reverse_complement" - - Return: - - processed_sequence (dict or str): Sequences, processed with desired tool. - """ - sequences = check_user_input(sequences, procedure) - processed_sequences = {} - for sequence in sequences: - processed_sequences[sequence] = DNA_RNA_PROCEDURES_FUNCTIONS[procedure]( - sequence - ) - if len(processed_sequences) == 1: - return list(processed_sequences.values())[0] - return processed_sequences diff --git a/src/fastq_tools.py b/src/fastq_tools.py index 5adf85e..2105f71 100644 --- a/src/fastq_tools.py +++ b/src/fastq_tools.py @@ -1,50 +1,12 @@ import dictionaries -def run_fastq_filter( - seqs: dict[str, tuple[str] | list[str]] = { - "@SRX079804:1:SRR292678:1:1101:21885:21885": ( - "ACAGCAACATAAACATGATGGGATGGCGTAAGCCCCCGAGATATCAGTTTACCCAGGATAAGAGATTAAATTATGAGCAACATTATTAA", - "FGGGFGGGFGGGFGDFGCEBB@CCDFDDFFFFBFFGFGEFDFFFF;D@DD>C@DDGGGDFGDGG?GFGFEGFGGEF@FDGGGFGFBGGD", - ), - "@SRX079804:1:SRR292678:1:1101:30161:30161": ( - "GAACGACAGCAGCTCCTGCATAACCGCGTCCTTCTTCTTTAGCGTTGTGCAAAGCATGTTTTGTATTACGGGCATCTCGAGCGAATC", - "DFFFEGDGGGGFGGEDCCDCEFFFFCCCCCB>CEBFGFBGGG?DE=:6@=>AD?D8DCEE:>EEABE5D@5:DDCA;EEE-DCD", - ), - }, - gc_bounds: (int | float | tuple[int | float] | list[int | float]) = (0, 100), - length_bounds: (tuple[int]) = (0, 2**32), - quality_threshold: (int) = 0, -) -> dict: - """ - Filter out fastq reads by several parameters: - - GC content - - length - - mean nucleotide quality - Please provide GC content bounds in percentages.\n - Please provide threshold quality in Phred-33 scale.\n - For GC content and length bounds: if one number is provided, bounds from 0 to number are considered.\n - - Arguments: - - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered - - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds - - length_bounds (int | tuple[int] | list[int]): read length thresholds - - quality_thresholds: read Phred-33 scaled quality thresholds - - Return: - - seqs_filtered (dict): similar dictionary as input, bad reads are filtered out - """ - seqs, gc_bounds, length_bounds, quality_threshold = check_user_input( - seqs, gc_bounds, length_bounds, quality_threshold - ) - return fastq_filter(seqs, gc_bounds, length_bounds, quality_threshold) - - def check_user_input( seqs: dict[str, tuple[str] | list[str]], gc_bounds: (int | float | tuple[int | float] | list[int | float]), length_bounds: (int | tuple[int] | list[int]), quality_threshold: int, + verbose, ): """ Check if user input can be correctly processed\n @@ -67,7 +29,9 @@ def check_user_input( raise ValueError("Invalid read sequence given") if not set(seqs[seq_name][1]).issubset(dictionaries.QUALITY_SYMBOLS): raise ValueError("Invalid quality sequence given") - return seqs, gc_bounds, length_bounds, quality_threshold + if verbose != True and verbose != False: + raise ValueError("Invalid *verbose* argument given") + return seqs, gc_bounds, length_bounds, quality_threshold, verbose def fastq_filter( @@ -75,6 +39,7 @@ def fastq_filter( gc_bounds: (int | float | tuple[int | float] | list[int | float]), length_bounds: (int | tuple[int] | list[int]), quality_threshold: (int | float), + verbose, ) -> dict: """ Parse checked input and filter out bad reads.\n @@ -92,9 +57,9 @@ def fastq_filter( for seq_name in seqs.keys(): seq = seqs[seq_name][0] seq_qual = seqs[seq_name][1] - gc_result = is_gc_good(seq, gc_bounds) - len_result = is_len_good(seq, length_bounds) - qual_result = is_qual_good(seq_qual, quality_threshold) + gc_result = is_gc_good(seq, gc_bounds, verbose) + len_result = is_len_good(seq, length_bounds, verbose) + qual_result = is_qual_good(seq_qual, quality_threshold, verbose) if gc_result and len_result and qual_result: seqs_filtered[seq_name] = seqs[seq_name] return seqs_filtered @@ -104,8 +69,7 @@ def fastq_filter( def is_gc_good( - seq: str, - gc_bounds: (int | float | tuple[int | float] | list[int | float]), + seq: str, gc_bounds: (int | float | tuple[int | float] | list[int | float]), verbose ) -> bool: """ Check GC content of a given read.\n @@ -123,12 +87,15 @@ def is_gc_good( gc_bounds = (0, gc_bounds) gc_content = seq.count("G") + seq.count("C") gc_content = gc_content / len(seq) * 100 - print(f"Read: {seq}") - print(f"GC Content: {round(gc_content, 4)}") + if verbose: + print(f"Read: {seq}") + print(f"GC Content: {round(gc_content, 4)}") return gc_bounds[0] <= gc_content <= gc_bounds[1] -def is_len_good(seq: str, length_bounds: (int | tuple[int] | list[int])) -> bool: +def is_len_good( + seq: str, length_bounds: (int | tuple[int] | list[int]), verbose +) -> bool: """ Check length of a given read.\n If one number is provided, bounds from 1 to number are considered.\n @@ -142,11 +109,12 @@ def is_len_good(seq: str, length_bounds: (int | tuple[int] | list[int])) -> bool """ if isinstance(length_bounds, int): length_bounds = (1, length_bounds) - print(f"Read Length: {round(len(seq), 4)}") + if verbose: + print(f"Read Length: {round(len(seq), 4)}") return length_bounds[0] <= len(seq) <= length_bounds[1] -def is_qual_good(seq_qual: str, quality_threshold: int | float) -> bool: +def is_qual_good(seq_qual: str, quality_threshold: int | float, verbose) -> bool: """ Check mean nucleotide sequencing quality for a given read.\n Reads with mean quality less then threshold are filtered out.\n @@ -160,5 +128,6 @@ def is_qual_good(seq_qual: str, quality_threshold: int | float) -> bool: - condition (bool): True - quality is above threshold, False - read is to be filtered """ mean_quality = sum(ord(symbol) - 33 for symbol in seq_qual) / len(seq_qual) - print(f"Mean Nucleotide Quality: {round(mean_quality, 4)}{NEW_LINE}") + if verbose: + print(f"Mean Nucleotide Quality: {round(mean_quality, 4)}{NEW_LINE}") return mean_quality > quality_threshold diff --git a/src/protein_tools.py b/src/protein_tools.py index 995dc2d..df92cef 100644 --- a/src/protein_tools.py +++ b/src/protein_tools.py @@ -287,53 +287,3 @@ def check_and_parse_user_input( procedure_arguments["nucl_acids"] = kwargs["nucl_acids"] procedure_arguments["sequences"] = sequences return procedure_arguments, procedure - - -def run_protein_tools(sequences: (str, tuple[str] or list[str]), **kwargs: str) -> dict: - """ - Process protein sequence by one of the developed tools.\n - Run one procedure at a time: - - Search for conserved amino acids residues in protein sequence - - Search for alternative frames in a protein sequences - - Convert protein sequences to RNA or DNA sequences - - Reverse the protein sequences from one-letter to three-letter format and vice-versa - - Define molecular weight of the protein sequences - - All functions except *search_for_alt_frames* are letter case sensitive\n - If only one sequence provided - *sequences* can be string.\n - If more - please provide *sequences* as list or tuple.\n - Provide protein sequence in one letter code.\n - You can obtain one letter code from three letter code with *three_one_letter_code*\n - If more information needed please see README or desired docstring - - Arguments: - - sequences (str, list[str] or tuple[str]): sequences to process - - procedure (str]: desired procedure: - - "search_for_motifs" - - "search_for_alt_frames" - - "convert_to_nucl_acids" - - "three_one_letter_code" - - "define_molecular_weight" - - For "search_for_motif" procedure provide: - - motif (str]: desired motif to check presense in every given sequence\n - Example: motif = "GA" - - overlapping (bool): count (True) or skip (False) overlapping matches. (Optional)\n - Example: overlapping = False - - For "search_for_alt_frames" procedure provide: - - alt_start_aa (str]: the name of an amino acid that is encoded by alternative start codon (Optional)\n - Example: alt_start_aa = "I" - - For "convert_to_nucl_acids" procedure provide: - - nucl_acids (str]: the nucleic acid to convert to\n - Example: nucl_acids = "RNA"\n - nucl_acids = "DNA"\n - nucl_acids = "both" - - Return: - - dict: Dictionary with processed sequences. Depends on desired tool\n - Please see Readme or desired docstring - """ - procedure_arguments, procedure = check_and_parse_user_input(sequences, **kwargs) - return PROTEINS_PROCEDURES_TO_FUNCTIONS[procedure](**procedure_arguments) From 9edb9431b4ff3c39343dd79d41f34c31eb766849 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Sun, 8 Oct 2023 18:55:24 +0300 Subject: [PATCH 08/27] Rename fastq_tools.py to fastq_filter.py, update bioinf_tools.py --- bioinf_tools.py | 6 +++--- src/{fastq_tools.py => fastq_filter.py} | 0 2 files changed, 3 insertions(+), 3 deletions(-) rename src/{fastq_tools.py => fastq_filter.py} (100%) diff --git a/bioinf_tools.py b/bioinf_tools.py index 1986edf..a240c04 100644 --- a/bioinf_tools.py +++ b/bioinf_tools.py @@ -1,5 +1,5 @@ from src import dna_rna_tools -from src import fastq_tools +from src import fastq_filter from src import protein_tools @@ -136,9 +136,9 @@ def run_fastq_filter( length_bounds, quality_threshold, verbose, - ) = fastq_tools.check_user_input( + ) = fastq_filter.check_user_input( seqs, gc_bounds, length_bounds, quality_threshold, verbose ) - return fastq_tools.fastq_filter( + return fastq_filter.fastq_filter( seqs, gc_bounds, length_bounds, quality_threshold, verbose ) diff --git a/src/fastq_tools.py b/src/fastq_filter.py similarity index 100% rename from src/fastq_tools.py rename to src/fastq_filter.py From 3e68daac7725301575a37fded90400ba7fc67d9d Mon Sep 17 00:00:00 2001 From: Vladimir Grigoriants Date: Sun, 8 Oct 2023 21:16:22 +0400 Subject: [PATCH 09/27] Update README.md --- README.md | 74 ++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 73 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 39f073e..5f0260c 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,74 @@ # IB_Bioinf_tools -This is the repo for the fifth homework of the BI Python 2023 course. It has few bioinformatis tools +## Tools to work with nucleic acids, protein sequences, fastq reads. + +Every bioinfomatician at some point of his/her career deals with nucleic acid sequences, protein sequences. Nowadays it's also hard to imagine bioinformatics without New Generation Sequecning methods and data, that needs to be analyzed with bioinformatician hands. For those reasons as well as for practical experience and better theorethical understanding of bioinformatics I create this repo. +`bioinf_tools.py` is an open-source program that facilitates working with bioinformatics data. + +## Usage and options +The programm is based on three main functions:\ + +**1.** `run_protein_tools` - facilitates work with protein sequences, has 5 procedures: + +- Search for conserved amino acids residues in protein sequence: **search_for_motifs** +- Search for alternative frames in a protein sequences: **search_for_alt_frames** +- Convert protein sequences to RNA or DNA sequences: **convert_to_nucl_acids** +- Reverse the protein sequences from one-letter to three-letter format and vice-versa: **three_one_letter_code** +- Define molecular weight of the protein sequences: **define_molecular_weight** + +Function works with *one-letter amino acid sequences*, a name of procedure and a relevant argument. If you have three-letter amino acids sequences please convert it with `three_one_letter_code` procedure before using any other procedures on them. + +To start with the program run the following command: + +`run_protein_tools(sequences, procedure: str ="procedure", ...)` + +Where: +- sequences: positional argument, a *str | list[str] | tuple[str]* of protein sequences +- procedure: keyword argument, a type of procedure to use that is inputed in *string* type +- ... - an additional keyword arguments. Those are: + For "search_for_motif" procedure: +- motif (str): desired motif to check presense in every given sequence. Example: motif = "GA" +- overlapping (bool): count (True) or skip (False) overlapping matches. Example: overlapping = False (Optional)\ + For "search_for_alt_frames" procedure: +- alt_start_aa (str): the name of an amino acid that is encoded by alternative start codon. Example: alt_start_aa = "I" (Optional)\ + For "convert_to_nucl_acids" procedure: +- nucl_acids (str): the nucleic acid to convert to. Example: nucl_acids = "both" | nucl_acids = "RNA" + +**2.** `run_dna_rna_tools` - facilitates work with nucleic acids sequences, has 4 procedures: +- Transcribe nucleic acid sequence: **transcribe** +- Convert nucleic acid sequence into its reverse counterpart: **reverse** +- Convert nucleic acid sequence into its complement counterpart: **complement** +- Convert nucleic acid sequence into its reverse-complement counterpart: **reverse_complement** + +Function works with nucleic acids sequences (DNA/RNA) and a name of procedure. + +To start with the program run the following command: + +`run_dna_rna_tools(sequences, procedure = "procedure")` + +Where: +- sequences *(str | list[str] | tuple[str])*: positional argument, a of nucleic acids sequences +- procedure *(str)*: keyword argument, a name of procedure to use that is inputed in *string* type + +**3.** `run_fastq_filter` - facilitates filtering of sequncing reads by several parameters: +- GC content: **gc_bounds** +- read length: **length_bounds** +- mean nucleotide quality: **quality_threshold** + +Function works with default arguments. +To start with the program run the following command: + +`run_fastq_filter(seqs, gc_bounds, length_bounds, quality_threshold, verbose=True)` + +Where: +- seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered. Default value: *an example dictionary* +- gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds. Default value: *(0:100)* +- length_bounds (int | tuple[int] | list[int]): read length thresholds. Default value: *(1,2**32)* +- quality_thresholds (int | float): read Phred-33 scaled quality thresholds. Default value: *0* +- verbose (bool): add detailed statistics for each read. Defaukt value: *True* + +Please provide GC content bounds in percentages.\ + +For *gc_bounds* and *length_bounds* if one number is provided bounds from 0 to number are considered.\ + +Please provide threshold quality in Phred-33 scale.\ + From c84a045e5d3aca96552b78564905783b7ea20e9d Mon Sep 17 00:00:00 2001 From: Vladimir Grigoriants Date: Sun, 8 Oct 2023 22:40:44 +0400 Subject: [PATCH 10/27] Update README.md --- README.md | 66 +++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 62 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 5f0260c..1571b44 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ Every bioinfomatician at some point of his/her career deals with nucleic acid se `bioinf_tools.py` is an open-source program that facilitates working with bioinformatics data. ## Usage and options -The programm is based on three main functions:\ +The programm is based on three main functions: **1.** `run_protein_tools` - facilitates work with protein sequences, has 5 procedures: @@ -66,9 +66,67 @@ Where: - quality_thresholds (int | float): read Phred-33 scaled quality thresholds. Default value: *0* - verbose (bool): add detailed statistics for each read. Defaukt value: *True* -Please provide GC content bounds in percentages.\ +Please provide GC content bounds in percentages. -For *gc_bounds* and *length_bounds* if one number is provided bounds from 0 to number are considered.\ +For *gc_bounds* and *length_bounds* if one number is provided bounds from 0 to number are considered. -Please provide threshold quality in Phred-33 scale.\ +Please provide threshold quality in Phred-33 scale. +## Examples +```python +### run_protein_tools + + +## three_one_letter_code +run_protein_tools(['met-Asn-Tyr', 'Ile-Ala-Ala'], procedure='three_one_letter_code') # ['mNY', 'IAA'] +run_protein_tools(['mNY','IAA'], procedure='three_one_letter_code') # ['met-Asn-Tyr', 'Ile-Ala-Ala'] + + +## define_molecular_weight +run_protein_tools(['MNY','IAA'], procedure='define_molecular_weight') # {'MNY': 426.52, 'IAA': 273.35} + + +## check_for_motifs +run_protein_tools(['mNY','IAA'], procedure='search_for_motifs', motif='NY') +#Sequence: mNY +#Motif: NY +#Motif is present in protein sequence starting at positions: 1 + +#Sequence: IAA +#Motif: NY +#Motif is not present in protein sequence + +{'mNY': [1], 'IAA': []} + + +## search_for_alt_frames +run_protein_tools(['mNYQTMSPYYDMId'], procedure='search_for_alt_frames') # {'mNYQTMSPYYDMId': ['MSPYYDMId']} +run_protein_tools(['mNYTQTSP'], procedure='search_for_alt_frames', alt_start_aa='T') # {'mNYTQTSP': ['TQTSP']} + + +## convert_to_nucl_acids +run_protein_tools('MNY', procedure='convert_to_nucl_acids', nucl_acids = 'RNA') # {'RNA': ['AUGAACUAU']} +run_protein_tools('MNY', procedure='convert_to_nucl_acids', nucl_acids = 'DNA') # {'DNA': ['TACTTGATA']} +run_protein_tools('MNY', procedure='convert_to_nucl_acids', nucl_acids = 'both') # {'RNA': ['AUGAACUAU'], 'DNA': ['TACTTGATA']} + +### run_dna_rna_tools +run_dna_rna_tools("ATGC", procedure='transcribe') # 'AUGC' +run_dna_rna_tools("ATGC", procedure='reverse') # 'CGTA' +run_dna_rna_tools("ATGC", procedure='complement') # 'TACG' +run_dna_rna_tools("ATGC", procedure='reverse_complement') # 'GCAT' + +### run_fastq_filter +run_fastq_filter(gc_bounds=50) +#Read: ACAGCAACATAAACATGATGGGATGGCGTAAGCCCCCGAGATATCAGTTTACCCAGGATAAGAGATTAAATTATGAGCAACATTATTAA +#GC Content: 38.2022 +#Read Length: 89 +#Mean Nucleotide Quality: 36.1011 + +#Read: GAACGACAGCAGCTCCTGCATAACCGCGTCCTTCTTCTTTAGCGTTGTGCAAAGCATGTTTTGTATTACGGGCATCTCGAGCGAATC +#GC Content: 49.4253 +#Read Length: 87 +#Mean Nucleotide Quality: 33.2989 + +#{'@SRX079804:1:SRR292678:1:1101:21885:21885': ('ACAGCAACATAAACATGATGGGATGGCGTAAGCCCCCGAGATATCAGTTTACCCAGGATAAGAGATTAAATTATGAGCAACATTATTAA', 'FGGGFGGGFGGGFGDFGCEBB@CCDFDDFFFFBFFGFGEFDFFFF;D@DD>C@DDGGGDFGDGG?GFGFEGFGGEF@FDGGGFGFBGGD'), +# '@SRX079804:1:SRR292678:1:1101:30161:30161': ('GAACGACAGCAGCTCCTGCATAACCGCGTCCTTCTTCTTTAGCGTTGTGCAAAGCATGTTTTGTATTACGGGCATCTCGAGCGAATC', 'DFFFEGDGGGGFGGEDCCDCEFFFFCCCCCB>CEBFGFBGGG?DE=:6@=>AD?D8DCEE:>EEABE5D@5:DDCA;EEE-DCD')} +``` From 34b1e887e64babfa4ae9cfe285960cfa8fba3d17 Mon Sep 17 00:00:00 2001 From: Vladimir Grigoriants Date: Sun, 8 Oct 2023 23:05:07 +0400 Subject: [PATCH 11/27] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 1571b44..63a2499 100644 --- a/README.md +++ b/README.md @@ -110,7 +110,7 @@ run_protein_tools('MNY', procedure='convert_to_nucl_acids', nucl_acids = 'DNA') run_protein_tools('MNY', procedure='convert_to_nucl_acids', nucl_acids = 'both') # {'RNA': ['AUGAACUAU'], 'DNA': ['TACTTGATA']} ### run_dna_rna_tools -run_dna_rna_tools("ATGC", procedure='transcribe') # 'AUGC' +run_dna_rna_tools(("ATGC", "TGCA"), procedure='transcribe') # {'ATGC': 'AUGC', 'TGCA': 'UGCA'} run_dna_rna_tools("ATGC", procedure='reverse') # 'CGTA' run_dna_rna_tools("ATGC", procedure='complement') # 'TACG' run_dna_rna_tools("ATGC", procedure='reverse_complement') # 'GCAT' From c2cac33f5dcd70c1872f34a322a7658e99e47f05 Mon Sep 17 00:00:00 2001 From: Vladimir Grigoriants Date: Sun, 8 Oct 2023 23:05:57 +0400 Subject: [PATCH 12/27] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 63a2499..ddad3f0 100644 --- a/README.md +++ b/README.md @@ -96,7 +96,7 @@ run_protein_tools(['mNY','IAA'], procedure='search_for_motifs', motif='NY') #Motif: NY #Motif is not present in protein sequence -{'mNY': [1], 'IAA': []} +#{'mNY': [1], 'IAA': []} ## search_for_alt_frames From 4b614bed366e9d447003bddb13a37dda0a1fb338 Mon Sep 17 00:00:00 2001 From: Vladimir Grigoriants Date: Sun, 8 Oct 2023 23:09:07 +0400 Subject: [PATCH 13/27] Update README.md --- README.md | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index ddad3f0..104cd4c 100644 --- a/README.md +++ b/README.md @@ -25,11 +25,14 @@ Where: - sequences: positional argument, a *str | list[str] | tuple[str]* of protein sequences - procedure: keyword argument, a type of procedure to use that is inputed in *string* type - ... - an additional keyword arguments. Those are: + For "search_for_motif" procedure: - motif (str): desired motif to check presense in every given sequence. Example: motif = "GA" -- overlapping (bool): count (True) or skip (False) overlapping matches. Example: overlapping = False (Optional)\ +- overlapping (bool): count (True) or skip (False) overlapping matches. Example: overlapping = False (Optional) + For "search_for_alt_frames" procedure: -- alt_start_aa (str): the name of an amino acid that is encoded by alternative start codon. Example: alt_start_aa = "I" (Optional)\ +- alt_start_aa (str): the name of an amino acid that is encoded by alternative start codon. Example: alt_start_aa = "I" (Optional) + For "convert_to_nucl_acids" procedure: - nucl_acids (str): the nucleic acid to convert to. Example: nucl_acids = "both" | nucl_acids = "RNA" @@ -54,10 +57,10 @@ Where: - read length: **length_bounds** - mean nucleotide quality: **quality_threshold** -Function works with default arguments. +Function works with default arguments as usage example.\ To start with the program run the following command: -`run_fastq_filter(seqs, gc_bounds, length_bounds, quality_threshold, verbose=True)` +`run_fastq_filter(gc_bounds=50)` Where: - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered. Default value: *an example dictionary* From f243add5b6a9bfc1497d364b52f830b5cdefd057 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Sun, 8 Oct 2023 22:12:20 +0300 Subject: [PATCH 14/27] Delete useless files --- bioinf_tools.py | 8 ++++---- src/__pycache__/dictionaries.cpython-311.pyc | Bin 2027 -> 0 bytes src/__pycache__/dna_rna_tools.cpython-311.pyc | Bin 4454 -> 0 bytes src/__pycache__/fastq_tools.cpython-311.pyc | Bin 7341 -> 0 bytes src/__pycache__/protein_tools.cpython-311.pyc | Bin 15212 -> 0 bytes 5 files changed, 4 insertions(+), 4 deletions(-) delete mode 100644 src/__pycache__/dictionaries.cpython-311.pyc delete mode 100644 src/__pycache__/dna_rna_tools.cpython-311.pyc delete mode 100644 src/__pycache__/fastq_tools.cpython-311.pyc delete mode 100644 src/__pycache__/protein_tools.cpython-311.pyc diff --git a/bioinf_tools.py b/bioinf_tools.py index a240c04..8a0c045 100644 --- a/bioinf_tools.py +++ b/bioinf_tools.py @@ -68,17 +68,17 @@ def run_protein_tools(sequences: (str, tuple[str] or list[str]), **kwargs: str) - "define_molecular_weight" For "search_for_motif" procedure provide: - - motif (str]: desired motif to check presense in every given sequence\n + - motif (str): desired motif to check presense in every given sequence\n Example: motif = "GA" - overlapping (bool): count (True) or skip (False) overlapping matches. (Optional)\n Example: overlapping = False For "search_for_alt_frames" procedure provide: - - alt_start_aa (str]: the name of an amino acid that is encoded by alternative start codon (Optional)\n + - alt_start_aa (str): the name of an amino acid that is encoded by alternative start codon (Optional)\n Example: alt_start_aa = "I" For "convert_to_nucl_acids" procedure provide: - - nucl_acids (str]: the nucleic acid to convert to\n + - nucl_acids (str): the nucleic acid to convert to\n Example: nucl_acids = "RNA"\n nucl_acids = "DNA"\n nucl_acids = "both" @@ -108,7 +108,7 @@ def run_fastq_filter( }, gc_bounds: (int | float | tuple[int | float] | list[int | float]) = (0, 100), length_bounds: (tuple[int]) = (0, 2**32), - quality_threshold: (int) = 0, + quality_threshold: (int | float) = 0, verbose=True, ) -> dict: """ diff --git a/src/__pycache__/dictionaries.cpython-311.pyc b/src/__pycache__/dictionaries.cpython-311.pyc deleted file mode 100644 index 4a24e0308a0c930fc1a68b13c469785f08ebe7e8..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 2027 zcmY+EOKcle6o%(`;?!{-b`xisrp-HT;x?_@rf-_mj^j?;I?ls(t40MAa$Lt!5=Y~- zQdm_&2*id33k2$}P+(KB>=GdaU4Yab)m?$ufsjV3U@tza9sgB=_J zqr4yNz5ID>uJjzEN>Yvk!oEqYz813XSu!G0J z7{|e5d>riI32=~4fEu3!6MPD6huBcm@n}8f@cPu!A#Tl&^uEuE*Gf8$1X0a2D+4>)>I&0rqna z9OQX0&kMZxDr*2^8txPZmmUUJ^mUDwAJIw^7YbUus3@iM)W)i!w9?bN4Mnr)H?|dJ zlx}X^Q*=$~i|Z9dru26e%_)7ZhIv_~XE!i@L+LqLm&5-G=FcmAzOthzFZ*mMTEu#* zu4qZtqrZ%CtE#A=^g_i_R8)F#-BS9^)eV@n4eVNRUC~WNw-nu0dZ|>x-AcGXN-R|} ztucIxu%C>j(^5#|`exd+;%=dok`*4LnnE>|LQ+%4OtF%gw#L25rtD(M&Zg{QrsKNR zkCHVgXalN$kKg(5+sw^mHT?6V?WyZ;B&(Cs9=(~YqVU%^zm=>$F@JqK_L!5^nAmsA z$?Cg*EIsY7tR<^cvaEv1?wo6wJRr5lkk7VJda#?n?WQFKt531JGd^W#OPNlQ61r;qWr}Bkd zs+e8KmzPVqjOzTiT)vgbnsdd1>dlwZId4Tdy)eI+%gkr;MQ?s&sg%lPi?_>#+w(IE zxx!c4OWad)7809VJBhVKsZp^SiTm~YtJTD8eeK?6WoxIA$j&4dcX!t7TjkS@%4TBj z*1582)v9%?wkqX@wU*%8+D@&$wQAKWjg#BEPUy;Jo$qZ_rmc1Cj~_Cv@h9GOTPGqr z&>ac9XxjIFx2?l|aV-e5o3&sFqx-}|p0|jHJ#Q6n^SoWW!}ExE)bst~ot}4z$2{*A zH#|Qe-sAZ}@m|jli68d7PrTpr0r5f4hs1|H9}yqb!jTs(9eCMs>1%s<&}vU~$NZo% zKWOYfLE(N-xE~bm2Zj4V;eJrK9~ABfh5JF_eo(j{6z&Iw`$6G;P`Doy{=cB;YY6A! zxfVrI!%u3sNew%QG4GKYc%*PQsbRN>yGaeZRoqVscas`+M8@5uhP_|hO={R(;%-vI z?iP2G8ukHkH>qJC6nB#v_91aMsbTktyGaduK-^7g*hAuOQo|k*N7Cqv7ChAm9x3|T q9-eER_JlS%_GRJg#y4HhCZG2$$V==XPSv$y;LGl$}V&*|DzW}AY4RKOCC72eU%1|h1t6TWAEXQ=8$Mzsx= zpzrAO2)dp$x>p$zRO<%-y#n(HW<6=l{E$&S4Q5nCs<+vC1e!c)XhXt9RI*93@jC2@je>=Q&u7g_ljXt+mQ-cLY7TP^aB+i*O1PY;s6T|tyPQlowUEzYs}&!y zF{*N1E=Z-S%wn}bq+=vfX(vvlGuMOf5LsHT$`c9%uXo>&Rp{NnR%00 zSdUF%QN08E9iK{54QHXK@1V*U_nW;in1PYZ%mcHdnwCgviA7iHayGY8G6~bi&?Emy*wkD%7!Rx(fL`CZ<a@|7s?= zBB@C}xd28~lB@D+b~!mN^Hm!qZ;U2q*3>()l=-fLSCVhM^GfCiL|m4Mm<1Ci6~ZS4 zDVrhisj{3?68SZ=gQvjAP*h|rRFDX)#lol1fbs#Vgps%NuhGl<(O37Puj6{mhVYW=Qi69e6so;-Advo|xf@RS@HPp&ahpo_;v3hX;)Cz-AjXE6D8{ zWvX=3dk_q5UE5AQOzB-0jIIlMu-^#w?+1tXg2M-as1dldAGo#`xTXh2jKIi0osJIQ zzYvfLa<*Ol!i~E6w7ysL&=n(eMGIUh`_F6*ZYLfl^w4=DbYAn*x1{Xe`JMDzNl&DV zL`u6z--A%E*7toqlr}f_!(>#kuD&nn^*wl8=(z4Y4EbbCyf~p`I#bG zMR}*-pe4AObNv^sZf3!-2{ut_0pHj9m*rkIU*N!tEK42hLrY*y49+}H3OMmu^*kRHQ zV=d=F0b-rUU!by$RNGkRSJYb$at8D=9n?Aw)_AK8Hn~%ks8RQ!_<6r^wtrncP*l+Yn<2 zUIpL8oLZZmr`}&Sjs7J}V$dIlPr*PwKwm_W zw|#T*$M4^N9~>9FzZ~h?j|}WZ21@+yB|VZfBFX*8)xF5or_QHydL(T`(wn|b-)E;` zEBhlF8+i^qqo4z1JK(SsY#%wB37X}!f2caGEP>6`azO$!=Cefo zzuT3ss&brB;i?&_>u|_`Q=k6_9?}DwCFkhT+uR&XXaSP+!@z5Fka|BjdN&G$x5nYO zv^D-9rLkw9bBrAz@N2w0C@^Xx#Y5=5MqO~wRRt$_7ft8pHXM-R(AaV0X-(nCAS9^` zDRzlPN7IZY8i|fPt&-LuecWIkq)sER#yc;t$FE)8O)XCR`J2FpmscP^qS+K=T;Cpo zK#d&ESM{g|hh`KonY4+lP!M-nA%Jwi;gU`urbkqA@;%6ITvh=Iixm?(n;XjjMzfM& zx^q^btEVm1`^MMgBQJrxH{nx$2?Sz9w=G66uIuHiBR>ro-sGoF-Fvwl?ExPzhk7cA zb-fMA1ErhpvcL1^sh^~D|2e~dj-~`}Gvz?vuX~Ny8-MIJVl(@(n|raFdTh~%E$V?q zEzk$=XLb%fFjMUGidgH$8->dAePM%9ZTK!8`lD34 zWtjCZzX^J3UOe19#BU2lHBB;2O=@PZ=PWi=k4GBZCKLcxM5^B`XHOK{a?ym2vXLan z%EG5y1Y&_HHne+9=Pn!EWsU7S;4W%|X`LH2xKU`9{rwu-Z`n7_T9=ijM^z&#Djuj} z-Vb$?rT=t{{kkBHFE#} diff --git a/src/__pycache__/fastq_tools.cpython-311.pyc b/src/__pycache__/fastq_tools.cpython-311.pyc deleted file mode 100644 index fe9d9859c3ce5fd9e359d90fd86b0684df40bc83..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 7341 zcmeHLTWs6b89tOmiMm*dE#I9uOkF3I?}?qXa2&gdlctU@O`5btGXyRzQnnVEa!4wT zU3v9_KKLOIPLU!C)&T03v?{s+-mr)4VOxu>K(Qw=AW(rouman|vL|7|fFe)(&!H}) z?6|?!0(%%m{p7j*=jS=+|NifskE^OY2(Fe-o)K1p2>lf&iWez?#G`L1gl-}p@l*_* zN4!m-1e;(NY)!?QNN|`_%Pxc>rc&N6(1H`^LaV}Cxh=_eY@N^Jd}y<|wl#;R1vUH_~086&{E1S;a;8L*V$-RIrBN3i`De*p1DZo zNveZ)cj4VSx1P#&Ft$|Emb6D4unfK2A7R=(EL>yrLii#Zjj#z>kl3iWm{8a-C$h5w z8;(no5LRMK?4lG83$iTm9uH}BRGLdH2%;hnkz@~xmowYJN5jgGWJT&`mBeC9z^@nB zBpZv$B>TdJ&=4EpWaSbo2^=r8N*s3+iN+K`0^Qi09=7BYY)4d7pm`(~=kU|eY04|< zHs?VXC&$2aQ33ndlKaR7?G`sH_E6Gn)`*S8X0k?;&8jxJs9b1p8hn=T85m&YFc%Ye zwouWmfmrj3pd=)*G=$uOz;eZbWjWcHjBYlJX9(0p;!<%+M3@YsV!7{hOyFbz&(P&4 zFW`wt2qH`tw2s8tx1!2C%dvQ3qH&RvmX=$_#mfLKe5ocY;*RH{mjzM6NHO97p;Je7 zYD}k2>eRV%6EYmd9*iarF?=ky7L^u1)mrJ?i_2`WM(1TgnvKgsXX#rFMIn#wiONw? zRyZ(UcVLrs&ud&PA)Jt;xTHHS3QMwXhcyL4P;@%6xF|?E6P4w}EY9$izA`Gvdf?2f zM<>V6y)kq4jj3ZVPo51qCG0a2{yGVxpoAescNNUh88CSU4~2~L@ls^j$J2caqS6=c zdlh^`?z0Y*XWa_a35Li zGfF%jlY1AJ^c`d_XRxcy5Lc6$p)+if+yDc)ivG11`RdpC+s9S!9?iQak0^SQdf;!& zHl0%a6PkY_%S_~I8*h%kH-3HM#zY=D=}{uqmr1xmSXTJ`1-Vb)F!EPKzw;Gt# z0+ZjZNDZ`Wf%dNhz4rsXpCz@viE_hHV!;u$W<;wQ$@)g}cFH~ajmy9O?A3Qx-$^}} z7VfunWm~$)2Od`a!mANIh9ANBAl)euciB1){u654L{YRt zLkz>!H(l0P0WoZ5#R-^bLU5icqxp)fqRd3VRv5(;C6lfsWsP_%ZgZ4{a$7~YM=`~w zB$KpPlzWqIDUkF?RY|Ypwc3%Slg?rIwrg)DX~1zi_m?7$pXMZhfg)BIIL?WH8U$V0 zIKcWWfEfYa-v7gRj4wun?;FM;63!Gunh_EqkjT-6XpECgKx_i|(!2m9mct>QGJ^30 zl_ej{ddMMR27!Tc2EycAJkIN0c!vmOBz%y--$PO%-&TkxHDC=@fOcgwrAC|wsCvnO zL9Bw1Q{bM77{D;@pz!V``i+E)jUq@{5EAwv-DjYlF$mH=$ZChb{62_Pl;4Z$8dj$) zAm*YE<`IZBCg-nNfBu@7I-~h_ui10Kh7w)@z}2><~uDKu7vFTuXL9H9q z>IT<559*&uKcDRxQR_#w`q4FS&RdJI>tL#xJThML$nEXSJby=An@V-4zP4|3-o~t> zQECCp9#bCa8F+?GzKyV-;`F1RfSXnDY-$^xU1wp(l!6hya9Mp_a`04=ir7}3$i4Xr z#{49X&fBewZyEYqu_x_HiDx2|L=|;_d-INEh=x$pCLJ&{T^^A1XXtiC?-a7saVk~C z?*i?1-evu|t=cCkGtMUMVcQkkMNA@4+zVw)8H0^n(H6E*2(Zr0eOAOegR`1wUOd8b z28|#92l$2;WCb-cc3QxTSl|s|V29kCAQ!5RN7%S1uwr6i7PDZv$fx*jw$Lmh#TVFq zOj!$+u&)Bml!0{t-P;aQ*$yBrMK%cv2uQP0WuIK&mc-hpVZN38)HDQTrVA!7v2V>s z!}A52ZAKMiWUR~+zJiZpe#myr0#gA?&q)aZsJKZz3<>6vFnGpzbh9sVF$g$421kbC zEy%oqT=?We!mIEh4`SK*3Xtj{cG(Bo#|)Z2#6H9xVqlz0!fQ&m1F6?(!u@qR44Y&{ zw*w66^rA#cB)o41K*T<4kU9cv9>B!{`IoUyybcH|Nv*gd7AN5^{}}|(tyaX;tcPx& zRb6{D7ol5ich7IQe|J!AKdQAK1Tbw%`SZv|AAG=6t%>XMnz$ya%r1@D zm1TB8Ue-Tw`_=SN=90Q=zqV^XB;i(>L5&&AGJ^nN>u0VFr-riLgK*y`GSjrN|Cc+B zcOyURI7kr8f?Rk%PLS(w@bdJ?Wx^wy)MMr`yv@v7wnG&_BG^ZiPaU3s0*OyS zA50vr0v~8jvf@Kw-?lY2Jf*|$|3Mwzu&{9E^PQW~!uR|Tj$TmxpE)@u0n1g`tJz6% zhy~l=F32O^AEJpnOWio5AB+81^nnQ34d*4BG3>lDUa%#9-Hq#E@h1G`pMn5?-HH9R z=WbASwQH_+?4Dh@re@eHGRpJs$B}jSVU5C;-NY6Vn#AUPPxeKKyPG`Dd>#k>SHTS^h}B8zBF0A?SsmOL zR?k@w-~dQWgkwTnF&!#I#kta120o0~BIM##U_=U*C@GL@h;uSThv=iKrn~h%<%zZpXv<<@~mM%QE1Mq(|1Y!bFur8D6Lj4M5Wv14Z$fIQl z)~}(P6osxekOOPm!#$v2Iz;LAxWtESx~;cgcgjl(vv8z=kBbW7^|5<+P7sAFi_*yQ zZu2grcQ_V@A(lsaOBJp{Gr0%E-{H54vh|~IXI^|?P93`b_KmmGuc!_Cj7+m6&bQn& zh1ceEd2OCC4$2Fd6QW@olUN}5QidCynLhFQ%;fm=2_OvEeaYF6u@my)bxM{C9Jg0c z1{PPvA_0Fn0sbb>(vI1t2>V8X{eWXqPoUZd)od%r&{}LTFH#vSyF5J^DbP* zpvIQG8^0l|I7aQHnzpk5il52?sP{w*a2bP|cBX_|-n<*9e8}m^SK&8g6(h*EW0lD} z>#0LIw=ZwU&*oie`c_wl&V)Bkd~#9i7*d;u?w$L*?JN3=F73#)+B{9_eM7tRG}O-{ zt2j_Y4HI4X*}5lvBGY!~l-k<6QM(cTO_SO>l($0(X${XjvWm_fRCj71kKnN}aIfxH Qhso<#$I0V~f~ zH%?m^3$(F3vyBznp!z1J4C-5-S>LACw?q9h*1}rZGS#n63|cncnkH6 z51?A~JpYs+v{VaAAZ0AnDO_o5{}lTFOyz+jZhr5Z@>J+=3#}g}$?eC-kSSAvqrz2E zuce{A1MYsE)TPXY`h_co)Ye@#CXFyd^Mu$9;O945_0+&#~l`pf?>ZnfheVLjX| z+i^|mpKndB>TVj7NSPSuYWP&kC^*isS6Pl{sl+Tx&2#Yt8;Mb>ilL^t_#73FvHmEV zNU$80h(kq=Rqvt=omv?jf7c0qD7Zs!Y65<+TBOZunBFOJaxI+d8k#?0X-Lq z#rnoC1mu8-}EXOT8jy)GEFqqQ8 zp`l&EqWqflb2?BvS;A+}nj?HPevRd_CVGCJjWJm}6T$OCa}hAWw&2+4_(Ujpa`e<> zwk%A~M-p^2vczVMax24I@klIdAgfqYKc~7>P0L^{n9vyU6=295S3) zR|L%jETWLQ%CeUHWwvt5+pJ0cAhI@^4&hE+^or;I$mDDwQqU-D;tm#;&n`h^`&RjSY8sQ={aV|o`LhItWaF;S8 zp+r0$xpYQ}x0|RcE@Y zbJMd|@*H?nv(bGb-F-srJ|)5F8Q$=`oc6pddd4Kr*b}|3y5=haid>-n3zVA#Hyl75-%QqgkQrxY!ny+Jd6}gk(P<;PdmkmIr~m zfpxc7=a=gID=%ki_T&hgqvi=JR&|^8%@6Gob!3AYN>f836_lu;SbsvQKe17NCS89f zQ&n?k?AF-(6SpUT#|9Vh*x*{VOtw0?qxXCfTUieN5( zGV{@l5Ev5Mj!SLFMSD=P2L(GwD44M|3WmlnVe$!$V~U(8?3n%)T-+fdieq}|jb$A} zF6*&7>H~LFpFmBO<1N+Q1BViQm0J$2Gf*)8OAZ2@-JAP?UXc@))M?Euqgw*!h|w>B ztFTmyS|p0xtE8SW{1kfqOzk-d&hr59aPv)UMvb%uUSYLyik*(cSZXdFWy1?mnxn3< zk(t>91(s`D%POm~=QV!lYn}pwy`{T|Y7l6OUjWlRP0i6f&w^2!rdT>WOVJ6i=%DlY zL?p40h{R?nI!}gGYzcdPo@1dALxC@_5Q`*uDojVCaL>@u1T>}EPtiW5G8bO}O~xS4 zU>2@8mPyR=^60)RBc-(=szzSi&cJ~=Mk*G2Cph1bUCE zMzD`MeDpxz$Pv{T!h_!4z>&ktSQgl9e1TT%7cNwwRpT>q$otN^3YH>tVR&@p%+zGo zxBxbV%j%;n*a$2}sErLlWnSs?A$b|`Ts5Lnk$H1OEIL5 zHeb0N-X2B`-h%>^atq7wzLK6Dt!LOYQxWDa!|A_r}(PVuSyYORVn|Exjuhn{6HI z6OUgN+lHjJAyD;e?V_bsu(W1uU=JF&M#QjjE}n=?XKisXeo?TquqMMu4EY0%4{MH} zfovg)0NMhyhXKLbOMkm;WOU0W0Ji3Ysw>bJoizT0xL@+MlZ2*@AxUmm2flIRRf7AW zl}H-5hlEhWq}2<{CICjJD_Bo(Wb8oFyuG%9cEIIlDyJlo)N!b-8S5)k%#}20^i>H| zlxy7_V=UW)H^%xdku;wwf!W$Vmn})lPl%)*jeXgokw_VOvsG#_m!_6YZmXpVBQYlp zNh8*df0r}>$gtCMFm3qLm}F@VbPsZP@nrLqAL}5NyNX<9R8Oh@ zCx4|yw7~6<>48?GXetqzQ*eP&C&Cv)2Jn$%d?XRMs#Mxe`N~XSuN^VV05_n~io9!t zbmSq(+?k(gGvM&7Wn%SMUtXfCOTS+z^O~#N^ZW@lx$m+yd)EaEHGVnWx%M(@`GX7c@#H-I{E_~T$(z_tenXVnb^GH#4v`(9)0%n?w4b$3>-7~#c=b<^Sy z%3oMJ0tK9+rA@N536?hexOvSelsCgES}4gv2^K12bFQ8dz--$ZaypBnJ>#JsIPN-9 z?D~1p(<^y;SB)9BXYItT{?-0Wl{e$`ZaAs5lS);K&V7<|pWxh=af45xuWX0cYp=fl z((3T)@aI+D`;)2a4=&!jC|2!T9nMt(>70x3Ha@W4waWAR$OpgTo&gC;&w%6^SQ*XK z?-e`;QnOINxiY$Gugo|v{`I2Zya*?N7e_e&Dn~hhKBf4rg>bk8*U017p@37gpONfm z1pApxOUH_1byl=A3zp`L%_A5*vRh&YUpItjCPdR=Y}!G<=un|d7H%+P-6mI76Uj8dH^37q$D9GC%X>32VJM2<<>(BbQD*kJ!t}? zt7nYGeaKD%;-LVCQmV^QOnIw2a`P;5uWm9)9T;#E{VN!DsR0MEuw0Brb`*>nJty00 zdW*GIKtRD1qi?Dxf>&yeFDUUeIz|;s9*m(J>Vo>J6QT4gN92&2#nCuwjs zMP_1gFohTrMrmFNjA|U9Vh2K|@NF4?#foZ1r)8Jc=J)PY7xR~ISiLN)Fb)O#2H_@6-mih92i8FN?Ncm-=}i_x1Spa4z1&nH;w z){cmly$}{lb*Ijy4&Dm@O0z;#2#S@LGM+lng?7&x1YR0fUddE=*WOH5_=F1I`o+fs z8~ta~{b$Agb5j4gwEvvoKL<*A?e$bNHJ6%`8V`!jgR9nzvn6#za&~Swd(+O|N28+i znB+VrIFJ1PVhiCY9y+g4h{#6i;qOZily`{Y2+sD3XZzHa7Evtg zAc1B8s#}g-&@W0HG32|%_Wf<2MA0hRr(-EC$J#p*lCLw6>%OOh=oT^!|Av#TMbW6j zz)Is%kP;*ry`tnBK*8&aR5=RiZaI?1(%x?#1sVp0#%aK^XqL7djFBn9vOCl>GnOrw z7vf{jBK0^I2iq~tVIMysTP}>4XBXz^SOA9x`S=1C zW&vl;B}69ApYFse(3+5#|9AKRn9DCh?j{dK(w%XOIqJxr#y z=neEht0_TS@1oYUPC<l?Vt_mj7c%*LwFly;^yv_Ezot-rL@_;DeF7BOi?3 z8x@>wIeorfp+(^e8$fU%^e|ajJXk6rvI}HiNg0ZnI`azTmeyE8EJJy96D5 z0)IS;yxt&wyO*eMdf>n7|G8ai9~SFQOLeCOL+$5=+6_Zf+R(IVvHjTbBgdxQS&-Xi z?~v4ccBA*Tbnk0o@9R?U>sZc^P^v@_t+E`h6tMP0Z?KkqML_XHZ!(&`BA|dYAQa_9 z6Sb}rsypDkF^=7&GGkAUfCAhkSKWqdZ`!q2c;SfXIx4x2u9!2n>b2^$oKdvbtsRo= zdp7LtX?y#6t!VF-?A-!BGmhG|Lum&kIH*ib~P|V#7S8Wg;X@(K3WUNTl$iLEAZWlfMSwJ6v1At_2;UR#P%_6tyWluczu$%O7Lz? zT$|8(D~wL@2P%s5Z%oM+s~?cavKDMK@(>d!y`mhu@#kRco3@V;@%h&c&g9c<2UhCQNEQiD>DD#Zg|^)uOdIYASvw%&`Lsk zw@bNLc9hF5`u~P<0Z$#a;x(g8bjn^BQ{hTGI`GsBk`x(h(g3-mzfT&A=XD3!nF=%u zS6b>?Y2getY!)o04ZBjmvo(|NZx<`qs^wIB2YI)mZ)&GZKJ!!>DdZw-IZ9Cvul=@d zT0VY2-uw?OC))SlP9LywZx7@rO|UNP$`|%J;{YFO7xTSEJC$oW#2V01(14IRbfDk` zQ?YCMCdLV^Dp0G~FF>DGXsF509|G@{3h!qb*EgW;2HGx6`wunxo@88wezap^s;*#^ zU!s?*wYBibi)g6*OaQMn3a_I~?XG#%t`yHUz10>tX1rTwrtax2fWKNKha>@4_4o=_ z#54-4Hqh6>$YM|bBu>gjrYQ&?a*(Q+hy8?M2*6!Nc|#5XspujCiLl!nwqwZX;w;){ zSU`qyR4l(85GVNZSxelIS>W<$E?+6X`>q5cs9OdL8wer6R;UEfNMS{ZB<|mTp@AN2*1FB;! z;sIeQK+6#>pIJk|2ZsXXBw6k7+zs)fi3sDLQD!0@~TLPI|Fn+)nOdOjly+t3Qa~qewJ6Y?@-@Ient4a(MM<)ctPU!1YUS)_#wvVZdKESRAssquh;FY1-N%>pzloFnh>jAkg8r-v1ZKHmEXJZ{TtuE-?|q0AaF0R9!%Hyg*yLc zxkvE!h~>Rfd9PsU%{XgT-dxmt76ULk9PhXYb;|3>SHbn6MR^_0|l3Lc9U z!xEJCVaYxWz1Bfg0*aM#=(QSB*mh(86I4L#p|NA5@j$xqfY{h0HTFQdCS$7-JY8v9 zmtgC9bn@|$M-z`GU_bH0)`wB4?bt@!$#mOEv29pt8x~tmOD(5G`)R@6gliN&)Xml| z%+&F8>v6I5gw%Qhf71{uX+uGiHWV~zMOUqI_yi~bYiM)s-n4!1!^w5`hp&J9x@Zqb z_JCjyY`SY!uits+);p=NSi4`U-7mU3BzMP#`(WCAP&hOqy3a`NGq5ro{hKwEu;-{) z(T??y3sM?n{ z?-R`XGB(#*+k=j~9Uu7b`32iPI5%AlLSw(^8jxHAf_)&f&%ZwN@%IHw{n`jhdD%QE z(2Mi2G^?CP;j~!Ash);c#vATpBf7ersv4xCaDLu3N)vbmfCJ=*RDQz`|jHxEZu#aYJm8+eS zty$PNDA>@6-}i~Ge#zA@*!we`J&y*Y&f{p{uZR|3X5Rsj-I`4_fDD>aR>9^I3_dxg zwgBAZ~d$H|Z1650kIoqFv4M`(*_k|LY7Qa}%n^yl>Oi(A5Pj?fmT$quq^YfJR-@X?$eesN3q%@Nw-_er9z zBh%QUoXv-GHq%}m#4Pvf1cJ)BPLkyQjMKe(ey#WRMWMPgb#Pt(@lm0|pVOmPK=#I) QE8(0G-Q Date: Mon, 23 Oct 2023 14:33:33 +0300 Subject: [PATCH 15/27] Add fixes from PR --- bioinf_tools.py | 7 ++++--- src/dna_rna_tools.py | 6 ++++-- src/fastq_filter.py | 20 +++++++++++++++----- src/protein_tools.py | 6 ++++-- 4 files changed, 27 insertions(+), 12 deletions(-) diff --git a/bioinf_tools.py b/bioinf_tools.py index 8a0c045..dd1cd66 100644 --- a/bioinf_tools.py +++ b/bioinf_tools.py @@ -109,7 +109,7 @@ def run_fastq_filter( gc_bounds: (int | float | tuple[int | float] | list[int | float]) = (0, 100), length_bounds: (tuple[int]) = (0, 2**32), quality_threshold: (int | float) = 0, - verbose=True, + verbose=False, ) -> dict: """ Filter out fastq reads by several parameters: @@ -124,8 +124,9 @@ def run_fastq_filter( - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds - length_bounds (int | tuple[int] | list[int]): read length thresholds - - quality_thresholds: read Phred-33 scaled quality thresholds - Examples are provided as default values for each argument + - quality_thresholds (int | float): read Phred-33 scaled quality thresholds + - verbose (bool): add detailed statistics for each read + Examples are provid ed as default values for each argument Return: - seqs_filtered (dict): similar dictionary as input, bad reads are filtered out diff --git a/src/dna_rna_tools.py b/src/dna_rna_tools.py index e87de3f..6a6f3ae 100644 --- a/src/dna_rna_tools.py +++ b/src/dna_rna_tools.py @@ -1,5 +1,7 @@ -import dictionaries - +if __name__ == "__main__": + import dictionaries +else: + from src import dictionaries def check_user_input(sequences: str or list[str] or tuple[str], procedure): """ diff --git a/src/fastq_filter.py b/src/fastq_filter.py index 2105f71..77f63da 100644 --- a/src/fastq_filter.py +++ b/src/fastq_filter.py @@ -1,4 +1,7 @@ -import dictionaries +if __name__ == "__main__": + import dictionaries +else: + from src import dictionaries def check_user_input( @@ -16,6 +19,7 @@ def check_user_input( - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds - length_bounds (int | tuple[int] | list[int]): read length thresholds - quality_thresholds: read Phred-33 scaled quality thresholds + - verbose (bool): add detailed statistics for each read Return: - same arguments as input, checked for correctness @@ -49,6 +53,7 @@ def fastq_filter( - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds - length_bounds (int | tuple[int] | list[int]): read length thresholds - quality_thresholds: read Phred-33 scaled quality thresholds + - verbose (bool): add detailed statistics for each read Return: - seqs_filtered (dict): similar dictionary as input, bad reads are filtered out @@ -58,10 +63,12 @@ def fastq_filter( seq = seqs[seq_name][0] seq_qual = seqs[seq_name][1] gc_result = is_gc_good(seq, gc_bounds, verbose) - len_result = is_len_good(seq, length_bounds, verbose) - qual_result = is_qual_good(seq_qual, quality_threshold, verbose) - if gc_result and len_result and qual_result: - seqs_filtered[seq_name] = seqs[seq_name] + if gc_result: + len_result = is_len_good(seq, length_bounds, verbose) + if len_result: + qual_result = is_qual_good(seq_qual, quality_threshold, verbose) + if qual_result: + seqs_filtered[seq_name] = seqs[seq_name] return seqs_filtered @@ -79,6 +86,7 @@ def is_gc_good( Arguments: - seq (str): read to check it's length - gc_bounds (int | float | tuple[int] | list[int]): GC content thresholds, by which reads are filtered + - verbose (bool): add detailed statistics for each read Return: - condition (bool): True - GC content is within bounds, False - read is to be filtered @@ -103,6 +111,7 @@ def is_len_good( Arguments: - seq (str): read to check it's length - length_bounds (int | tuple[int] | list[int]): length thresholds, by which reads are filtered + - verbose (bool): add detailed statistics for each read Return: - condition (bool): True - length is within bounds, False - read is to be filtered @@ -123,6 +132,7 @@ def is_qual_good(seq_qual: str, quality_threshold: int | float, verbose) -> bool Arguments: - seq_qual (str): quality sequence in Phred-33 scale for a read - quality_threshold (int | float): threshold, by which reads are filtered + - verbose (bool): add detailed statistics for each read Return: - condition (bool): True - quality is above threshold, False - read is to be filtered diff --git a/src/protein_tools.py b/src/protein_tools.py index df92cef..7a248cd 100644 --- a/src/protein_tools.py +++ b/src/protein_tools.py @@ -1,5 +1,7 @@ -import dictionaries - +if __name__ == "__main__": + import dictionaries +else: + from src import dictionaries def three_one_letter_code(sequences: (tuple[str] or list[str])) -> list: """ From 24afd1ad0032a3ead57fe0bba3415f15493e7387 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Mon, 23 Oct 2023 14:41:26 +0300 Subject: [PATCH 16/27] Add fasta file parsing to dict in fastq_filter.py --- src/fastq_filter.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/fastq_filter.py b/src/fastq_filter.py index 77f63da..2e86c6f 100644 --- a/src/fastq_filter.py +++ b/src/fastq_filter.py @@ -5,7 +5,7 @@ def check_user_input( - seqs: dict[str, tuple[str] | list[str]], + sequences_path: str, gc_bounds: (int | float | tuple[int | float] | list[int | float]), length_bounds: (int | tuple[int] | list[int]), quality_threshold: int, @@ -24,8 +24,17 @@ def check_user_input( Return: - same arguments as input, checked for correctness """ - if not isinstance(seqs, dict): - raise ValueError("Please provide sequences info with a dictionary") + seqs = {} + with open(sequences_path, "r") as seqs_file: + count = 0 + for line in seqs_file: + count += 1 + if count == 1 or count % 5 == 0: + seqs[line.strip()] = [] + if count % 2 == 0: + seqs[list(seqs)[-1]].append(line) + if count % 4 == 0: + seqs[list(seqs)[-1]].append(line) for seq_name in seqs.keys(): if not isinstance(seq_name, str): raise ValueError("Invalid sequence name given") From 419d18744ca089b3c7968aa124a21cf704f75419 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Mon, 23 Oct 2023 15:05:08 +0300 Subject: [PATCH 17/27] Fix fasta file parsing to dict in fastq_filter.py --- src/fastq_filter.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/fastq_filter.py b/src/fastq_filter.py index 2e86c6f..1cf0753 100644 --- a/src/fastq_filter.py +++ b/src/fastq_filter.py @@ -29,12 +29,12 @@ def check_user_input( count = 0 for line in seqs_file: count += 1 - if count == 1 or count % 5 == 0: + if count == 1 % (count-1) % 4 == 0: seqs[line.strip()] = [] - if count % 2 == 0: - seqs[list(seqs)[-1]].append(line) + if count % 2 == 0 and count % 4 != 0: + seqs[list(seqs)[-1]].append(line.strip()) if count % 4 == 0: - seqs[list(seqs)[-1]].append(line) + seqs[list(seqs)[-1]].append(line.strip()) for seq_name in seqs.keys(): if not isinstance(seq_name, str): raise ValueError("Invalid sequence name given") From 683a12a4093bbdb50bc5bb5a754137e026c5edac Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Mon, 23 Oct 2023 15:05:37 +0300 Subject: [PATCH 18/27] Fix fasta file parsing to dict in fastq_filter.py --- src/fastq_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fastq_filter.py b/src/fastq_filter.py index 1cf0753..10ba7c8 100644 --- a/src/fastq_filter.py +++ b/src/fastq_filter.py @@ -29,7 +29,7 @@ def check_user_input( count = 0 for line in seqs_file: count += 1 - if count == 1 % (count-1) % 4 == 0: + if count == 1 or (count-1) % 4 == 0: seqs[line.strip()] = [] if count % 2 == 0 and count % 4 != 0: seqs[list(seqs)[-1]].append(line.strip()) From b95b0f28dd6e0feed14adb52f2ebcb7bd18da339 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Mon, 23 Oct 2023 18:01:13 +0300 Subject: [PATCH 19/27] Add save_filtered_seqs function in fastq_filter.py --- bioinf_tools.py | 38 ++++++++++++++++++++--------------- src/fastq_filter.py | 49 ++++++++++++++++++++++++++++++++++++--------- 2 files changed, 62 insertions(+), 25 deletions(-) diff --git a/bioinf_tools.py b/bioinf_tools.py index dd1cd66..296b7e2 100644 --- a/bioinf_tools.py +++ b/bioinf_tools.py @@ -96,20 +96,14 @@ def run_protein_tools(sequences: (str, tuple[str] or list[str]), **kwargs: str) def run_fastq_filter( - seqs: dict[str, tuple[str] | list[str]] = { - "@SRX079804:1:SRR292678:1:1101:21885:21885": ( - "ACAGCAACATAAACATGATGGGATGGCGTAAGCCCCCGAGATATCAGTTTACCCAGGATAAGAGATTAAATTATGAGCAACATTATTAA", - "FGGGFGGGFGGGFGDFGCEBB@CCDFDDFFFFBFFGFGEFDFFFF;D@DD>C@DDGGGDFGDGG?GFGFEGFGGEF@FDGGGFGFBGGD", - ), - "@SRX079804:1:SRR292678:1:1101:30161:30161": ( - "GAACGACAGCAGCTCCTGCATAACCGCGTCCTTCTTCTTTAGCGTTGTGCAAAGCATGTTTTGTATTACGGGCATCTCGAGCGAATC", - "DFFFEGDGGGGFGGEDCCDCEFFFFCCCCCB>CEBFGFBGGG?DE=:6@=>AD?D8DCEE:>EEABE5D@5:DDCA;EEE-DCD", - ), - }, + sequences_path: str, gc_bounds: (int | float | tuple[int | float] | list[int | float]) = (0, 100), length_bounds: (tuple[int]) = (0, 2**32), quality_threshold: (int | float) = 0, - verbose=False, + verbose: bool = False, + save_filtered_seqs: bool = True, + save_to_dir: str = "./fastq_filtrator_resuls", + output_filename: str = "", ) -> dict: """ Filter out fastq reads by several parameters: @@ -121,12 +115,15 @@ def run_fastq_filter( For GC content and length bounds: if one number is provided, bounds from 0 to number are considered.\n Arguments: - - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered + - sequences_path (str): absolute or relative path to desired file, containing sequences in fasta format - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds - length_bounds (int | tuple[int] | list[int]): read length thresholds - quality_thresholds (int | float): read Phred-33 scaled quality thresholds - verbose (bool): add detailed statistics for each read - Examples are provid ed as default values for each argument + - save_filtered_seqs (bool): save filtered reads to fasta file + - save_to_dir (str): absolute or realtive path to directory to save to + - output_filename (str): output name of the filtered fasta file + For examples see default values for each argument Return: - seqs_filtered (dict): similar dictionary as input, bad reads are filtered out @@ -137,9 +134,18 @@ def run_fastq_filter( length_bounds, quality_threshold, verbose, - ) = fastq_filter.check_user_input( - seqs, gc_bounds, length_bounds, quality_threshold, verbose + output_filename, + ) = fastq_filter.parse_and_check_user_input( + sequences_path, + gc_bounds, + length_bounds, + quality_threshold, + verbose, + output_filename, ) - return fastq_filter.fastq_filter( + filtered_seqs = fastq_filter.fastq_filter( seqs, gc_bounds, length_bounds, quality_threshold, verbose ) + if save_filtered_seqs: + fastq_filter.save_filtered_seqs(filtered_seqs, output_filename, save_to_dir) + return filtered_seqs diff --git a/src/fastq_filter.py b/src/fastq_filter.py index 10ba7c8..3fa7e2d 100644 --- a/src/fastq_filter.py +++ b/src/fastq_filter.py @@ -2,34 +2,39 @@ import dictionaries else: from src import dictionaries +import os -def check_user_input( +def parse_and_check_user_input( sequences_path: str, gc_bounds: (int | float | tuple[int | float] | list[int | float]), length_bounds: (int | tuple[int] | list[int]), quality_threshold: int, - verbose, + verbose: bool, + output_filename: "str", ): """ - Check if user input can be correctly processed\n + Parse input fasta file to dictionary[sequence_name: [sequence, sequence_quality]] and check if input can be correctly processed\n Arguments: - - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered + - sequences_path(str): absolute or relative path to desired file, containing sequences in fasta format - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds - length_bounds (int | tuple[int] | list[int]): read length thresholds - quality_thresholds: read Phred-33 scaled quality thresholds - verbose (bool): add detailed statistics for each read + - output_filename (str): output name of the filtered fasta file Return: - same arguments as input, checked for correctness """ + if output_filename == "": + output_filename = os.path.basename(sequences_path) seqs = {} with open(sequences_path, "r") as seqs_file: count = 0 for line in seqs_file: count += 1 - if count == 1 or (count-1) % 4 == 0: + if count == 1 or (count - 1) % 4 == 0: seqs[line.strip()] = [] if count % 2 == 0 and count % 4 != 0: seqs[list(seqs)[-1]].append(line.strip()) @@ -44,21 +49,21 @@ def check_user_input( raise ValueError("Invalid quality sequence given") if verbose != True and verbose != False: raise ValueError("Invalid *verbose* argument given") - return seqs, gc_bounds, length_bounds, quality_threshold, verbose + return seqs, gc_bounds, length_bounds, quality_threshold, verbose, output_filename def fastq_filter( - seqs: dict[str, tuple[str] | list[str]], + seqs: dict[str, list[str]], gc_bounds: (int | float | tuple[int | float] | list[int | float]), length_bounds: (int | tuple[int] | list[int]), quality_threshold: (int | float), verbose, ) -> dict: """ - Parse checked input and filter out bad reads.\n + Filter out bad reads.\n Arguments: - - seqs (dict[str, tuple[str] | list[str]]): fastq reads to be filtered + - seqs (dict[str, list[str]]): fastq reads to be filtered - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds - length_bounds (int | tuple[int] | list[int]): read length thresholds - quality_thresholds: read Phred-33 scaled quality thresholds @@ -150,3 +155,29 @@ def is_qual_good(seq_qual: str, quality_threshold: int | float, verbose) -> bool if verbose: print(f"Mean Nucleotide Quality: {round(mean_quality, 4)}{NEW_LINE}") return mean_quality > quality_threshold + + +def save_filtered_seqs( + filtered_seqs: dict[str, list[str]], + output_filename: str, + save_to_dir: str, +): + """ + Save filtered reads in fasta format. + + Arguments: + - filtered_seqs (dict[str, list[str]]): filtered fastq_reads + - output_filename (str): output name of the filtered fasta file + - save_to_dir (str): absolute or realtive path to directory to save to + """ + os.makedirs(save_to_dir, exist_ok=True) + output_path = os.path.join(save_to_dir, output_filename) + for seq_name in filtered_seqs.keys(): + filtered_seqs[seq_name].append("+" + seq_name[1:]) + filtered_seqs[seq_name] = [x + "\n" for x in filtered_seqs[seq_name]] + filtered_seqs[list(filtered_seqs)[-1]][1] = filtered_seqs[list(filtered_seqs)[-1]][ + 1 + ].strip() # get rid of last \n, to avoid empty line in filtered fasta file + with open(output_path, "w") as output_file: + for key, value in filtered_seqs.items(): + output_file.write("%s\n%s%s%s" % (key, value[0], value[2], value[1])) From ec4623db191a929e1d6547ac554c097c0fcdbb97 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Mon, 23 Oct 2023 18:08:27 +0300 Subject: [PATCH 20/27] Update docstrings --- bioinf_tools.py | 7 ++++--- src/fastq_filter.py | 12 ++++++------ 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/bioinf_tools.py b/bioinf_tools.py index 296b7e2..670cd05 100644 --- a/bioinf_tools.py +++ b/bioinf_tools.py @@ -110,19 +110,20 @@ def run_fastq_filter( - GC content - length - mean nucleotide quality + Save filtered reads in fastq format. Please provide GC content bounds in percentages.\n Please provide threshold quality in Phred-33 scale.\n For GC content and length bounds: if one number is provided, bounds from 0 to number are considered.\n Arguments: - - sequences_path (str): absolute or relative path to desired file, containing sequences in fasta format + - sequences_path (str): absolute or relative path to desired file, containing sequences in fastq format - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds - length_bounds (int | tuple[int] | list[int]): read length thresholds - quality_thresholds (int | float): read Phred-33 scaled quality thresholds - verbose (bool): add detailed statistics for each read - - save_filtered_seqs (bool): save filtered reads to fasta file + - save_filtered_seqs (bool): save filtered reads to fastq file - save_to_dir (str): absolute or realtive path to directory to save to - - output_filename (str): output name of the filtered fasta file + - output_filename (str): output name of the filtered fastq file For examples see default values for each argument Return: diff --git a/src/fastq_filter.py b/src/fastq_filter.py index 3fa7e2d..828eaed 100644 --- a/src/fastq_filter.py +++ b/src/fastq_filter.py @@ -14,15 +14,15 @@ def parse_and_check_user_input( output_filename: "str", ): """ - Parse input fasta file to dictionary[sequence_name: [sequence, sequence_quality]] and check if input can be correctly processed\n + Parse input fastq file to dictionary[sequence_name: [sequence, sequence_quality]] and check if input can be correctly processed\n Arguments: - - sequences_path(str): absolute or relative path to desired file, containing sequences in fasta format + - sequences_path(str): absolute or relative path to desired file, containing sequences in fastq format - gc_bounds (int | float | tuple[int | float] | list[int | float]): GC content thresholds - length_bounds (int | tuple[int] | list[int]): read length thresholds - quality_thresholds: read Phred-33 scaled quality thresholds - verbose (bool): add detailed statistics for each read - - output_filename (str): output name of the filtered fasta file + - output_filename (str): output name of the filtered fastq file Return: - same arguments as input, checked for correctness @@ -163,11 +163,11 @@ def save_filtered_seqs( save_to_dir: str, ): """ - Save filtered reads in fasta format. + Save filtered reads in fastq format. Arguments: - filtered_seqs (dict[str, list[str]]): filtered fastq_reads - - output_filename (str): output name of the filtered fasta file + - output_filename (str): output name of the filtered fastq file - save_to_dir (str): absolute or realtive path to directory to save to """ os.makedirs(save_to_dir, exist_ok=True) @@ -177,7 +177,7 @@ def save_filtered_seqs( filtered_seqs[seq_name] = [x + "\n" for x in filtered_seqs[seq_name]] filtered_seqs[list(filtered_seqs)[-1]][1] = filtered_seqs[list(filtered_seqs)[-1]][ 1 - ].strip() # get rid of last \n, to avoid empty line in filtered fasta file + ].strip() # get rid of last \n, to avoid empty line in filtered fastq file with open(output_path, "w") as output_file: for key, value in filtered_seqs.items(): output_file.write("%s\n%s%s%s" % (key, value[0], value[2], value[1])) From ccca0f3ade2285bc998caea0eee22b749bd7bc83 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Mon, 23 Oct 2023 18:12:35 +0300 Subject: [PATCH 21/27] Update docstrings --- src/fastq_filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/fastq_filter.py b/src/fastq_filter.py index 828eaed..7e0e7db 100644 --- a/src/fastq_filter.py +++ b/src/fastq_filter.py @@ -11,7 +11,7 @@ def parse_and_check_user_input( length_bounds: (int | tuple[int] | list[int]), quality_threshold: int, verbose: bool, - output_filename: "str", + output_filename: str, ): """ Parse input fastq file to dictionary[sequence_name: [sequence, sequence_quality]] and check if input can be correctly processed\n From 93addf3a4189761000e18942b3536d6aa9627453 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Tue, 24 Oct 2023 02:20:44 +0300 Subject: [PATCH 22/27] Add bio_files_processor.py with 3 functions --- bio_files_processor.py | 133 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 bio_files_processor.py diff --git a/bio_files_processor.py b/bio_files_processor.py new file mode 100644 index 0000000..ded1ab1 --- /dev/null +++ b/bio_files_processor.py @@ -0,0 +1,133 @@ +import os + + +def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = ""): + """ + Converts a multiline FASTA file to a one-line FASTA file. + + Arguments: + - input_fasta (str): Path to the input FASTA file. + - output_fasta (str): Path to the output FASTA file.\\ + If output_fasta is not provided, output file will be saved in the current directory with the same name as input file but with "_oneline.fasta" suffix. + + Return: + None + + Example: + convert_multiline_fasta_to_oneline("input.fasta", "output.fasta")""" + if output_fasta == "": + output_fasta = os.path.basename(input_fasta) + output_fasta = output_fasta.replace(".fasta", "_oneline.fasta") + output_path = os.path.join("./", output_fasta) + with open(input_fasta, "r") as input: + with open(output_path, "w"): + read = [] + while True: + line = input.readline().strip() + if not line: + break + if line.startswith(">"): + line += "\n" + if read: + with open(output_path, "a") as output: + output.write("".join(read) + "\n") + read = [line] + else: + read.append(line) + with open(output_path, "a") as output: + output.write("".join(read)) + + +def select_genes_from_gbk_to_fasta( + input_gbk: str, + genes: tuple[str] or list[str], + n_before: int = 1, + n_after: int = 1, + output_fasta: str = "", +): + """ + Extract the translations of neighbour genes in amino acids from a GenBank (gbk) file for a list of desired genes and saves them in a FASTA file format ready for blasting. + + Arguments: + - input_gbk (str): Path to the input GenBank (gbk) file. + - genes (tuple[str] or list[str]): List of gene names to extract translations for. + - n_before (int, optional): Number of genes to include before the desired gene. Defaults to 1. + - n_after (int, optional): Number of genes to include after the desired gene. Defaults to 1. + - output_fasta (str, optional): Path to the output FASTA file.\\ + If output_fasts is not provided, output file will be saved in the current directory with the same name as input file but with "_trans_for_blast.fasta" suffix. + + Return: + None + + Example: + select_genes_from_gbk_to_fasta("input.gbk", ["gene1", "gene2"], n_before=2, n_after=2, output_fasta="output.fasta") + """ + if output_fasta == "": + output_fasta = os.path.basename(input_gbk) + output_fasta = output_fasta.replace(".gbk", "_trans_for_blast.fasta") + output_path = os.path.join("./", output_fasta) + with open(input_gbk, "r") as input: + translations = [] + ind = 0 + prev_gene = "undefined" + was_printed = set() + while True: + line = input.readline().strip() + if not line: + break + if line.startswith("/gene="): + prev_gene = line[7:-1] + if line.startswith("/translation="): + translation = line[14:] + while not translation.endswith('"'): + line = input.readline().strip() + translation += line + translation = translation.rstrip('"') + translations.append((prev_gene, translation)) + prev_gene = "undefined" + with open(output_path, "w") as out: + for gene in genes: + for ind, cur in enumerate(translations): + if gene in cur[0]: + for i in range( + max(0, ind - n_before), + min(len(translations), ind + n_after + 1), + ): + if i == ind or i in was_printed: + continue + was_printed.add(i) + out_gene, translation = translations[i] + out.write(f">{out_gene}\n{translation}\n") + + +def change_fasta_start_pos(input_fasta: str, shift: int, output_fasta: str = ""): + """ + Shift the starting position of each sequence in a fasta file by a specified number of nucleotides and saves the modified sequences in a new fasta file. + + Arguments: + - input_fasta (str): Path to the input fasta file. + - shift (int): Number of nucleotides to shift the starting position of each sequence. + - output_fasta (str, optional): Path to the output fasta file.\\ + If output_fasta is not provided, output file will be saved in the current directory with the same name as input file but with "_shifted.fasta" suffix. + + Return: + None + + Example: + change_fasta_start_pos("input.fasta", 3, "output.fasta") + """ + if output_fasta == "": + output_fasta = os.path.basename(input_fasta) + output_fasta = output_fasta.replace(".fasta", "_shifted.fasta") + output_path = os.path.join("./", output_fasta) + with open(input_fasta, "r") as input: + with open(output_path, "w"): + while True: + line = input.readline().strip() + print(line) + if not line: + break + if not line.startswith(">"): + line = line[shift:] + line[:shift] + with open(output_path, "a") as output: + output.write(line + "\n") From ca5abedfa36574de1c4454fbac859c056ed5d0d5 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Tue, 24 Oct 2023 02:26:57 +0300 Subject: [PATCH 23/27] Update bio_files_processor.py with 3 --- bio_files_processor.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/bio_files_processor.py b/bio_files_processor.py index ded1ab1..cba3ad7 100644 --- a/bio_files_processor.py +++ b/bio_files_processor.py @@ -11,10 +11,11 @@ def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = "") If output_fasta is not provided, output file will be saved in the current directory with the same name as input file but with "_oneline.fasta" suffix. Return: - None + None Example: - convert_multiline_fasta_to_oneline("input.fasta", "output.fasta")""" + convert_multiline_fasta_to_oneline("input.fasta", "output.fasta") + """ if output_fasta == "": output_fasta = os.path.basename(input_fasta) output_fasta = output_fasta.replace(".fasta", "_oneline.fasta") From 313ed97d1137ece7ac6fc0d0122d5bfde2fab675 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Tue, 24 Oct 2023 02:27:35 +0300 Subject: [PATCH 24/27] Update fastq_filter.py --- src/fastq_filter.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/fastq_filter.py b/src/fastq_filter.py index 7e0e7db..28becd3 100644 --- a/src/fastq_filter.py +++ b/src/fastq_filter.py @@ -169,6 +169,9 @@ def save_filtered_seqs( - filtered_seqs (dict[str, list[str]]): filtered fastq_reads - output_filename (str): output name of the filtered fastq file - save_to_dir (str): absolute or realtive path to directory to save to + + Return: + None """ os.makedirs(save_to_dir, exist_ok=True) output_path = os.path.join(save_to_dir, output_filename) From f8bcddf78ad83dffdafa7f2017ce512a31ddc146 Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Tue, 24 Oct 2023 02:35:17 +0300 Subject: [PATCH 25/27] Update bio_files_processor.py --- bio_files_processor.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/bio_files_processor.py b/bio_files_processor.py index cba3ad7..5f1b11f 100644 --- a/bio_files_processor.py +++ b/bio_files_processor.py @@ -41,7 +41,7 @@ def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = "") def select_genes_from_gbk_to_fasta( input_gbk: str, - genes: tuple[str] or list[str], + genes: str or tuple[str] or list[str], n_before: int = 1, n_after: int = 1, output_fasta: str = "", @@ -51,7 +51,7 @@ def select_genes_from_gbk_to_fasta( Arguments: - input_gbk (str): Path to the input GenBank (gbk) file. - - genes (tuple[str] or list[str]): List of gene names to extract translations for. + - genes (str or tuple[str] or list[str]): Gene name or list of gene names to extract neigbhour genes translations for. - n_before (int, optional): Number of genes to include before the desired gene. Defaults to 1. - n_after (int, optional): Number of genes to include after the desired gene. Defaults to 1. - output_fasta (str, optional): Path to the output FASTA file.\\ @@ -67,6 +67,8 @@ def select_genes_from_gbk_to_fasta( output_fasta = os.path.basename(input_gbk) output_fasta = output_fasta.replace(".gbk", "_trans_for_blast.fasta") output_path = os.path.join("./", output_fasta) + if isinstance(genes, str): + genes = [genes] with open(input_gbk, "r") as input: translations = [] ind = 0 @@ -101,6 +103,14 @@ def select_genes_from_gbk_to_fasta( out.write(f">{out_gene}\n{translation}\n") +select_genes_from_gbk_to_fasta( + "/mnt/c/users/vovag/Documents/IB/Python_1sem/HW5_Grigoriants/example_gbk.gbk", + genes="ybgD_2", + n_before=10, + n_after=10, +) + + def change_fasta_start_pos(input_fasta: str, shift: int, output_fasta: str = ""): """ Shift the starting position of each sequence in a fasta file by a specified number of nucleotides and saves the modified sequences in a new fasta file. From 49756b1eefa078b75f9b065f3dddc5571f53b7fe Mon Sep 17 00:00:00 2001 From: VovaGrig Date: Tue, 24 Oct 2023 02:46:01 +0300 Subject: [PATCH 26/27] Update bio_files_processor.py --- bio_files_processor.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/bio_files_processor.py b/bio_files_processor.py index 5f1b11f..5ba06dd 100644 --- a/bio_files_processor.py +++ b/bio_files_processor.py @@ -103,14 +103,6 @@ def select_genes_from_gbk_to_fasta( out.write(f">{out_gene}\n{translation}\n") -select_genes_from_gbk_to_fasta( - "/mnt/c/users/vovag/Documents/IB/Python_1sem/HW5_Grigoriants/example_gbk.gbk", - genes="ybgD_2", - n_before=10, - n_after=10, -) - - def change_fasta_start_pos(input_fasta: str, shift: int, output_fasta: str = ""): """ Shift the starting position of each sequence in a fasta file by a specified number of nucleotides and saves the modified sequences in a new fasta file. From 21dce21d789bcba7a196dfa6ac304a413afefbe4 Mon Sep 17 00:00:00 2001 From: Vladimir Grigoriants Date: Tue, 24 Oct 2023 04:00:04 +0400 Subject: [PATCH 27/27] Update README.md --- README.md | 45 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 43 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 104cd4c..bbf21e9 100644 --- a/README.md +++ b/README.md @@ -2,10 +2,12 @@ ## Tools to work with nucleic acids, protein sequences, fastq reads. Every bioinfomatician at some point of his/her career deals with nucleic acid sequences, protein sequences. Nowadays it's also hard to imagine bioinformatics without New Generation Sequecning methods and data, that needs to be analyzed with bioinformatician hands. For those reasons as well as for practical experience and better theorethical understanding of bioinformatics I create this repo. -`bioinf_tools.py` is an open-source program that facilitates working with bioinformatics data. +`bioinf_tools.py` and 'bio_files_processor.py' are an open-source program that facilitate working with bioinformatics data. ## Usage and options -The programm is based on three main functions: + +### Bioinf_tools.py +`bioinf_tools.py` program contains three main functions: **1.** `run_protein_tools` - facilitates work with protein sequences, has 5 procedures: @@ -75,6 +77,33 @@ For *gc_bounds* and *length_bounds* if one number is provided bounds from 0 to n Please provide threshold quality in Phred-33 scale. +### Bio Files Processor + +`bio_files_processor.py` program contains functions for processing biological sequence files in FASTA and GenBank (gbk) formats. + +#### Functions: + +- convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = ""): Converts a multiline FASTA file to a one-line FASTA file. + - Arguments: + - input_fasta (str): Path to the input FASTA file. + - output_fasta (str): Path to the output FASTA file. If not provided, output file will be saved in the current directory with the same name as input file but with "_oneline.fasta" suffix. + - Returns: None + +- select_genes_from_gbk_to_fasta(input_gbk: str, genes: str or tuple[str] or list[str], n_before: int = 1, n_after: int = 1, output_fasta: str = ""): Extracts the translations of neighbour genes in amino acids from a GenBank (gbk) file for a list of desired genes and saves them in a FASTA file format ready for blasting. + - Arguments: + - input_gbk (str): Path to the input GenBank (gbk) file. + - genes (str or tuple[str] or list[str]): Gene name or list of gene names to extract neigbhour genes translations for. + - n_before (int, optional): Number of genes to include before the desired gene. Defaults to 1. + - n_after (int, optional): Number of genes to include after the desired gene. Defaults to 1. + - output_fasta (str, optional): Path to the output FASTA file. If not provided, output file will be saved in the current directory with the same name as input file but with "_trans_for_blast.fasta" suffix. + - Returns: None + +- change_fasta_start_pos(input_fasta: str, shift: int, output_fasta: str = ""): Shifts the starting position of each sequence in a FASTA file by a specified number of nucleotides and saves the modified sequences in a new FASTA file. + - Arguments: + - input_fasta (str): Path to the input FASTA file. + - shift (int): Number of nucleotides to shift the starting position of each sequence. + - output_fasta (str, optional): Path to the output FASTA file. If not provided, output file will be saved in the current directory with the same name as input file but with "_shifted.fasta" suffix. + - Returns: None ## Examples ```python ### run_protein_tools @@ -132,4 +161,16 @@ run_fastq_filter(gc_bounds=50) #{'@SRX079804:1:SRR292678:1:1101:21885:21885': ('ACAGCAACATAAACATGATGGGATGGCGTAAGCCCCCGAGATATCAGTTTACCCAGGATAAGAGATTAAATTATGAGCAACATTATTAA', 'FGGGFGGGFGGGFGDFGCEBB@CCDFDDFFFFBFFGFGEFDFFFF;D@DD>C@DDGGGDFGDGG?GFGFEGFGGEF@FDGGGFGFBGGD'), # '@SRX079804:1:SRR292678:1:1101:30161:30161': ('GAACGACAGCAGCTCCTGCATAACCGCGTCCTTCTTCTTTAGCGTTGTGCAAAGCATGTTTTGTATTACGGGCATCTCGAGCGAATC', 'DFFFEGDGGGGFGGEDCCDCEFFFFCCCCCB>CEBFGFBGGG?DE=:6@=>AD?D8DCEE:>EEABE5D@5:DDCA;EEE-DCD')} + +### convert_multiline_fasta_to_oneline +convert_multiline_fasta_to_oneline("input.fasta", "output.fasta") + + +### select_genes_from_gbk_to_fasta +select_genes_from_gbk_to_fasta(example_gbk.gbk, genes="ybgD_2", n_before=10,n_after=10) + + +### change_fasta_start_pos +change_fasta_start_pos("input.fasta", 3, "output.fasta") ``` +