diff --git a/README.md b/README.md index f3e345c..16c6301 100644 --- a/README.md +++ b/README.md @@ -1 +1,261 @@ -# Junior_Bioinformatic_tool \ No newline at end of file +# My junior biochemical tool +This package consists of 2 mini-tools +`tools_for_bioinformatics` designed to work with nucleotide and amino acid sequences and FASTQ files +`bio_files_processor` provides opportunity to manipulate with FASTA files + + +# My junior biochemical tool +This package for junior bioinformaticians is designed to work with nucleotide and amino acid sequences and FASTQ files, and consists of three parts. Each module `amino_acid_tools`, `dna_rna_tools` and `fastq_filtration` presents its own possibilities for sequence processing. + +# `amino_acids_tools` +This tool is designed to work with amino acid sequences consisting of _22 proteinogenic amino acid_ residues (including pyrrolizine and selenocysteine) recorded in a standard one-letter format. It is not intended to process sequences with post-translational and other amino acid modifications + +## Usage +You call the `amino_acid_tools` function, which takes as input an arbitrary number of arguments with amino-acid sequences (str), as well as the name of the procedure to be executed (it is always the last argument, str). After that the command performs the specified action on all the given sequences. If one sequence is submitted, a string with the result is returned. If several sequences are submitted, a list of strings is returned. +Input sequences can contain both uppercase and lowercase letters, but the last argument with the function name must correspond to the listed functions. + + +## Options +The following options for aminoacid sequence processing are available at the moment: + +- **molecular_weight**: calculate the molecular weight of the amino acid chain in Da, according to the average amino acid residues molecular masses rounded to 1 or 2 decimal places. +- **three_letter_code**: converts standard single letter translations to three letter translations +- **show_length**: count the overall number of amino acids in the given +- **sequence folding**: count the number of amino acids characteristic separately for alpha helixes and beta sheets,and give out what will be the structure of the protein more. This function has been tested on proteins such as 2M3X, 6DT4 (PDB ID) and MHC, CRP. The obtained results corresponded to reality. +- **seq_charge**: evaluates the overall charge of the aminoacid chain in neutral aqueous solution (pH = 7), according to the pKa of amino acid side chains, lysine, pyrrolizine and arginine contribute +1, while asparagine and glutamic amino acids contribute -1. The total charge of a protein is evaluated as positive, negative, or neutral as the sum of these contributions + +### Remarks +If sequense contains symbols differ from IUPAC 1-letter code for 22 proteinogenic amino acids, +the result *for this sequence* will be `'unexpected symbols in sequence'` regardless of function +If action is not in the function the message `'unexpected action'` will occur + +## Examples +Below is an example of processing an amino acid sequence for different input data +```Python +amino_acid_tools('DNA', 'molecular_weight') +``` +Output: 300.28 + +```Python +amino_acid_tools('DNA', 'russco', 'LOLkEk', 'azaz', 'three_letter_code') +``` +Output: ['AspAsnAla', 'ArgSecSerSerCysPyl', 'LeuPylLeuLysGluLys', 'unexpected symbols in sequence'] + +```Python +amino_acid_tools('DNA', 'russco', 'LOLkEk', 'azaz', 'letter_code') +``` +ValueError: Unexpected action + +### Using the function for molecular weight calculation + +```Python +amino_acid_tools('EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'molecular_weight') +``` + +Input: 'EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'molecular_weight' +Output: '[1228.66, 1447.8400000000001, 1224.6399999999999]' + +### Using the function to convert one-letter translations to three-letter translations + +```Python +amino_acid_tools('EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'three_letter_code') +``` + +Input: 'EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'three_letter_code' +Output: '['GluGlyValIleMetSerGluLeuLysLeuLys', 'ProLeuProLysValGluLeuProProAspPheValAsp', 'AspValIleGlyIleSerIleLeuGlyLysGluVal']' + +### Using the function to counts the number of amino acids in the given sequence + +```Python +amino_acid_tools('EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'show_length') +``` + +Input: 'EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'show_length' +Output: '[11, 13, 12]' + +### Using the function to determine the predominant secondary structure + +```Python +amino_acid_tools('EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'folding') +``` +Input: 'EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'folding' +Output: '['alfa_helix', 'equally', 'equally']' + +### Using the function to estimate relative charge + +```Python +amino_acid_tools('EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'seq_charge') +``` + +Input: 'EGVIMSELKLK', 'PLPKvelPPDFVD', 'DVIGISILGKEV', 'seq_charge' +Output: '['neutral', 'negative', 'negative']' + +# `DNA_RNA_tools` + +This program is designed to work with nucleic acid sequences. We expect RNA and DNA chains in the standard NNNNNNNN form, without phosphate groups, in uppercase or lowercase letters + +### Usage +You call the `dna_rna_tools` function, which takes as input anumber of arguments with amino-acid sequences (str), as well as the name of the procedure to be executed (it is always the last argument, str). + +The nucleotide sequence can consist of both uppercase and lowercase letters. +Input example: +```Python +dna_rna_tools('GCGGT','auccuc','GCTatGc','complement') +``` + +### Functions: +Folowing options are avaliable now: + + - **complement:** makes complement sequence to seq + - **reverse**: reverses sequence, (from 5'-3' to 3' -5' there or back) + - **reverse_complement**: makes complement sequence and reverse it + - **transcribe**: make transcript of DNA sequence. in the case of RNA + sequence, no changes will be returned + +*Returns:* +If one sequence is supplied, a string with the result is returned. +If several are submitted, a list of strings is returned. + +#### Remarks +If seq contains both U(u) and T(t) the result *for this sequence* will be `"U and T in one seq"` regardless of function + +If seq contains symbols differ from standardized oligonucleotide notation, the result *for this sequence* will be `'unexpected symbols in sequence'` regardless of function +If action is not in the function the message `'unexpected action'` will occur + +### Examples: +```Python +dna_rna_tools('ATTC', 'CGcGc', 'AZA', 'atuc', 'Au', 'transcribe') +``` +Output: ['AUUC', 'CGcGc', 'unexpected symbols in sequence', 'T and U in one seq', 'Au'] + + +```Python +dna_rna_tools('ATTC', 'CGcGc', 'AZA', 'atuc', 'Au', 'revese') +``` +ValueError: Unexpected action + +```Python +dna_rna_tools('UuCG','complement') +``` +Output: 'AaGC' + +# `FASTQ_filtration tool` + +This function provides you the opportunity to filter the FASTQ list to select sequences according to requirements on three parameters: length, GC composition, and quality of the reed + +### Usage +It is required to input the list of sequences in dictionary format, as well as values for filtering parameters: interval by GC-composition, interval by sequence length, threshold value of average quality of reed. You also need input file and output file name (optional) + **input_fastq:** FASTQ file +**output_fastq** name of output file, you can find it in fastq_filtrator_resuls directory, if name is not defined, it will be the same of the input file in `fastq_filtration_result` directory + +```Python +fastq_filtration('example_fastq.fastq', , gc_bounds=(20, 40), length_bounds=(0, 2 *8), quality_treshold=10, 'outfastq_gc' ) +``` + +#### Filtering parameters + + - **gc_bounds:** interval for the of acceptable GC content, in %, *Default + value = (0,100)* + - **length_bounds**: interval for the of acceptable sequense length in + number of nucleotide, *Default value = (2,2**32)* + - **quality_treshold**: threshold value for average quality per nucleotide + (phred33 scale), *Default value = 0* + + +### Result: +New FASTQ file consists of selected sequences after 3-step filtration +#### Remarks +After running without specifying the filtering parameters, all sequences will be be selected as appropriate +You also can specify only the upper limit for GC-content and length filtering + +### Examples +```Python +fastq_filtration('example_fastq.fastq', gc_bounds= (80, 100)) +``` +@SRX079804:1:SRR292678:1:1101:1120907:1120907 1:N:0:1 BH:ok +CATGGTGGCG ++SRX079804:1:SRR292678:1:1101:1120907:1120907 1:N:0:1 BH:ok +HHHHHFHHGG +@SRX079804:1:SRR292678:1:1101:1175112:1175112 1:N:0:1 BH:failed +AGGCC ++SRX079804:1:SRR292678:1:1101:1175112:1175112 1:N:0:1 BH:failed +EC8EE +@SRX079804:1:SRR292678:1:1101:1269735:1269735 1:N:0:1 BH:failed +C ++SRX079804:1:SRR292678:1:1101:1269735:1269735 1:N:0:1 BH:failed +G +```Python +fastq_filtration('example_fastq.fastq', length_bounds= (8, 12),output_fastq = 'outfastq_length') +``` +@SRX079804:1:SRR292678:1:1101:1119783:1119783 1:N:0:1 BH:failed +TGTTGTGTCAAG ++SRX079804:1:SRR292678:1:1101:1119783:1119783 1:N:0:1 BH:failed +C@ABABEEEEAC +@SRX079804:1:SRR292678:1:1101:1120614:1120614 1:N:0:1 BH:failed +TGCTTTGCTTT ++SRX079804:1:SRR292678:1:1101:1120614:1120614 1:N:0:1 BH:failed +DFDD8FE@FFE +@SRX079804:1:SRR292678:1:1101:1120907:1120907 1:N:0:1 BH:ok +CATGGTGGCG ++SRX079804:1:SRR292678:1:1101:1120907:1120907 1:N:0:1 BH:ok +HHHHHFHHGG +@SRX079804:1:SRR292678:1:1101:1130921:1130921 1:N:0:1 BH:failed +AAGGGTCGA ++SRX079804:1:SRR292678:1:1101:1130921:1130921 1:N:0:1 BH:failed +D@CDDBEEB +@SRX079804:1:SRR292678:1:1101:1156698:1156698 1:N:0:1 BH:failed +GAAAGAAC ++SRX079804:1:SRR292678:1:1101:1156698:1156698 1:N:0:1 BH:failed +E=FFFFFF +```Python +fastq_filtration('example_fastq.fastq') +``` +input file (no filtratation with default parameters) +```Python +fastq_filtration('example_fastq.fastq', gc_bounds=(100, 100), length_bounds=(0, 0 ), quality_treshold=0, output_fastq = 'outfastq_none') +``` +empty file + + +# `bio_files_processor` + +works with FASTA files, provides 2 functions + +## `convert_multiline_fasta_to_oneline +This function covert multiple string FASTA-file to one-long-string format + + - **input_fasta**: path to the FASTA file + - **output_fasta**: name for output FASTA-file, if not defined it will be saved to the directory `Result` as `'Result'/input_fasta` + +Example input: +```Python +convert_multiline_fasta_to_oneline('example_multiline_fasta.fasta', 'output2'): +``` + +## `change_fasta_start_pos` + + +This function moves the position to a new starting nucleotide and needs 3 parameters + + - **input_fasta:** data for processing (1-string single fasta) + - **n:** position of new start nucleotide (started with 0) + - **output_fasta:** name for output fasta file + + + +Example input: +```Python +change_fasta_start_pos('example_fasta.fasta', 3, 'output2'): +``` + +Examples: +ACTGGA initial sequence +TGGAAC n =2 +AACTGG n =-1 +AACTGG n = 5 + + + +##### Contacts +Maria Lukina +maria.v.luk@gmail.com \ No newline at end of file diff --git a/bio_files_processor.py b/bio_files_processor.py new file mode 100644 index 0000000..b65fe39 --- /dev/null +++ b/bio_files_processor.py @@ -0,0 +1,68 @@ +import os + + +def convert_multiline_fasta_to_oneline(input_fasta, output_fasta = ''): + """ + This function covert multiple string FASTA-file to one-long-string format + :param input_fasta: the pass to the FASTA file + :param output_fasta: name for output FASTA-file, if not defined it will be saved to the directory + 'Result'/input_fasta + Example input: convert_multiline_fasta_to_oneline('example_multiline_fasta.fasta', 'output2'): + """ + if not os.path.isdir('multiple_to_online_results'): + os.mkdir('multiple_to_online_results') + if output_fasta == '': + output_fasta = os.path.join('multiple_to_online_results', os.path.basename(input_fasta)) + else: + output_fasta = os.path.join('multiple_to_online_results', output_fasta + ".fasta") + with open(input_fasta) as input_file, open(output_fasta, 'w') as output_file: + current = [] + output_file.write(input_file.readline()) + while True: + line = input_file.readline() + current.append(line.strip()) + if line.startswith('>'): + output_file.write(''.join(current) + '\n') + output_file.write(line) + current = [] + break + + for line in input_file: + if line.startswith('>'): + output_file.write(''.join(current) + '\n') + output_file.write(line) + current = [] + else: + current.append(line.strip()) + output_file.write(''.join(current) + '\n') + + +def change_fasta_start_pos(input_fasta, n: int, output_fasta=''): + """ + This function moves the position to a new starting nucleotide. + :param input_fasta: data for processing (1-string single fasta) + :param n: position of new start nucleotide (started with 0) + :param output_fasta: name for output fasta file + """ + if not os.path.isdir('shifted_fasta_results'): + os.mkdir('shifted_fasta_results') + if output_fasta == '': + output_fasta = os.path.join('shifted_fasta_results', os.path.basename(input_fasta)) + else: + output_fasta = os.path.join('shifted_fasta_results', output_fasta + ".fasta") + with open(input_fasta) as input_file: + name = input_file.readline().strip() + string = input_file.readline().strip() + s1 = '' + s2 = '' + if n <= 0: + n = len(string) + n + for i in range(0, len(string)): + if i < n: + s2 = s2 + string[i] + else: + s1 = s1 + string[i] + new_fasta = s1 +s2 + with open(output_fasta, 'w') as output_file: + output_file.write(name +'\n') + output_file.write(new_fasta[0] + '\n') diff --git a/requirments.txt b/requirments.txt new file mode 100644 index 0000000..04e9ef2 --- /dev/null +++ b/requirments.txt @@ -0,0 +1,2 @@ +Python 3.12.0 +biopython==1.83 \ No newline at end of file diff --git a/tool_for_bioinformatics.py b/tool_for_bioinformatics.py new file mode 100644 index 0000000..f2acf42 --- /dev/null +++ b/tool_for_bioinformatics.py @@ -0,0 +1,190 @@ +import os +from abc import ABC, abstractmethod +from Bio import SeqIO +from Bio import SeqUtils + + +class BiologicalSequence(ABC, str): + + def __init__(self, seq): + self.seq = seq + + #да, Никита уже сказал, что так лучше не делать, а оставить сдесь только абстрактные методы, + # но я зацепилась за идею не повторять проверки, раз суть для всех одна, только алфавиты разные, + # а переделать уже не успеваю + def check_alphabet(self): + if set(self.seq) <= self.alphabet: + raise UnexpectedSymbolInSeqError + + # не поняла, не доделала + @abstractmethod + def to_print_seq(self): + return self.seq + pass + + +class UnexpectedSymbolInSeqError(ValueError): + pass + + +class AminoAcidSequence(BiologicalSequence): + alphabet = {'M', 'O', 'v', 'D', 'f', 'N', 'c', 'A', 'R', 'W', 'I', 'm', 'L', 's', 'H', 'q', 'w', 'V', 'n', 'i', + 'g', 'F', 'S', 'e', 'l', 'U', 'P', 'Q', 'K', 'Y', 'u', 'y', 'd', 'h', 'k', 'r', 't', 'G', 'o', 'E', + 'p', 'T', 'C', 'a'} + masses = {'A': 71.08, 'R': 156.2, 'N': 114.1, 'D': 115.1, 'C': 103.1, 'E': 129.1, 'Q': 128.1, 'G': 57.05, 'H': 137.1, + 'I': 113.2, 'L': 113.2, 'K': 128.2, 'M': 131.2, 'F': 147.2, 'P': 97.12, 'S': 87.08, 'T': 101.1, + 'W': 186.2, 'Y': 163.2, 'V': 99.13, 'U': 168.05, 'O': 255.3, 'a': 71.08, 'r': 156.2, 'n': 114.1, 'd': 115.1, + 'c': 103.1, 'e': 129.1, 'q': 128.1, 'g': 57.05, 'h': 137.1, 'i': 113.2, 'l': 113.2, 'k': 128.2, 'm': 131.2, + 'f': 147.2, 'p': 97.12, 's': 87.08, 't': 101.1, 'w': 186.2, 'y': 163.2, 'v': 99.13, 'u': 168.05, 'o': 255.3} + + def __init__(self, seq): + super().__init__(seq) + self.seq = seq + self.check_alphabet() + def molecular_weight(self) -> float: + """ + Function calculates molecular weight of the amino acid chain + Input Parameters: + each letter refers to one-letter coded proteinogenic amino acids + Returns: + (float) Molecular weight of tge given amino acid chain in Da + """ + m = 0 + for acid in str(self.seq): + m += self.masses[acid] + return m + +class NucleicAcidSequence(BiologicalSequence): + def __init__(self, seq): + super().__init__(seq) + self.check_alphabet() + + def gc_content(self) -> float: + """ + Function counts GC-content in sequence, and returns result in % + """ + n = 0 + for nucl in self.seq: + if nucl == 'c' or nucl == 'g' or nucl == 'C' or nucl == 'G': + n += 1 + return 100 * n / len(self.seq) + + def complement(self): + """ + Function return complement sequence + """ + complementary_dna = self.seq.maketrans(self.dict_comp) + res = self.seq.translate(complementary_dna) + return NucleicAcidSequence(res) + +class DNASequence(BiologicalSequence): + alphabet = {'A', 'g', 't', 'G', 'T', 'a', 'c', 'C'} + dict_trans = {'A': 'A', 'C': 'C', 'T': 'U', 'G': 'G', 'a': 'a', 'c': 'c', 't': 'u', 'g': 'g'} + dict_comp = {'A': 'T', 'C': 'G', 'T': 'A', 'G': 'C', 'a': 't', 'c': 'g', 't': 'a', 'g': 'c'} + def __init__(self, seq): + super().__init__(seq) + + + def transcribe(self): + """ + Function return transcript of DNA sequence + """ + transcribe = self.seq.maketrans(self.dict_trans) + res = self.seq.translate(transcribe) + + return RNASequence(res) + + +class RNASequence(BiologicalSequence): + alphabet = {'U', 'A', 'g', 'G', 'a', 'c', 'C', 'u'} + dict_comp = {'A': 'U', 'C': 'G', 'U': 'A', 'G': 'C', 'a': 'u', 'c': 'g', 'u': 'a', 'g': 'c'} + def __init__(self, seq): + super().__init__(seq) + + + #фильтратор +def filter_gc(records, gc_bounds_both_side=(0, 100)) -> list: + """ + This function selects sequences with the GC content of your interest + :parameters: + records: records from fastq parced by SeqIO + gc_bound: interval for the of acceptable GC content, in % + :return:(dict) new dictionary consists of selected sequences + """ + new_records = [] + for record in records: + if (gc_bounds_both_side[1]/100 >= SeqUtils.gc_fraction(record.seq) >= gc_bounds_both_side[0]/100): + new_records.append(record) + + return new_records + + +def filter_length(records, length_bounds_both_side=(0, 2 ** 32)) -> list: + """ + This function selects sequences with the length of your interest + :parameters: + records: records from fastq parced by SeqIO + length_bound: interval for the of acceptable sequense length in number of nucleotide + :return:(dict) new dictionary consists of selected sequences + """ + new_records = [] + for record in records: + if (length_bounds_both_side[1] >= len(record.seq) >= length_bounds_both_side[0]): + new_records.append(record) + print(new_records) + return new_records + + + +def filter_quality(records, quality_threshold=0) -> list: + """ + This function selects FASTQ sequences with appropriate average nucleotide read quality + parameters: + seqs: dictionary of FASTQ sequences {name: (sequence, quality)} + quality_treshold: threshold value for average quality per nucleotide (phred33 scale) + :return:(dict) recordes for selected sequences + """ + new_records = [] + for record in records: + if (sum(record.letter_annotations["phred_quality"])/len(record.seq) >= quality_threshold): + new_records.append(record) + print(new_records) + return new_records + + +def fastq_filtration(input_fastq, gc_bounds=(0, 100), length_bounds=(0, 2 ** 32), quality_treshold=0, output_fastq=''): + """ + This function provides you the opportunity to filter the FASTQ file to select sequences + that fit requirements on 5 parameters: input and output(optional) files, length, GC composition, + and quality of the reed + :parameters + input_fastq: path to fastq file + gc_bounds: (tuple) interval for the of acceptable GC content, in % + length_bounds: (tuple) interval for the of acceptable sequense length in number of nucleotide + quality_treshold: (float) threshold value for average quality per nucleotide (phred33 scale) + output_fastq = name of output file, ./fastq_filtrator_resuls/output_fastq, if it is not defined, + it will be the same of the input file + + """ + if not os.path.isdir('fastq_filtrator_resuls'): + os.mkdir('fastq_filtrator_resuls') + if output_fastq == '': + output_fastq = os.path.join('fastq_filtrator_resuls', os.path.basename(input_fastq)) + else: + output_fastq = os.path.join('fastq_filtrator_resuls', output_fastq + ".fasta") + if type(gc_bounds) == float or type(gc_bounds) == int: + gc_bounds_both_side = (0, gc_bounds) + else: + gc_bounds_both_side = gc_bounds + if type(length_bounds) == int: + length_bounds_both_side = (0, length_bounds) + else: + length_bounds_both_side = length_bounds + records = list(SeqIO.parse(input_fastq, "fastq")) + filtered_records_l = filter_length(records, length_bounds_both_side) + filtered_records_gc = filter_gc(filtered_records_l, gc_bounds_both_side) + filtered_records_q = filter_quality(filtered_records_gc, quality_treshold) + SeqIO.write((record for record in filtered_records_q), output_fastq, "fastq") + return + +#fastq_filtration('example_fastq.fastq', length_bounds= (10, 30), gc_bounds=(20, 30), output_fastq='')