diff --git a/README.md b/README.md index 39f073e..bbf21e9 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,176 @@ # 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` and 'bio_files_processor.py' are an open-source program that facilitate working with bioinformatics data. + +## Usage and options + +### Bioinf_tools.py +`bioinf_tools.py` program contains 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 as usage example.\ +To start with the program run the following command: + +`run_fastq_filter(gc_bounds=50)` + +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. + +### 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 + + +## 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", "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' + +### 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')} + +### 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") +``` + diff --git a/bio_files_processor.py b/bio_files_processor.py new file mode 100644 index 0000000..5ba06dd --- /dev/null +++ b/bio_files_processor.py @@ -0,0 +1,136 @@ +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: str or 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 (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 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) + if isinstance(genes, str): + genes = [genes] + 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") diff --git a/bioinf_tools.py b/bioinf_tools.py new file mode 100644 index 0000000..670cd05 --- /dev/null +++ b/bioinf_tools.py @@ -0,0 +1,152 @@ +from src import dna_rna_tools +from src import fastq_filter +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( + 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: 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: + - 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 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 fastq file + - save_to_dir (str): absolute or realtive path to directory to save to + - output_filename (str): output name of the filtered fastq file + For examples see 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, + output_filename, + ) = fastq_filter.parse_and_check_user_input( + sequences_path, + gc_bounds, + length_bounds, + quality_threshold, + verbose, + output_filename, + ) + 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/dictionaries.py b/src/dictionaries.py new file mode 100644 index 0000000..7efc1e3 --- /dev/null +++ b/src/dictionaries.py @@ -0,0 +1,110 @@ +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" +} +TRANSLATION_RULE = { + "F": "UUU", + "L": "CUG", + "I": "AUU", + "M": "AUG", + "V": "GUG", + "P": "CCG", + "T": "ACC", + "A": "GCG", + "Y": "UAU", + "H": "CAU", + "Q": "CAG", + "N": "AAC", + "K": "AAA", + "D": "GAU", + "E": "GAA", + "C": "UGC", + "W": "UGG", + "R": "CGU", + "S": "AGC", + "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, +} +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/src/dna_rna_tools.py b/src/dna_rna_tools.py new file mode 100644 index 0000000..6a6f3ae --- /dev/null +++ b/src/dna_rna_tools.py @@ -0,0 +1,120 @@ +if __name__ == "__main__": + import dictionaries +else: + from src import dictionaries + +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. + + 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: + transcript = sequence.replace("T", "U").replace("t", "u") + return transcript + + +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: 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): + dictionaries.NUCL_ACIDS_COMPLEMENT_RULE["A"] = "U" + complement_strand = "" + for letter in sequence: + 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: str) -> str: + """ + Convert nucleic acid sequence into its reverse-complement counterpart.\n + Letter case sensitive.\n + + Arguments: + - sequence (str): sequence to convert + + Return: + - reverse_complement_strand (str): reverse-comlepement sequence + """ + reverse_complement_strand = complement(sequence) + reverse_complement_strand = reverse(reverse_complement_strand) + return reverse_complement_strand + + +DNA_RNA_PROCEDURES_FUNCTIONS = { + "transcribe": transcribe, + "reverse": reverse, + "complement": complement, + "reverse_complement": reverse_complement, +} diff --git a/src/fastq_filter.py b/src/fastq_filter.py new file mode 100644 index 0000000..28becd3 --- /dev/null +++ b/src/fastq_filter.py @@ -0,0 +1,186 @@ +if __name__ == "__main__": + import dictionaries +else: + from src import dictionaries +import os + + +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: bool, + output_filename: str, +): + """ + 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 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 fastq 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: + seqs[line.strip()] = [] + 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.strip()) + 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") + if verbose != True and verbose != False: + raise ValueError("Invalid *verbose* argument given") + return seqs, gc_bounds, length_bounds, quality_threshold, verbose, output_filename + + +def fastq_filter( + 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: + """ + Filter out bad reads.\n + + Arguments: + - 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 + - verbose (bool): add detailed statistics for each read + + 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, verbose) + 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 + + +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]), verbose +) -> 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 + - verbose (bool): add detailed statistics for each read + + 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 + 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]), verbose +) -> 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 + - verbose (bool): add detailed statistics for each read + + Return: + - condition (bool): True - length is within bounds, False - read is to be filtered + """ + if isinstance(length_bounds, int): + length_bounds = (1, length_bounds) + 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, verbose) -> 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 + - verbose (bool): add detailed statistics for each read + + 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) + 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 fastq format. + + Arguments: + - 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) + 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 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])) diff --git a/src/protein_tools.py b/src/protein_tools.py new file mode 100644 index 0000000..7a248cd --- /dev/null +++ b/src/protein_tools.py @@ -0,0 +1,291 @@ +if __name__ == "__main__": + import dictionaries +else: + from src import dictionaries + +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 + + 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: + 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: + 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: (tuple[str] or list[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: (tuple[str] or list[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: (tuple[str] or list[str]), 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 = 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) + if nucl_acids == "RNA": + nucl_acid_seqs["RNA"].append(rna_seq) + elif nucl_acids == "DNA": + dna_seq = rna_seq.replace("U", "T").replace("u", "t") + nucl_acid_seqs["DNA"].append(dna_seq) + 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 + + +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, + "three_one_letter_code": three_one_letter_code, + "define_molecular_weight": define_molecular_weight, +} + + +def check_and_parse_user_input( + 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 + + Return: + - string: procedure name + - dictionary: a collection of procedure arguments and their values + """ + 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 PROTEINS_PROCEDURES_TO_FUNCTIONS.keys(): + raise ValueError("Wrong procedure") + allowed_inputs = set(dictionaries.AMINO_ACIDS.keys()) + allowed_inputs = allowed_inputs.union( + set(k.lower() for k in dictionaries.AMINO_ACIDS.keys()) + ) + 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()) + 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: + if not set(sequence).issubset(allowed_inputs_seq): + 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