From 91e416f641f8bb68fd6f1caed9e80d52975a2445 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sat, 7 Oct 2023 23:46:41 +0300 Subject: [PATCH 01/30] Innitial commit for dna_rna_analysis.py --- modules_for_BSAT/dna_rna_analysis.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 modules_for_BSAT/dna_rna_analysis.py diff --git a/modules_for_BSAT/dna_rna_analysis.py b/modules_for_BSAT/dna_rna_analysis.py new file mode 100644 index 0000000..e69de29 From 52772ece7db27afbc75eecc7b617dbfad3afc266 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sat, 7 Oct 2023 23:47:29 +0300 Subject: [PATCH 02/30] Innitial commit for fastq_analysis.py --- modules_for_BSAT/fastq_analysis.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 modules_for_BSAT/fastq_analysis.py diff --git a/modules_for_BSAT/fastq_analysis.py b/modules_for_BSAT/fastq_analysis.py new file mode 100644 index 0000000..e69de29 From edf755beee4cc7cc97235be5dd856dc943478f27 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sat, 7 Oct 2023 23:48:42 +0300 Subject: [PATCH 03/30] Innitial commit for protein_analysis --- modules_for_BSAT/protein_analysis.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 modules_for_BSAT/protein_analysis.py diff --git a/modules_for_BSAT/protein_analysis.py b/modules_for_BSAT/protein_analysis.py new file mode 100644 index 0000000..e69de29 From b2471156463e35edc7edea6dbef0cd782403a3ac Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sat, 7 Oct 2023 23:52:07 +0300 Subject: [PATCH 04/30] Add python script with all functions for this module --- modules_for_BSAT/dna_rna_analysis.py | 96 ++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) diff --git a/modules_for_BSAT/dna_rna_analysis.py b/modules_for_BSAT/dna_rna_analysis.py index e69de29..1e32797 100644 --- a/modules_for_BSAT/dna_rna_analysis.py +++ b/modules_for_BSAT/dna_rna_analysis.py @@ -0,0 +1,96 @@ +def transcribe(seq: str) -> str: + """ + Function return return transcribed sequence. + + :param seq: DNA sequence + :type seq: str + :rtype: str + :return: transcribed sequence + """ + nucleotide_type = set(seq.upper()) + if 'U' in nucleotide_type: + raise ValueError(f'Sequence {seq} is not DNA!') + rna_seq = seq.replace('T', 'U').replace('t', 'u') + return rna_seq + + +def reverse(seq: str) -> str: + """ + Function return return reversed sequence. + + :param seq: DNA or RNA sequence + :type seq: str + :rtype: str + :return: reversed sequence + """ + reverse_seq = seq[::-1] + return reverse_seq + + +def complement(seq: str) -> str: + """ + Function return return complement sequence. + + :param seq: DNA or RNA sequence + :type seq: str + :rtype: str + :return: complement sequence + """ + complement_dict = {'A': 'T', 'C': 'G', + 'G': 'C', 'T': 'A', 'U': 'A', 'a': 't', + 'c': 'g', 'g': 'c', 't': 'a', 'u': 'a'} + complement_seq = [] + length = len(seq) + for i in range (length): + if seq[i] in complement_dict: + complement_seq.append(complement_dict[seq[i]]) + return ''.join(complement_seq) + + +def reverse_complement(seq: str) -> str: + """ + Function return return reverse complement sequence. + + :param seq: DNA or RNA sequence + :type seq: str + :rtype: str + :return: reverse complement sequence + """ + seq = complement(seq) + reverse_complement_seq = reverse(seq) + return reverse_complement_seq + + +def gc_calculate(seq: str) -> float: + """ + Function return sequence GC-content in percent. + + :param seq: DNA or RNA sequence + :type seq: str + :rtype: float + :return: GC-contentn percent + """ + length = len(seq) + gc_content = 0.0 + seq_up = seq.upper() + c = seq_up.count("C") + g = seq_up.count("G") + gc_content = round(((c+g)/length*100),2) + return gc_content + + +def is_na(seq: str) -> None: + """ + Function return None if sequence is DNA or RNA. + + :param seq: some sequence + :type seq: str + :rtype: None + :return: None + """ + nucleotide_set = {'A', 'G', 'C', 'T', 'U'} + seq_nucleotide_type = set(seq.upper()) + if seq_nucleotide_type.issubset(nucleotide_set) != True: + raise ValueError(f'Sequence {seq} is not nucleic acid!') + elif 'T' in seq_nucleotide_type and 'U' in seq_nucleotide_type: + raise ValueError(f'Sequence {seq} is not nucleic acid!') From 78dc3b59be41cfd80b97307e441aaf80886d7e5a Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sat, 7 Oct 2023 23:53:03 +0300 Subject: [PATCH 05/30] Add python script with all functions for this module --- modules_for_BSAT/fastq_analysis.py | 59 ++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/modules_for_BSAT/fastq_analysis.py b/modules_for_BSAT/fastq_analysis.py index e69de29..0dba952 100644 --- a/modules_for_BSAT/fastq_analysis.py +++ b/modules_for_BSAT/fastq_analysis.py @@ -0,0 +1,59 @@ +def is_nucleotide(seq: str) -> bool: + """ + Function return 'True' if sequence is DNA or RNA. + + :param seq: some sequence + :type seq: str + :rtype: bool + :return: 'True' if sequence is DNA or RNA + """ + seq_alphabet = {'A', 'T', 'U', 'G', 'C'} + if set(seq.upper()) <= seq_alphabet: + return True + + +def analyse_gc(seq: str) -> float: + """ + Return GC-content of DNA/RNA sequence. + + :param seq: DNA/RNA sequence + :type seq: str + :rtype: float + :return: User sequence GC-content + """ + length = len(seq) + gc_content = 0.0 + seq_up = seq.upper() + c = seq_up.count("C") + g = seq_up.count("G") + gc_content = round(((c+g)/length*100),2) + return gc_content + + +def analyse_length(seq: str) -> int: + """ + Return length of DNA/RNA sequence. + + :param seq: DNA/RNA sequence + :type seq: str + :rtype: int + :return: User sequence length + """ + length = len(seq) + return length + + +def analyse_quality(quality: str) -> float: + """ + Return quality score of read, that coding by ASCII code. + + :param seq: quality symbols for each nucleotide + :type seq: str + :rtype: float + :return: User sequence quality + """ + q_score = 0 + for char in quality: + q_score += ord(char) - 33 + q_score_sum = q_score/len(quality) + return round(q_score_sum,2) From c6ade48290a29223cab12aa97b28bf0319a32f55 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sat, 7 Oct 2023 23:54:07 +0300 Subject: [PATCH 06/30] Add python script with all functions for protein module --- modules_for_BSAT/protein_analysis.py | 316 +++++++++++++++++++++++++++ 1 file changed, 316 insertions(+) diff --git a/modules_for_BSAT/protein_analysis.py b/modules_for_BSAT/protein_analysis.py index e69de29..3e6f98c 100644 --- a/modules_for_BSAT/protein_analysis.py +++ b/modules_for_BSAT/protein_analysis.py @@ -0,0 +1,316 @@ +from typing import List, Union + +RESIDUES_NAMES = {'ALA': 'A', + 'ARG': 'R', + 'ASN': 'N', + 'ASP': 'D', + 'CYS': 'C', + 'GLN': 'Q', + 'GLU': 'E', + 'GLY': 'G', + 'HIS': 'H', + 'ILE': 'I', + 'LEU': 'L', + 'LYS': 'K', + 'MET': 'M', + 'PHE': 'F', + 'PRO': 'P', + 'SER': 'S', + 'THR': 'T', + 'TRP': 'W', + 'TYR': 'Y', + 'VAL': 'V' + } + +RESIDUES_NAMES_THREE = {'A': 'ALA', + 'R': 'ARG', + 'N': 'ASN', + 'D': 'ASP', + 'C': 'CYS', + 'Q': 'GLN', + 'E': 'GLU', + 'G': 'GLY', + 'H': 'HIS', + 'I': 'ILE', + 'L': 'LEU', + 'K': 'LYS', + 'M': 'MET', + 'F': 'PHE', + 'P': 'PRO', + 'S': 'SER', + 'T': 'THR', + 'W': 'TRP', + 'Y': 'TYR', + 'V': 'VAL' + } + +# first value is hydrophobicity index, second is pKa (pKa1, pKa2, pKa3 respectively), third is molecular mass in Da +RESIDUES_CHARACTERISTICS = {'A': [1.8, [2.34, 9.69, 0], 89], + 'R': [-4.5, [2.17, 9.04, 12.48], 174], + 'N': [-3.5, [2.02, 8.80, 0], 132], + 'D': [-3.5, [1.88, 9.60, 3.65], 133], + 'C': [2.5, [1.96, 10.28, 8.18], 121], + 'Q': [-3.5, [2.17, 9.13, 0], 146], + 'E': [-3.5, [2.19, 9.67, 4.25], 147], + 'G': [-0.4, [2.34, 9.60, 0], 75], + 'H': [-3.2, [1.82, 9.17, 6.00], 155], + 'I': [4.5, [2.36, 9.60, 0], 131], + 'L': [3.8, [2.36, 9.60, 0], 131], + 'K': [-3.9, [2.18, 8.95, 10.53], 146], + 'M': [1.9, [2.28, 9.21, 0], 149], + 'F': [2.8, [1.83, 9.13, 0], 165], + 'P': [-1.6, [1.99, 10.60, 0], 115], + 'S': [-0.8, [2.21, 9.15, 0], 105], + 'T': [-0.7, [2.09, 9.10, 0], 119], + 'W': [-0.9, [2.83, 9.39, 0], 204], + 'Y': [-1.3, [2.20, 9.11, 0], 181], + 'V': [4.2, [2.32, 9.62, 0], 117]} + +# amino acid with corresponding degenerate codon/codons +AMINO_ACID_TO_MRNA = {'A': 'GCN', + 'R': '(CGN/AGR)', + 'N': 'AAY', + 'D': 'GAY', + 'C': 'UGY', + 'Q': 'CAR', + 'E': 'GAR', + 'G': 'GGN', + 'H': 'CAY', + 'I': 'AUH', + 'L': '(CUN/UUR)', + 'K': 'AAR', + 'M': 'AUG', + 'F': 'UUY', + 'P': 'CCN', + 'S': '(UCN/AGY)', + 'T': 'ACN', + 'W': 'UGG', + 'Y': 'UAY', + 'V': 'GUN'} + +def save_register(seq: str) -> List[int]: + """ + Additional function for correct change_residues_encoding return + :param seq: protein seq (str) + :return: List of type of register character in protein seq (List[num]) + """ + register = [] + sep_seq = [] + if ' ' in seq: + sep_seq = seq.split() + else: + sep_seq = list(seq) + for residue in sep_seq: + if residue.isupper(): + register.append(1) + else: + register.append(0) + return register + + +def change_residues_encoding(seq: str, encoding = 'one') -> str: + """ + Transfer amino acids from 3-letter to 1-letter code. By default, converts all seq into 1-letter + :param seq: protein seq (str) + :param encoding: specify target encoding (str) + :return: same protein seq in another encoding (str) + """ + encode_seq = [] + encode_seq_registered = [] + if ' ' in seq: + sep_seq = seq.upper().split() + else: + sep_seq = list(seq) + residue_code_length = len(sep_seq[0]) + for residue in sep_seq: + if len(residue) != residue_code_length: + raise ValueError (f'Wrong sequence format in {seq}!') + if encoding == 'one': + if len(sep_seq[0]) == 1: + encode_seq.append(residue) + if len(sep_seq[0]) == 3: + encode_seq.append(RESIDUES_NAMES[residue]) + if encoding == 'three': + if len(sep_seq[0]) == 3: + encode_seq.append(residue) + if len(sep_seq[0]) == 1: + encode_seq.append(RESIDUES_NAMES_THREE[residue]) + for residue, reg in zip(encode_seq, save_register(seq)): + if (reg == 1): + encode_seq_registered += residue.upper() + elif (reg == 0): + encode_seq_registered += residue.lower() + if encoding == 'one': + return ''.join(encode_seq_registered) + if encoding == 'three': + fin_seq = [] + for i, residue in enumerate(encode_seq_registered): + fin_seq.append(residue) + if ((i+1) % 3 == 0): + fin_seq.append(' ') + return ' '.join(fin_seq) + + +def is_protein(seq: str) -> bool: + """ + Check if sequence is protein or not by identify invalid seq elements, which are not presented in dicts above. + :param seq: protein seq in 1-letter encoding (str) + :return: if seq is correct protein seq or not (bool) + """ + for residue in seq.upper(): + if residue not in RESIDUES_NAMES.values(): + return False + return True + + +def get_seq_characteristic(seq: str) -> dict: + """ + Count entry of each residue type in your seq. Get description of amino acid composition. + + :param seq: protein seq in 1-letter encoding (str) + :return: each residue type in seq in 3-letter code and its amount in current seq (dict) + """ + seq = seq.upper() + residue_count = {} + for residue in seq: + residue_count[[tl_code for tl_code in RESIDUES_NAMES if RESIDUES_NAMES[tl_code] == residue][0]] = 0 + for residue in seq: + residue_count[[tl_code for tl_code in RESIDUES_NAMES if RESIDUES_NAMES[tl_code] == residue][0]] += 1 + return residue_count + + +def find_residue(seq: str, residue_of_interest: str) -> str: + """ + Find all positions of certain residue in your seq. + + :param seq: protein seq in 1-letter encoding (str) + :param residue_of_interest: specify the residue of interest (str) + :return: positions of specified residue in your seq (str) + """ + residue_of_interest = residue_of_interest.upper() + seq = seq.upper() + if len(residue_of_interest) == 3: + residue_of_interest = RESIDUES_NAMES[residue_of_interest] + residue_of_interest_position = [] + for ind, residue in enumerate(seq, 1): + if residue == residue_of_interest: + residue_of_interest_position.append(ind) + return f'{residue_of_interest} positions: {residue_of_interest_position}' + + +def find_site(seq: str, site: str) -> str: + """ + Find if seq contains certain site and get positions of its site. + + :param seq: protein seq in 1-letter encoding (str) + :param site: specify site of interest (str) + :return: positions of residues for each certain site in seq (str) + """ + site = change_residues_encoding(site).upper() + seq = seq.upper() + if not is_protein(site): + return f'Site {site} is not a protein!' + if site in seq: + site_full_position = [] + site_count = seq.count(site) + site_start_position = [(coordinate + 1) for coordinate in range(len(seq)) if seq.startswith(site, coordinate)] + site_end_position = [(coordinate + len(site)) for coordinate in site_start_position] + for counter in range(len(site_start_position)): + site_full_position.append(f'{site_start_position[counter]}:{site_end_position[counter]}') + return f'Site entry in sequence = {site_count}. Site residues can be found at positions: {site_full_position}' + else: + return f'{site} site is not in sequence!' + + +def calculate_protein_mass(seq: str) -> float: + """ + Get mass of residues in your seq in Da. + + :param seq: protein seq in 1-letter encoding (str) + :return: mass in Da (float) + """ + total_mass = 0 + for res in seq.upper(): + total_mass += RESIDUES_CHARACTERISTICS[res][2] + return total_mass + + +def calculate_average_hydrophobicity(seq: str) -> float: + """ + Get hydrophobicity index for protein seq as sum of index for each residue + in your seq divided by its length. + + :param seq: protein seq in 1-letter encoding (str) + :return: average hydrophobicity (float) + """ + sum_hydrophobicity_ind = 0 + for res in seq.upper(): + sum_hydrophobicity_ind += RESIDUES_CHARACTERISTICS[res][0] + average_hydrophobicity = round(sum_hydrophobicity_ind / len(seq),2) + return average_hydrophobicity + + +def get_mrna(seq: str) -> List[str]: + """ + Get encoding mRNA nucleotides for your seq. + + :param seq: protein seq in 1-letter encoding (str) + :return: potential encoding mRNA sequence with multiple choice for some positions (str) + """ + mrna_seq = [] + for residue in seq.upper(): + mrna_seq.append(AMINO_ACID_TO_MRNA[residue]) + return mrna_seq + + +def calculate_isoelectric_point(seq: str) -> float: + """ + Find isoelectrinc point as sum of known pI for residues in your seq. + + :param seq: protein seq in 1-letter encoding (str) + :return: isoelectric point (float) + """ + sum_pka = 0 + pka_amount = 0 + for ind, res in enumerate(seq.upper(), 1): + if ind == 1: + sum_pka += RESIDUES_CHARACTERISTICS[res][1][1] + pka_amount += 1 + elif RESIDUES_CHARACTERISTICS[res][1][2] != 0: + sum_pka += RESIDUES_CHARACTERISTICS[res][1][2] + pka_amount += 1 + elif ind == len(seq): + sum_pka += RESIDUES_CHARACTERISTICS[res][1][0] + pka_amount += 1 + pi = round(sum_pka / pka_amount, 2) + return pi + + +def analyze_secondary_structure(seq: str) -> List[float]: + """ + Calculate the percentage of amino acids found in the three main + types of protein secondary structure: beta-turn, beta-sheet and alpha-helix. + + :param seq: protein seq in 1-letter encoding (str) + :return: percentage of amino acids belonging to three types of secondary structure (list[str]) + """ + b_turn_set = {'G', 'P', 'N', 'D'} + b_sheet_set = {'F', 'Y', 'I', 'V', 'C', 'W'} + alpha_helix_set = {'M', 'A', 'L', 'E', 'K'} + protein_length = len(seq) + result = [] + counter_b_turn = 0 + counter_b_sheet = 0 + counter_alpha_helix = 0 + for residue in seq.upper(): + if residue in b_turn_set: + counter_b_turn += 1 + elif residue in b_sheet_set: + counter_b_sheet += 1 + elif residue in alpha_helix_set: + counter_alpha_helix += 1 + result.append(round(counter_b_turn / protein_length * 100,2)) + result.append(round(counter_b_sheet / protein_length * 100,2)) + result.append(round(counter_alpha_helix / protein_length * 100,2)) + return result + From caee4ae05ae01730fa84ce33e14a24181353f425 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sat, 7 Oct 2023 23:57:03 +0300 Subject: [PATCH 07/30] Initail commit for main script of Bio_Seq_Analysis_Tool --- Bio_Seq_Analysis_Tool.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 Bio_Seq_Analysis_Tool.py diff --git a/Bio_Seq_Analysis_Tool.py b/Bio_Seq_Analysis_Tool.py new file mode 100644 index 0000000..e69de29 From 10e061b28591d58f9a0cc2f4594fbca6d68b07d9 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sat, 7 Oct 2023 23:58:42 +0300 Subject: [PATCH 08/30] Add forder with required modules for Bio_Seq_Analysis_Tool.py --- .../.ipynb_checkpoints/dna_rna_analysis-checkpoint.py | 0 modules_for_BSAT/.ipynb_checkpoints/fastq_analysis-checkpoint.py | 0 .../.ipynb_checkpoints/protein_analysis-checkpoint.py | 0 3 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 modules_for_BSAT/.ipynb_checkpoints/dna_rna_analysis-checkpoint.py create mode 100644 modules_for_BSAT/.ipynb_checkpoints/fastq_analysis-checkpoint.py create mode 100644 modules_for_BSAT/.ipynb_checkpoints/protein_analysis-checkpoint.py diff --git a/modules_for_BSAT/.ipynb_checkpoints/dna_rna_analysis-checkpoint.py b/modules_for_BSAT/.ipynb_checkpoints/dna_rna_analysis-checkpoint.py new file mode 100644 index 0000000..e69de29 diff --git a/modules_for_BSAT/.ipynb_checkpoints/fastq_analysis-checkpoint.py b/modules_for_BSAT/.ipynb_checkpoints/fastq_analysis-checkpoint.py new file mode 100644 index 0000000..e69de29 diff --git a/modules_for_BSAT/.ipynb_checkpoints/protein_analysis-checkpoint.py b/modules_for_BSAT/.ipynb_checkpoints/protein_analysis-checkpoint.py new file mode 100644 index 0000000..e69de29 From 1de06fd04ae1f46e709b03ad5e45bdbeb3729f1a Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sun, 8 Oct 2023 00:01:14 +0300 Subject: [PATCH 09/30] Add python script with all functions into Bio_Seq_Analysis_Tool --- Bio_Seq_Analysis_Tool.py | 159 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 159 insertions(+) diff --git a/Bio_Seq_Analysis_Tool.py b/Bio_Seq_Analysis_Tool.py index e69de29..285c1b7 100644 --- a/Bio_Seq_Analysis_Tool.py +++ b/Bio_Seq_Analysis_Tool.py @@ -0,0 +1,159 @@ +from typing import Union, Tuple, List, Set, Dict +from modules_for_BSAT import fastq_analysis as fq +from modules_for_BSAT import dna_rna_analysis as na +from modules_for_BSAT import protein_analysis as pa + + +def dna_rna_analysis(*args: str, operation: str) -> Union[List[float], List[str]]: + """ + This function performs a number of operations on DNA or RNA. + + Operations supported by this functions: + -transcribe - return transcribed sequence + -reverse - return reverse sequence + -complement - return complement sequence + -reverse_complement - return reverse complement sequence + -gc_calculate - return sequence GC-content in percent + + :param args: nucleic acid sequence + :type param seqs: str + + :param operation: type of operation required + :type param operation: str + + :return: Analysis of nucleic acid sequence + :rtype: List[str] + + :raises ValueError: if sequence not RNA or DNA, also if the operation value out of OPERATION_DICT + """ + OPERATION_DICT = {'transcribe': na.transcribe, + 'reverse': na.reverse, + 'reverse_complement': na.reverse_complement, + 'complement': na.complement, + 'gc_calculate': na.gc_calculate} + analysis = [] + for seq in args: + na.is_na(seq) + if operation in OPERATION_DICT.keys(): + analysis.append(OPERATION_DICT[operation](seq)) + else: + raise ValueError(f'Wrong operation!') + return analysis + + +def analyse_fastq(seqs: dict, + gc_bounds: Union[int, float, Tuple [int], Tuple [float]] = (0, 100), + length_bounds: Union[int, Tuple [int]] = (0, 2**32), + quality_threshold: float = 0.0) -> Dict[str,str]: + """ + This function help analyze a set of reads obtained from next-generation sequencing. + + The function allow to filter the desired reads according to three parameters: + GC-content, length and reading quality. + + :param seqs: + A dictionary consisting of fastq sequences. The structure is as follows: Key - string, sequence name. + The value is a tuple of two strings: sequence and quality. The sequence is RNA or DNA. + :type seqs: Dict[str] + + :param gc_bounds: + Boundary parameters for filtering sequences by GC-content. Save only reads with a GC-content between boundaries + or lower than one boundary. Lower boundary cannot be less than 0 and upper boundary cannot be greater than 100. + gc_bounds default value is (0,100). + :type param gc_bounds: Union[int, float, Tuple [int], Tuple [float]] + + :param length_bounds: + Boundary parameters for filtering sequences by length. Works the same as gc_bounds. Lower boundary cannot be less + than 0 and upper boundary cannot be greater than 2^32. length_bounds default value is (0,2^32) + :type param length_bounds: Union[int, Tuple [int] + + :param quality_threshold: + Threshold for quality of each nucleotide in read. Quality incodes by ASCII codes. The threshold cannot be more + than 40. quality_threshold default value is 0 + :type param quality_threshold: float + + :return: + New dictionaries with fastq sequence.The first one consisting of filtered fastq sequences and the other one with + sequences that did not pass filters. + :rtype: Dict[str] + + :raises ValueError: if sequence not RNA or DNA, also if the argument values are outside the allowed ones + """ + if type(gc_bounds) == float or type(gc_bounds) == int: + gc_bounds = (0,gc_bounds) + if gc_bounds[0] < 0 or gc_bounds[1] > 100: + raise ValueError(f'Wrong boundaries!') + if type(length_bounds) == int: + length_bounds = (0,length_bounds) + if length_bounds[0] < 0 or length_bounds[1] > 2**32: + raise ValueError(f'Wrong boundaries!') + if quality_threshold > 40: + raise ValueError(f'Wrong quality threshold!') + analysed_seq = {} + error_seq = {} + for seq in seqs.items(): + if fq.is_nucleotide(seq[1][0]) != True: + raise TypeError(f'Wrong sequence format') + if fq.analyse_gc(seq[1][0]) > gc_bounds[0] and fq.analyse_gc(seq[1][0]) < gc_bounds[1]: + if fq.analyse_length(seq[1][0]) > length_bounds[0] and fq.analyse_length(seq[1][0]) < length_bounds[1]: + if fq.analyse_quality(seq[1][1]) > quality_threshold: + analysed_seq[seq[0]] = (seq[1]) + else: + error_seq[seq[0]] = (seq[1]) + else: + error_seq[seq[0]] = (seq[1]) + else: + error_seq[seq[0]] = (seq[1]) + return analysed_seq, error_seq + + +def run_protein_analysis(*args: str) -> Union[List[str], List[float]]: + """ + Launch desired operation with proteins sequences. Pass comma-separated sequences, + additional argument (if certain function requires it) and specify function name you want to apply to all + sequences. + Pass arguments strictly in this order, otherwise it won't be parsed. + If the input is a sequence of amino acids written in a three-letter code, then the amino acids must be separated + by a space. If the input is a sequence of amino acids written in a single-letter code, then the amino acids may + not be separated by a space. + + :param args: + - seq (str): amino acids sequences for analysis in 1-letter or 3-letter code (as many as you wish) + - additional arg (str): necessary parameter for certain functions (for example, specify target protein site) + - operation name (str): specify procedure you want to apply + :type param args: str + + :return: the result of procedure in list or str format + :rtype: Union[List[str], List[float] + :raises ValueError: if sequence not protein or sequence in wrong format + """ + # first value is function name, second is real function, third is number of function arguments + function_names = {'change_residues_encoding': [pa.change_residues_encoding, 2], + 'is_protein': [pa.is_protein, 1], + 'get_seq_characteristic': [pa.get_seq_characteristic, 1], + 'find_residue': [pa.find_residue, 2], + 'find_site': [pa.find_site, 2], + 'calculate_protein_mass': [pa.calculate_protein_mass, 1], + 'calculate_average_hydrophobicity': [pa.calculate_average_hydrophobicity, 1], + 'get_mrna': [pa.get_mrna, 1], + 'calculate_isoelectric_point': [pa.calculate_isoelectric_point, 1], + 'analyze_secondary_structure': [pa.analyze_secondary_structure, 1]} + + procedure = args[-1] + + processed_result = [] + + seqs = [pa.change_residues_encoding(seq) for seq in args[:-1 * (function_names[procedure][1])]] + for idx, seq in enumerate(seqs): + if not pa.is_protein(seq): + raise ValueError(f'Sequence {seq} is not protein!') + continue + if function_names[procedure][1] == 1: + processed_result.append(function_names[procedure][0](seq)) + elif function_names[procedure][1] == 2: + add_arg = args[-2] + processed_result.append(function_names[procedure][0](seq, add_arg)) + if len(processed_result) == 1: + return processed_result[0] + return processed_result + From be10c5763210c75935396ce94fc857f4f027fc1a Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Sun, 8 Oct 2023 00:04:23 +0300 Subject: [PATCH 10/30] Add README for Bio_Seq_Analysis_Tool module with detailed description of this one --- README.md | 487 +++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 485 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 602c1e1..4c49bb8 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,485 @@ -# BSAT -BSAT - Biological Sequences Analysis Toolbox. Repo contains tools which can help work with nucleic acid or protein sequences and with NGS-reads +# BSAT - Biological Sequences Analysis Toolbox + +This repository contains tools which helps you work with nucleic acid or protein sequences and with NGS-reads. It is capable to process multiple sequences, that makes analysis faster. + +## Installation + +To use this toolbox one need to clone repository + +```shell +git clone git@github.com:grishchenkoira/BSAT.git +cd BSAT +``` + +### System requirements: + +Key packages and programs: +- [Python](https://www.python.org/downloads/) (version >= 3.9) + +## Usage + +```python +# import main function +from Bio_Seq_Analysis_Tool import Bio_Seq_Analysis_Tool +``` + +## Works with main functions + +This section contains description of Bio_Seq_Analysis_Tool main functions. + +### dna_rna_analysis(*args: str, operation: str) + +This function performs a number of operations on DNA or RNA. +Operations supported by this functions: +- transcribe - return transcribed sequence +- reverse - return reverse sequence +- complement - return complement sequence +- reverse_complement - return reverse complement sequence +- gc_calculate - return sequence GC-content in percent + +**Parameters:** +- **args**: *str* + +Nucleic acid sequence +- **operation**: *str* + +Type of operation required + +**Returns**: +- **analysis**: *str* + +Analysis of nucleic acid sequence + +### analyse_fastq(seqs, gc_bounds, length_bounds, quality_threshold) + +Apply one of the operations described below to fastq sequences. + +**Parameters:** +- **seqs**: *dict* + +A dictionary consisting of fastq sequences. The structure is as follows: Key - string, sequence name. The value is a tuple of two strings: sequence and quality. The sequence is RNA or DNA. + +- **gc_bounds**: *Union[int, float, Tuple [int], Tuple [float]]* + +Boundary parameters for filtering sequences by GC-content. Save only reads with a GC-content between boundaries or lower than one boundary. Lower boundary cannot be less than 0 and upper boundary cannot be greater than 100. gc_bounds default value is (0,100). + +- **length_bounds** : *Union[int, Tuple [int]]* + +Boundary parameters for filtering sequences by length. Works the same as gc_bounds. Lower boundary cannot be less than 0 and upper boundary cannot be greater than 2^32. length_bounds default value is (0,2^32). + +- **quality_threshold** : *float* + +Threshold for quality of each nucleotide in read. Quality incodes by ASCII codes. The threshold cannot be more than 40. quality_threshold default value is 0 + +**Returns**: +- **analysed_seq**: *Dict[str]* + +New dictionary with fastq sequence.This one consists of filtered fastq sequences and + +- **analysed_seq**: *Dict[str, str]* + +New dictionary with fastq sequence. This one consists of sequences that did not pass filters. + +### run_protein_analysis(*args: str) + +Apply operations described below to any number of sequences with any case. + +**Parameters:** +**\*args**: +- **sequences**: *str* + +input coma-separated sequences in 1-letter or 3-letter code with any case (as many as you wish) +- **add_arg**: *str* + +necessary parameter for certain functions (for example, specify target protein site) +- **procedure** : *str* + +specify procedure you want to apply + +**Returns**: +- **operation_result**: str or list + +result of function work in list or str format (dependent on number of input sequences) + +**Note!** +- Operation name always must be the last argument +- Additional argument must be always before operation name + +## Modules + +This section contains description of modules using by main functions you can find in our library. + +- DNA & RNA analysis tool(#title1) +- FASTQ analysis tool(#title2) +- Amino acid sequences analysis tool(#title3) + +### DNA & RNA analysis + +This module performs a number of operations on DNA or RNA. + +#### Operations + +##### transcribe(seq) + +Function return return transcribed sequence. + +**Parameters:** +- **seq**: *str* + +DNA sequence + +**Returns:** +- **rna_seq**: *str* + +**Example** +```python +run_dna_rna_tools('ATG', 'transcribe') # 'AUG' +``` + + +##### reverse(seq) + +Function return return reversed sequence. + +**Parameters:** +- **seq**: *str* + +DNA or RNA sequence + +**Returns:** +- **reverse_seq**: *str* + +**Example** +```python +run_dna_rna_tools('ATG', 'reverse') # 'GTA' +``` + +##### complement(seq) + +Function return return complement sequence. + +**Parameters:** +- **seq**: *str* + +DNA or RNA sequence + +**Returns:** +- **complement_seq**: *str* + +**Example** +```python +run_dna_rna_tools('AtG', 'complement') # 'TaC' +``` + +##### reverse_complement(seq) + +Function return return reverse complement sequence. + +**Parameters:** +- **seq**: *str* + +DNA or RNA sequence + +**Returns:** +- **reverse_complement_seq**: *str* + +**Example** +```python +run_dna_rna_tools('ATg', 'reverse_complement') # 'cAT' +``` + +##### gc_calculate(seq) + +Function return sequence GC-content in percent. + +**Parameters:** +- **seq**: *str* + +DNA or RNA sequence + +**Returns:** +- **gc_content**: *str* + +**Example** +```python +run_dna_rna_tools ('GTAccca','gc_calculate') # '57.14' +``` + + +### FASTQ analysis tool + +This module contains functions for FASTQ sequnces filtration. The function allow to filter the desired reads according to three parameters: GC-content, length and reading quality. + +#### Operations + +##### analyse_gc(seq) + +Return GC-content of DNA/RNA sequence. + +**Parameters:** + +- **seq**: *str* + +DNA/RNA sequence + +**Returns:** +- **gc_content**: *float* + +##### analyse_length(seq) + +Return length of DNA/RNA sequence + +**Parameters:** + +- **seq**: *str* + +DNA/RNA sequence + +**Returns:** +- **length**: *int* + +##### analyse_quality(seq) + +Return quality score of read, that coding by ASCII code + +**Parameters:** + +- **seq**: *str* + +quality symbols for each nucleotide + +**Returns:** +- **q_score_sum**: *float* + + +### Amino acid sequences analysis tool + +This module contains functions for protein sequences analysis. You can reencode peptides sequences: 1-letter to 3-letter code and vice versa, calculate physical features, find specific sites, get predicted mRNA that coding your protein. + +#### Operations + +##### change_residues_encoding(seq, query='one') + +Transfer amino acids from 3-letter to 1-letter code and vice versa. + +**Parameters:** + +- **seq**: *str* + +Input protein seq in any encoding and case. If the input is a sequence of amino acids written in a three-letter code, then the amino acids must be separated by a space. If the input is a sequence of amino acids written in a single-letter code, then the amino acids may not be separated by a space. + +- **encoding**: {'one', 'three'}, default: 'one' + +specify target encoding + +**Returns:** +- **encode_seq_registered**: *str* + +same protein seq in another encoding + +**Example** +```python +seq = 'AAA' +change_residues_encoding(seq, 'one', 'change_residues_encoding') # 'AAA' + +seq = 'ALA ALA ALA' +change_residues_encoding(seq, 'one', 'change_residues_encoding') # 'AAA' + +seq = 'AAA' +change_residues_encoding(seq, 'three', 'change_residues_encoding') # 'ALA ALA ALA' +``` + +##### is_protein(seq) + +Check if sequence is protein or not by identify invalid seq elements, which are not presented in dicts above. + +**Parameters:** +- **seq**: *str* + +input protein seq in 1-letter encoding and upper case + +**Returns:** +- **verification_result**: *bool* + +if seq is correct protein seq or not + +**Example** +```python +seq = 'AAA' +is_protein(seq) #True +``` + +##### get_seq_characteristic(seq) + +Count entry of each residue type in your sequence + +**Parameters:** +- **seq**: *str* + +input protein seq in 1-letter encoding and upper case + +**Returns:** +- **res_count**: *dict* + +each residue type in seq in 3-letter code and its amount in current seq + +**Example** +```python +seq = 'AAA' +get_seq_characteristic(seq) #{'ALA': 3} +``` + +##### find_res(seq, res_of_interest) + +Find all positions of certain residue in your seq + +**Parameters:** +- **seq**: *str* + +input protein seq in 1-letter encoding and upper case +- **res_of_interest**: *str* + +residue of interest in 1-letter encoding and upper case + +**Returns:** +- **res_positions**: *str* + +positions of specified residue in your seq + +**Example** +```python +seq = 'AAA' +res = 'A' +find_res(seq, res) # 'A positions: [1, 2, 3]' +``` + +##### find_site(seq, site) + +Find if seq contains certain site and get positions of its site + +**Parameters:** +- **seq**: *str* + +input protein seq in 1-letter encoding and upper case + +- **site**: *str* + +specify site of interest + +**Returns:** +- **site_positions**: *str* + +the range of values for amino acid positions of specified site in your seq in which the last number is excluded + +**Example** +```python +seq = 'AAADDDF' +site = 'AAA' +find_site(seq, site) # "Site entry in sequence = 1. Site residues can be found at positions: ['1:4']" +``` + +##### calculate_protein_mass(seq) + +Get sum of residues masses in your seq in Da + +**Parameters:** +- **seq**: *str* + +input protein seq in 1-letter encoding and upper case + +**Returns:** +- **total_mass**: *float* + +mass of all residues in seq in Da + +**Example** +```python +seq = 'AAA' +calculate_protein_mass(seq) #267 +``` + +##### calculate_average_hydrophobicity(seq) + +Get average hydrophobicity index for protein seq as sum of index for each residue in your seq divided by its length + +**Parameters:** +- **seq**: *str* + +input protein seq in 1-letter encoding and upper case + +**Returns:** +- **average_hydrophobicity_idx**: *float* + +average hydrophobicity index for your seq + +**Example** +```python +seq = 'AAA' +calculate_average_hydrophobicity(seq) #1.8 +``` + +##### get_mrna(seq) + +Get encoding mRNA nucleotides for your seq + +**Parameters:** +- **seq**: *str* + +input protein seq in 1-letter encoding and upper case + +**Returns:** +- **mrna_seq**: *str* + +potential encoding mRNA sequences with multiple choice for some positions + +**Example** +```python +seq = 'AAA' +get_mrna(seq) # ['GCN', 'GCN', 'GCN'] +``` + +##### calculate_isoelectric_point(seq) + +Find isoelectrinc point as sum of known pI for residues in your seq + +**Parameters:** +- **seq**: *str* + +input protein seq in 1-letter encoding and upper case + +**Returns:** +- **pi**: *float* + +isoelectric point for your seq + +**Example** +```python +seq = 'AAA' +calculate_isoelectric_point(seq) # 6.01 +``` +##### analyze_secondary_structure(seq) + +Calculates the percentage of amino acids found in the three main types of protein secondary structure: beta-turn, beta-sheet and alpha-helix in your seq + +**Parameters:** +- **seq**: *str* + +input protein seq in 1-letter encoding and upper case + +**Returns:** +- **result**: *list* + +percentage of amino acids belonging to three types of secondary structure for seq + +**Example** +```python +seq = 'AAA' +analyze_secondary_structure(seq) # [0.0, 0.0, 100.0] +``` + +## Contact + +*This is the repo for the 5th homework of the BI Python 2023 course* + +Author: +- *Grishenko Irina* \ No newline at end of file From 247f2abc28736884d7d2a18ff06e5fd24ab9f07f Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Thu, 12 Oct 2023 00:13:10 +0300 Subject: [PATCH 11/30] Add def for read FASTQ-seq and def for creating file with filtered FASTQ-seq --- modules_for_BSAT/fastq_analysis.py | 70 +++++++++++++++++++++++++----- 1 file changed, 60 insertions(+), 10 deletions(-) diff --git a/modules_for_BSAT/fastq_analysis.py b/modules_for_BSAT/fastq_analysis.py index 0dba952..2c792b9 100644 --- a/modules_for_BSAT/fastq_analysis.py +++ b/modules_for_BSAT/fastq_analysis.py @@ -1,15 +1,34 @@ -def is_nucleotide(seq: str) -> bool: +from typing import Union +import os + +def read_fastq(path: str) -> Union[dict, str]: """ - Function return 'True' if sequence is DNA or RNA. - - :param seq: some sequence - :type seq: str - :rtype: bool - :return: 'True' if sequence is DNA or RNA + The function takes as input a file with DNA sequences in FASTQ format and creates a dictionary. The key + is the sequence name, the value is the nucleotide sequence and the quality of reading + :param path: path to the file with sequences in FASTQ format + :type path: str + :rtype: Union [dict, Str] + :return: Dictionary with FASTQ-sequences and path to the file with sequences in FASTQ format. The last one is + necessary for the tool to work correctly. """ - seq_alphabet = {'A', 'T', 'U', 'G', 'C'} - if set(seq.upper()) <= seq_alphabet: - return True + personal_path = path + fastq_dict = {} + name = '' + seq_and_quality = [] + with open (path) as seq_fastq: + for line in seq_fastq: + if line.startswith('@SRX079804'): + name = line.split(' ')[0] + elif set(line.strip()) <= SEQ_ALPHABET: + seq_and_quality += [line.strip()] + elif line.startswith('+'): + continue + else: + seq_and_quality += [line.strip()] + fastq_dict[name] = seq_and_quality + name = '' + seq_and_quality = [] + return fastq_dict, personal_path def analyse_gc(seq: str) -> float: @@ -57,3 +76,34 @@ def analyse_quality(quality: str) -> float: q_score += ord(char) - 33 q_score_sum = q_score/len(quality) return round(q_score_sum,2) + + +def write_fastq(seqs, personal_path, new_file_name): + """ + The function writes filtered FASTQ reads into new file and save it into folder "fastq_filtrator_resuls". + :param seqs: path to the file with sequences in FASTQ format + :type seqs: dict + :param personal_path: path to the file with original sequences in FASTQ format + :type personal_path: str + :param new_file_name: name for file with filtered sequences + :type new_file_name: str + :rtype: None + :return: None + """ + os.mkdir("fastq_filtrator_resuls") + with open (personal_path) as seq_fastq: + res = [] + count = 0 + for line in seq_fastq: + if line.split(' ')[0] in fastq_dict: + count += 1 + res += [line] + continue + count += 1 + res += [line] + if count == 4: + count = 0 + with open (os.path.join('fastq_filtrator_resuls', new_file_name), mode='w') as new_seq_fastq: + for el in res: + new_seq_fastq.write(el) + From f24c0924638633380ba9b2fa2bdefa2035adb31a Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Thu, 12 Oct 2023 00:16:15 +0300 Subject: [PATCH 12/30] Add into def analyze_fastq reading data from a file and writing return to new file --- Bio_Seq_Analysis_Tool.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/Bio_Seq_Analysis_Tool.py b/Bio_Seq_Analysis_Tool.py index 285c1b7..1f1e900 100644 --- a/Bio_Seq_Analysis_Tool.py +++ b/Bio_Seq_Analysis_Tool.py @@ -41,10 +41,10 @@ def dna_rna_analysis(*args: str, operation: str) -> Union[List[float], List[str] return analysis -def analyse_fastq(seqs: dict, +def analyse_fastq(input_path: str, gc_bounds: Union[int, float, Tuple [int], Tuple [float]] = (0, 100), length_bounds: Union[int, Tuple [int]] = (0, 2**32), - quality_threshold: float = 0.0) -> Dict[str,str]: + quality_threshold: float = 0.0, filtered_file_name: Union[None, str] = None) -> Dict[str,str]: """ This function help analyze a set of reads obtained from next-generation sequencing. @@ -52,9 +52,8 @@ def analyse_fastq(seqs: dict, GC-content, length and reading quality. :param seqs: - A dictionary consisting of fastq sequences. The structure is as follows: Key - string, sequence name. - The value is a tuple of two strings: sequence and quality. The sequence is RNA or DNA. - :type seqs: Dict[str] + Path to the file with FASTQ-sequences in the format. + :type seqs: str :param gc_bounds: Boundary parameters for filtering sequences by GC-content. Save only reads with a GC-content between boundaries @@ -79,6 +78,8 @@ def analyse_fastq(seqs: dict, :raises ValueError: if sequence not RNA or DNA, also if the argument values are outside the allowed ones """ + seqs, path_to_file = fq.read_fastq(input_path) + file_name = path_to_file.split("\\")[-1] if type(gc_bounds) == float or type(gc_bounds) == int: gc_bounds = (0,gc_bounds) if gc_bounds[0] < 0 or gc_bounds[1] > 100: @@ -92,8 +93,6 @@ def analyse_fastq(seqs: dict, analysed_seq = {} error_seq = {} for seq in seqs.items(): - if fq.is_nucleotide(seq[1][0]) != True: - raise TypeError(f'Wrong sequence format') if fq.analyse_gc(seq[1][0]) > gc_bounds[0] and fq.analyse_gc(seq[1][0]) < gc_bounds[1]: if fq.analyse_length(seq[1][0]) > length_bounds[0] and fq.analyse_length(seq[1][0]) < length_bounds[1]: if fq.analyse_quality(seq[1][1]) > quality_threshold: @@ -104,7 +103,12 @@ def analyse_fastq(seqs: dict, error_seq[seq[0]] = (seq[1]) else: error_seq[seq[0]] = (seq[1]) - return analysed_seq, error_seq + if filtered_file_name == None: + new_file_name = file_name + else: + new_file_name = filtered_file_name + fq.write_fastq(analysed_seq, path_to_file, new_file_name) + return error_seq def run_protein_analysis(*args: str) -> Union[List[str], List[float]]: From 8222e3f4fb4b29f8c03bb53eab980fbea3812eae Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Thu, 12 Oct 2023 01:12:03 +0300 Subject: [PATCH 13/30] Fixs bag in def read_fastq and write_fastq --- modules_for_BSAT/fastq_analysis.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/modules_for_BSAT/fastq_analysis.py b/modules_for_BSAT/fastq_analysis.py index 2c792b9..cd1ea85 100644 --- a/modules_for_BSAT/fastq_analysis.py +++ b/modules_for_BSAT/fastq_analysis.py @@ -11,6 +11,7 @@ def read_fastq(path: str) -> Union[dict, str]: :return: Dictionary with FASTQ-sequences and path to the file with sequences in FASTQ format. The last one is necessary for the tool to work correctly. """ + seq_alphabet = {'A', 'T', 'G', 'C'} personal_path = path fastq_dict = {} name = '' @@ -19,7 +20,7 @@ def read_fastq(path: str) -> Union[dict, str]: for line in seq_fastq: if line.startswith('@SRX079804'): name = line.split(' ')[0] - elif set(line.strip()) <= SEQ_ALPHABET: + elif set(line.strip()) <= seq_alphabet: seq_and_quality += [line.strip()] elif line.startswith('+'): continue @@ -78,10 +79,10 @@ def analyse_quality(quality: str) -> float: return round(q_score_sum,2) -def write_fastq(seqs, personal_path, new_file_name): +def write_fastq(fastq_dict, personal_path, new_file_name): """ The function writes filtered FASTQ reads into new file and save it into folder "fastq_filtrator_resuls". - :param seqs: path to the file with sequences in FASTQ format + :param seqs: dictionary with filtered sequences :type seqs: dict :param personal_path: path to the file with original sequences in FASTQ format :type personal_path: str From 8333b6634e5174a9270d5eda94d1018d32ff9237 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Thu, 12 Oct 2023 01:13:04 +0300 Subject: [PATCH 14/30] Fix bags in def analyse_fastq --- Bio_Seq_Analysis_Tool.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Bio_Seq_Analysis_Tool.py b/Bio_Seq_Analysis_Tool.py index 1f1e900..b9db4b5 100644 --- a/Bio_Seq_Analysis_Tool.py +++ b/Bio_Seq_Analysis_Tool.py @@ -79,7 +79,7 @@ def analyse_fastq(input_path: str, :raises ValueError: if sequence not RNA or DNA, also if the argument values are outside the allowed ones """ seqs, path_to_file = fq.read_fastq(input_path) - file_name = path_to_file.split("\\")[-1] + file_name = path_to_file.split("/")[-1] if type(gc_bounds) == float or type(gc_bounds) == int: gc_bounds = (0,gc_bounds) if gc_bounds[0] < 0 or gc_bounds[1] > 100: From c7b0739672479175483bcfaac985639bca5eb515 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 20:47:44 +0300 Subject: [PATCH 15/30] Include boundaries for analysis in fastq_analysis in Bio_seq_Analysis_Tool.py --- Bio_Seq_Analysis_Tool.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Bio_Seq_Analysis_Tool.py b/Bio_Seq_Analysis_Tool.py index b9db4b5..3f7e194 100644 --- a/Bio_Seq_Analysis_Tool.py +++ b/Bio_Seq_Analysis_Tool.py @@ -93,8 +93,8 @@ def analyse_fastq(input_path: str, analysed_seq = {} error_seq = {} for seq in seqs.items(): - if fq.analyse_gc(seq[1][0]) > gc_bounds[0] and fq.analyse_gc(seq[1][0]) < gc_bounds[1]: - if fq.analyse_length(seq[1][0]) > length_bounds[0] and fq.analyse_length(seq[1][0]) < length_bounds[1]: + if fq.analyse_gc(seq[1][0]) >= gc_bounds[0] and fq.analyse_gc(seq[1][0]) <= gc_bounds[1]: + if fq.analyse_length(seq[1][0]) >= length_bounds[0] and fq.analyse_length(seq[1][0]) <= length_bounds[1]: if fq.analyse_quality(seq[1][1]) > quality_threshold: analysed_seq[seq[0]] = (seq[1]) else: From 6c5a702f1157970e3d3b425400f503f7486b5ce5 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 20:55:01 +0300 Subject: [PATCH 16/30] Initial commit for Bio_Files_Processor.py --- Bio_Files_Processor.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) 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..e69de29 From 8665b413a117f1130f05278307ae77d7d5903cb5 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 20:59:36 +0300 Subject: [PATCH 17/30] Add import of standard modules in Bio_Files_Processor.py --- Bio_Files_Processor.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Bio_Files_Processor.py b/Bio_Files_Processor.py index e69de29..d96bef2 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -0,0 +1,3 @@ +import os +from typing import List + From f9af641534ad1de42daae609b940719d69970a52 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 21:04:53 +0300 Subject: [PATCH 18/30] Add function 'convert_multiline_fasta_to_oneline' into Bio_Files_Processor.py --- Bio_Files_Processor.py | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/Bio_Files_Processor.py b/Bio_Files_Processor.py index d96bef2..f330c00 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -1,3 +1,42 @@ import os from typing import List +def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = None) -> str: + """ + Function conver multiline fasta file into fasta with name line marked by '>' and the other line with sequence. + Results save into new folder "Converted_data". If the folder already exists, new data is written to it. + :param input_fasta: Path to fasta-file with your seqs with. + It is necessary to indicate the name along with extensions (.fasta) + :type input_fasta: str + :param output_fasta: Name of new fasta-file with your seqs in one line, default value is None + It is necessary to indicate the name along with extensions (.fasta) + :type output_fasta: str + :rtype: str + :return: Script completion message + """ + if os.path.exists(os.path.join('.', 'Converted_data')) == False: + os.mkdir(os.path.join('.', 'Converted_data')) + if output_fasta == None: + out_name = 'one_line_' + input_fasta + else: + out_path = os.path.abspath(output_fasta) + counter = 0 + with open (input_fasta) as seq_fasta: + for line in seq_fasta: + if counter == 0: + if line.startswith('>'): + with open (os.path.join('.', 'Converted_data', out_name), mode = 'w') as new_seq_fasta: + new_seq_fasta.write(line) + counter += 1 + continue + else: + raise ValueError('Wrong start of FASTA') + if counter != 0: + if line.startswith('>'): + with open (os.path.join('.', 'Converted_data', out_name), mode = 'a') as new_seq_fasta: + new_seq_fasta.write('\n' + line) + else: + with open (os.path.join('.', 'Converted_data', out_name), mode = 'a') as new_seq_fasta: + new_seq_fasta.write(line.strip('\n')) + continue + return 'All sequences processed!' \ No newline at end of file From 16fa99a2f31b56561be7b5acc6a9cf3f8c48920a Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 21:10:07 +0300 Subject: [PATCH 19/30] Add function 'select_genes_from_gbk_to_fasta' into Bio_Files_Processor.py --- Bio_Files_Processor.py | 72 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 71 insertions(+), 1 deletion(-) diff --git a/Bio_Files_Processor.py b/Bio_Files_Processor.py index f330c00..6407fb3 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -39,4 +39,74 @@ def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = Non with open (os.path.join('.', 'Converted_data', out_name), mode = 'a') as new_seq_fasta: new_seq_fasta.write(line.strip('\n')) continue - return 'All sequences processed!' \ No newline at end of file + return 'All sequences processed!' + + +def select_genes_from_gbk_to_fasta(input_gbk: str, genes: List[str], n_before: int = 1, n_after: int = 1, + output_fasta: str = None) -> str: + ''' + Function help to search neighbours of GOI (gene of interest). Function writes neighbours of GOI in new FASTA-file + as: name of gene, protein sequence. Results save into new folder "Analyzed_data". If the folder already exists, + new data is written to it. + + :input_gbk: Path to gbk-file with your seqs with. + !! It is necessary to indicate the name along with extensions (.gbk) !! + :type input_gbk: str + :param genes: Gene of interest names + :type genes: str + :param n_before: number of genes before GOI (>0), default value = 1 + :type n_before: int + :param n_after: number of genes after GOI (>0), default value = 1 + :type n_after: int + :output_fasta: Name of FASTA-file with neighbours of GOI (names and seqs), default value = None + :type output_fasta: str + :rtype: str + :return: Script completion message + + ''' + if os.path.exists(os.path.join('.', 'Analyzed_data')) == False: + os.mkdir(os.path.join('.', 'Analyzed_data')) + genes_for_search = genes + genes_gbk = [] + genes_for_search_in_gbk = [] + neighbour_genes = dict() + if output_fasta == None: + output_fasta = 'output_for_gbk.fasta' + with open (input_gbk) as gbk: + for line in gbk: + if '/gene' in line: + genes_gbk += [line.strip().split('=')[1]] + for el in genes_for_search: + genes_for_search_in_gbk += [gn for gn in genes_gbk if el in gn] + for gene in genes_for_search_in_gbk: + gene_index = genes_gbk.index(gene) + if gene_index >= 0 and gene_index < (len(genes_gbk) - 1): + for i in range(1, n_before + 1): + neighbour_genes[(genes_gbk[gene_index - i])] = 0 + for i in range(1, n_after + 1): + neighbour_genes[(genes_gbk[gene_index + i])] = 0 + else: + for i in range(1, n_before + 1): + neighbour_genes[(genes_gbk[gene_index - i])] = 0 + for i in range(1, n_after): + neighbour_genes[(genes_gbk[0 + i])] = 0 + with open (input_gbk) as gbk: + counter = 0 + gene_name = '' + for line in gbk: + for key in neighbour_genes: + if key in line: + counter = 1 + gene_name = key + continue + if counter != 0 and '/translation' in line: + counter = 0 + neighbour_genes[gene_name] = line.strip().split('=')[1] + continue + with open (os.path.join('.', 'Analyzed_data', output_fasta), mode = 'w') as fasta: + for name, seq in neighbour_genes.items(): + name = '>'+name.replace('\"','') + fasta.write(name + '\n') + fasta.write(seq.replace('\"','') + '\n') + return 'All sequences processed!' + From 52eceb89a456dcf8e7424cbdc8a7b944d24441ec Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 21:13:52 +0300 Subject: [PATCH 20/30] Add corrections to the description of the functions into Bio_Files_Processor.py --- 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 6407fb3..460141c 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -5,10 +5,11 @@ def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = Non """ Function conver multiline fasta file into fasta with name line marked by '>' and the other line with sequence. Results save into new folder "Converted_data". If the folder already exists, new data is written to it. + :param input_fasta: Path to fasta-file with your seqs with. - It is necessary to indicate the name along with extensions (.fasta) + !! It is necessary to indicate the name along with extensions (.fasta) !! :type input_fasta: str - :param output_fasta: Name of new fasta-file with your seqs in one line, default value is None + :param output_fasta: Name of new fasta-file with your seqs in one line, default value = None It is necessary to indicate the name along with extensions (.fasta) :type output_fasta: str :rtype: str From 473504966c0a7efc05acf1a3ec5254b0bd8cfc30 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 21:20:36 +0300 Subject: [PATCH 21/30] Fix output_fasta parametr in 'Convert...' function into Bio_Files_Processor.py --- Bio_Files_Processor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Bio_Files_Processor.py b/Bio_Files_Processor.py index 460141c..ff987ce 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -20,7 +20,7 @@ def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = Non if output_fasta == None: out_name = 'one_line_' + input_fasta else: - out_path = os.path.abspath(output_fasta) + out_name = output_fasta counter = 0 with open (input_fasta) as seq_fasta: for line in seq_fasta: From fea880dd8257ffb4c9d07f48f7ee9be21d0839a1 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 21:34:04 +0300 Subject: [PATCH 22/30] Fix input parametr in 'Convert...' function into Bio_Files_Processor.py --- Bio_Files_Processor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Bio_Files_Processor.py b/Bio_Files_Processor.py index ff987ce..d7a4d0c 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -18,7 +18,7 @@ def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = Non if os.path.exists(os.path.join('.', 'Converted_data')) == False: os.mkdir(os.path.join('.', 'Converted_data')) if output_fasta == None: - out_name = 'one_line_' + input_fasta + out_name = 'one_line_' + input_fasta.strip("/")[-1] else: out_name = output_fasta counter = 0 From 0aea07f8ee8f665ea3e8dff408a0a2d9f0abfbe8 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 21:58:06 +0300 Subject: [PATCH 23/30] Add data format check in select_genes_from_gbk_to_fasta --- Bio_Files_Processor.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Bio_Files_Processor.py b/Bio_Files_Processor.py index d7a4d0c..aef42f1 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -54,7 +54,7 @@ def select_genes_from_gbk_to_fasta(input_gbk: str, genes: List[str], n_before: i !! It is necessary to indicate the name along with extensions (.gbk) !! :type input_gbk: str :param genes: Gene of interest names - :type genes: str + :type genes: List[str] :param n_before: number of genes before GOI (>0), default value = 1 :type n_before: int :param n_after: number of genes after GOI (>0), default value = 1 @@ -65,6 +65,8 @@ def select_genes_from_gbk_to_fasta(input_gbk: str, genes: List[str], n_before: i :return: Script completion message ''' + if input_gbk.find('.gbk') == 0: + raise ValueError(f'Wrong file format in input!') if os.path.exists(os.path.join('.', 'Analyzed_data')) == False: os.mkdir(os.path.join('.', 'Analyzed_data')) genes_for_search = genes From 640a3911132143f2676f9d17cfe2d16ddf707319 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 21:59:44 +0300 Subject: [PATCH 24/30] Add data format check in convert_multiline_fasta_to_oneline --- Bio_Files_Processor.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Bio_Files_Processor.py b/Bio_Files_Processor.py index aef42f1..5b0166f 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -15,6 +15,8 @@ def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = Non :rtype: str :return: Script completion message """ + if input_fasta.find('.fasta') == 0: + raise ValueError(f'Wrong file format in input!') if os.path.exists(os.path.join('.', 'Converted_data')) == False: os.mkdir(os.path.join('.', 'Converted_data')) if output_fasta == None: From 180e726b528d88f99c187b5facaacaa14496df2c Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Tue, 17 Oct 2023 22:26:49 +0300 Subject: [PATCH 25/30] Add information about Bio_Files_Processor.py into README.md --- README.md | 83 +++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 81 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 4c49bb8..f81dc14 100644 --- a/README.md +++ b/README.md @@ -19,13 +19,92 @@ Key packages and programs: ## Usage ```python -# import main function +# import main function Bio_Seq_Analysis_Tool from Bio_Seq_Analysis_Tool import Bio_Seq_Analysis_Tool + +# import main function Bio_Files_Processor +from Bio_Files_Processor import Bio_Files_Processor ``` ## Works with main functions -This section contains description of Bio_Seq_Analysis_Tool main functions. +This section contains description of two main scripts: Bio_Files_Processor and Bio_Seq_Analysis_Tool + +## Bio_Files_Processor + +Functions from this script help you to analyze genomic data from different database saved into standard format of these storages. + +### convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = None) -> str + +Function conver multiline fasta file into fasta with name line marked by '>' and the other line with sequence. Results save in new file into new folder "Converted_data". If the folder already exists, new data is written to it. + +**Parameters:** + +-**input_fasta**: *str* + +Path to fasta-file with your seqs with. !! It is necessary to indicate the name along with extensions (.fasta) !! In function exist check for correct format. + +-**output_fasta**: *str* + +Name of new fasta-file with your seqs in one line, default value = None. It is necessary to indicate the name along with extensions (.fasta) + +**Returns:** + +Script completion message + +**Example** + +```python +convert_multiline_fasta_to_oneline('./data_example/example_multiline_fasta.fasta') # 'All sequences processed!' +``` +Structure of new file with FASTA-seq: +'>_name: +sequence' + +### select_genes_from_gbk_to_fasta(input_gbk: str, genes: List[str], n_before: int = 1, n_after: int = 1, output_fasta: str = None) -> str + +Function help to search neighbours of GOI (gene of interest). Function writes neighbours of GOI in new FASTA-file as: name of gene, protein sequence. Results save into new folder "Analyzed_data". If the folder already exists, new data is written to it. + +**Parameters:** + +-**input_gbk**: *str* + +Path to gbk-file with your seqs with. !! It is necessary to indicate the name along with extensions (.gbk) !! In function exist check for correct format. + +-**genes**: *List[str]* + +Gene of interest names. + +Example of input: +```python +['pndA'] +``` +-**n_before**: *int* + +Number of genes before GOI (>0), default value = 1 + +-**n_after**: *int* + +number of genes after GOI (>0), default value = 1 + +-**output_fasta**: *str* + +Name of new fasta-file with your seqs, default value = None. It is necessary to indicate the name along with extensions (.fasta) + +**Returns:** +Script completion message + +**Example** +```python +select_genes_from_gbk_to_fasta('.\\example_data\\example_gbk.gbk', ['pndA'], 2, 2) # 'All sequences processed!' +``` +Structure of new file with FASTA-seq: +'> gene_name: +protein_sequence' + +## Bio_Seq_Analysis_Tool + +Main functions from this script help you to analyze different types of bio sequences: DNA, RNA, NGS-reads, protein ### dna_rna_analysis(*args: str, operation: str) From db4d2ea556a4d2581e93cd39b28e65ff553abe1b Mon Sep 17 00:00:00 2001 From: grishchenkoira <143939574+grishchenkoira@users.noreply.github.com> Date: Tue, 17 Oct 2023 22:52:10 +0300 Subject: [PATCH 26/30] Delete modules_for_BSAT/.ipynb_checkpoints/protein_analysis-checkpoint.py --- .../.ipynb_checkpoints/protein_analysis-checkpoint.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 modules_for_BSAT/.ipynb_checkpoints/protein_analysis-checkpoint.py diff --git a/modules_for_BSAT/.ipynb_checkpoints/protein_analysis-checkpoint.py b/modules_for_BSAT/.ipynb_checkpoints/protein_analysis-checkpoint.py deleted file mode 100644 index e69de29..0000000 From fe687105d2f7acde7b43caff8d7b20936b29267d Mon Sep 17 00:00:00 2001 From: grishchenkoira <143939574+grishchenkoira@users.noreply.github.com> Date: Tue, 17 Oct 2023 22:52:34 +0300 Subject: [PATCH 27/30] Delete modules_for_BSAT/.ipynb_checkpoints/dna_rna_analysis-checkpoint.py --- .../.ipynb_checkpoints/dna_rna_analysis-checkpoint.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 modules_for_BSAT/.ipynb_checkpoints/dna_rna_analysis-checkpoint.py diff --git a/modules_for_BSAT/.ipynb_checkpoints/dna_rna_analysis-checkpoint.py b/modules_for_BSAT/.ipynb_checkpoints/dna_rna_analysis-checkpoint.py deleted file mode 100644 index e69de29..0000000 From 9e52d572d997e63839a40601d80bcc15b27697ed Mon Sep 17 00:00:00 2001 From: grishchenkoira <143939574+grishchenkoira@users.noreply.github.com> Date: Tue, 17 Oct 2023 22:52:49 +0300 Subject: [PATCH 28/30] Delete modules_for_BSAT/.ipynb_checkpoints/fastq_analysis-checkpoint.py --- modules_for_BSAT/.ipynb_checkpoints/fastq_analysis-checkpoint.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 modules_for_BSAT/.ipynb_checkpoints/fastq_analysis-checkpoint.py diff --git a/modules_for_BSAT/.ipynb_checkpoints/fastq_analysis-checkpoint.py b/modules_for_BSAT/.ipynb_checkpoints/fastq_analysis-checkpoint.py deleted file mode 100644 index e69de29..0000000 From a1d5be10008fa6ddfc7f0d3e3f130edfc4d9b841 Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Wed, 18 Oct 2023 22:52:54 +0300 Subject: [PATCH 29/30] Rewritten code to add complete protein sequence in select_genes_from_gbk_to_fasta --- Bio_Files_Processor.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/Bio_Files_Processor.py b/Bio_Files_Processor.py index 5b0166f..72fcbe9 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -96,18 +96,25 @@ def select_genes_from_gbk_to_fasta(input_gbk: str, genes: List[str], n_before: i for i in range(1, n_after): neighbour_genes[(genes_gbk[0 + i])] = 0 with open (input_gbk) as gbk: - counter = 0 + gene_name_read = 0 + protein_read = 0 gene_name = '' + protein = '' for line in gbk: - for key in neighbour_genes: - if key in line: - counter = 1 - gene_name = key - continue - if counter != 0 and '/translation' in line: - counter = 0 - neighbour_genes[gene_name] = line.strip().split('=')[1] - continue + if '/gene' in line: + gene_name = line.strip().split('=')[1] + if gene_name in neighbour_genes: + gene_name_read = 1 + if protein_read == 1: + protein += line.strip('\n').strip(' ') + if '"' in line: + protein_read = 0 + gene_name_read = 0 + neighbour_genes[gene_name] = protein + protein = '' + if gene_name_read == 1 and '/translation' in line: + protein_read = 1 + protein += line.strip().split('=')[1] with open (os.path.join('.', 'Analyzed_data', output_fasta), mode = 'w') as fasta: for name, seq in neighbour_genes.items(): name = '>'+name.replace('\"','') From dca7a004202832423878150675e89262b3c7297b Mon Sep 17 00:00:00 2001 From: Irina Grishchenko Date: Thu, 19 Oct 2023 01:40:34 +0300 Subject: [PATCH 30/30] Rewrite function to generate a FASTA-file --- Bio_Files_Processor.py | 36 +++++++++++++++++++----------------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/Bio_Files_Processor.py b/Bio_Files_Processor.py index 72fcbe9..3872474 100644 --- a/Bio_Files_Processor.py +++ b/Bio_Files_Processor.py @@ -23,26 +23,28 @@ def convert_multiline_fasta_to_oneline(input_fasta: str, output_fasta: str = Non out_name = 'one_line_' + input_fasta.strip("/")[-1] else: out_name = output_fasta - counter = 0 + gene = 0 + seq = 0 + gene_name = '' + prot_seq = '' + gene_and_seq = dict() with open (input_fasta) as seq_fasta: for line in seq_fasta: - if counter == 0: - if line.startswith('>'): - with open (os.path.join('.', 'Converted_data', out_name), mode = 'w') as new_seq_fasta: - new_seq_fasta.write(line) - counter += 1 - continue + if gene == 1: + if line.startswith('>') == False: + gene_and_seq[gene_name] += line.strip('\n') else: - raise ValueError('Wrong start of FASTA') - if counter != 0: - if line.startswith('>'): - with open (os.path.join('.', 'Converted_data', out_name), mode = 'a') as new_seq_fasta: - new_seq_fasta.write('\n' + line) - else: - with open (os.path.join('.', 'Converted_data', out_name), mode = 'a') as new_seq_fasta: - new_seq_fasta.write(line.strip('\n')) - continue - return 'All sequences processed!' + gene_name = line.strip('\n') + gene_and_seq[gene_name] = '' + if gene == 0 and line.startswith('>'): + gene_name = line.strip('\n') + gene_and_seq[gene_name] = '' + gene = 1 + with open (os.path.join('.', 'Converted_data', out_name), mode = 'w') as new_seq_fasta: + for key, item in gene_and_seq.items(): + new_seq_fasta.write(key + '\n') + new_seq_fasta.write(item + '\n') + return 'All sequences processed!' def select_genes_from_gbk_to_fasta(input_gbk: str, genes: List[str], n_before: int = 1, n_after: int = 1,