diff --git a/README.md b/README.md index 16c6301..24140ea 100644 --- a/README.md +++ b/README.md @@ -1,261 +1,37 @@ -# 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 first bioinformatics tools +This repository contains homeworks assignments from the Python course within the "Bioinformatics for Biologists" annual program at the Institute of Bioinformatics (2023-2024). It contains python realizations of a number of training tasks -# 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). +### Contents: -The nucleotide sequence can consist of both uppercase and lowercase letters. -Input example: -```Python -dna_rna_tools('GCGGT','auccuc','GCTatGc','complement') -``` +* ##### `main_script`: + * `run_genscan` function and `GenscanOutput` class: Retrieve data from the GENSCAN website (http://argonaute.mit.edu/GENSCAN.html). + * `telegram_logger` decorator: Sends a Telegram bot message with the function's completion time and any error messages. + * `fastq_filtration` function: Filters reads from FASTQ files based on length, GC content, and read quality, returning a file containing filtered sequences. + * `BiologicalSequence`, `NucleicAcidSequence`, `DNASequence`, `RNASequence`, `AminoAcidSequence` classes: Enable various operations on biological sequences. +* ##### `bio_files_processor`: + * `OpenFasta` class: Reads and processes multi-line FASTA files record by record, providing results in a user-friendly format using the `FastaRecord` dataclass. + * `convert_multiple_fasta_to_online` function: Converts a multi-line FASTA file to a single-line format. + * `change_fasta_start_pos` function: Shifts the starting position to a specified nucleotide. +* ##### `test_modules.py`: + Contains a set of tests for the `change_fasta_start_pos` function and the `DNASequence`, `NucleicAcidSequence` classes. +* ##### `custom_random_forest.py`: + An attempt to accelerate the creation of an ensemble of models using parallel programming techniques. +_Note: The `custom_random_forest` file is currently under development and will be added to the repository soon_. Now you can find draft notebook with realization and working examples -### Functions: -Folowing options are avaliable now: +#### Showcases - - **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 +The `Showcases.ipynb` notebook provides examples of how to use the following: -*Returns:* -If one sequence is supplied, a string with the result is returned. -If several are submitted, a list of strings is returned. +* `telegram_logger` +* `RunGenscan` +* `BiologicalSequences` classes -#### 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 +_Feedback and Contact:_ -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` +I welcome any questions, comments, or suggestions. Please feel free to reach out to me: + Maria Lukina -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 - - + Email: _maria.v.luk@gmail.com_ Telegram: _earthling_maria_ -##### Contacts -Maria Lukina -maria.v.luk@gmail.com \ No newline at end of file diff --git a/Showcases.ipynb b/Showcases.ipynb new file mode 100644 index 0000000..4b748bb --- /dev/null +++ b/Showcases.ipynb @@ -0,0 +1,271 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### В этом ноутбуке представлены примеры выполнения функций и работы классов" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 1. `run_genscan`" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from main_script import run_genscan" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "GenscanOutput🦠 :\n", + "Sequence name is: ---\n", + "📬Request status code is <200>\n", + "📈exon_list= {'Intr_1': (383, 474), 'Intr_2': (1215, 1290), 'Intr_3': (1927, 2101), 'Term_4': (2887, 3015), 'PlyA_5': (3036, 3041)}\n", + "📉intron_list= {'1': (1, 382), '2': (475, 1214), '3': (1291, 1926), '4': (2102, 2886), '5': (3016, 3035)}\n", + "🧬cds_list = ['nngccttaccaccaagctttttccacacagccatccagactgaggaagacccggaaacttaggggccacgtgagccacggccacggccgcataggcaagcaccggaagcaccccggcggccgcggtaatgctggtggtctgcatcaccaccggatcaacttcgacaaataccacccaggctactttgggaaagttggtatgaagcattaccacttaaagaggaaccagagcttctgcccaactgtcaaccttgacaaattgtggactttggtcagtgaacagacacgggtgaatgctgctaaaaacaagactggggctgctcccatcattgatgtggtgcgatcgggctactacaaagttctgggaaagggaaagctcccaaagcagcctgtcatcgtgaaggccaaattcttcagcagaagagctgaggagaagattaagagtgttgggggggcctgtgtcctggtggcttga']\n" + ] + } + ], + "source": [ + "from data.sequence import SEQ # в этой переменной записана последовательность гена рибосомального белка\n", + "a = run_genscan(sequence=SEQ, sequence_file=None, organism=\"Vertebrate\", exon_cutoff=0.5, sequence_name=\"\")\n", + "print(a)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "GenscanOutput🦠 :\n", + "Sequence name is: Arabidopsis thaliana chromosome 3\n", + "📬Request status code is <200>\n", + "📈exon_list= {'Intr_1': (124, 279), 'Intr_2': (570, 632), 'Intr_3': (844, 925), 'Intr_4': (994, 1112), 'Intr_5': (1196, 1291), 'Intr_6': (1415, 1525), 'Intr_7': (2018, 2065), 'Intr_8': (2151, 2366), 'Term_9': (2444, 2560), 'PlyA_10': (2674, 2679)}\n", + "📉intron_list= {'1': (1, 123), '2': (280, 569), '3': (633, 843), '4': (926, 993), '5': (1113, 1195), '6': (1292, 1414), '7': (1526, 2017), '8': (2066, 2150), '9': (2367, 2443), '10': (2561, 2673)}\n", + "🧬cds_list = ['ttttcctttgaagaagagagagacgtcacaaaggaaagcactcaatcgtcatttcctcggtatctcctcctcgaaagaagaatgcctggaattagaggtccttcggaatactcgcaggaaccacctcgtcacccttctctcaaggtcaacgccaagctactccgtcacccttactggattgatccagaacccgagaaagctctttatcaaagacatcagtgtgccggtaacagaaggactgccatgagcaaagttaggaatgttagaggtgttggatgggatgtttctgctattggcaacgctgtctggggtggggcgaaactggccgatgttcttgagcttgtggggataccaaagctgactgcttctaccaatttaggagccagacatgttgagttcgttagtgttgatcgctgtaaggaggaaaatgggggcccttataaggcgtcaatcactctaagtcaagccacaaatcctgaagcggatgttctactcgcttatgagatgaatggagagaccctgaacagggatcacggatttccgttaagggtggttgtccctggtgtgattggtgctcgttcggtcaaatggcttgattccatcaatgtcatcgctgaagaaagccagagtgcaatctgctctgtggaggatgtgcaaatggtgaagcctggaaaggtaagtatcaaaggatatgcggtttcaggaggtggacgcgggatagaaagagtggacatatccctggatggaggcaaaaactgggtggaagcttctagaacgcaggaaccaggaaagcagtacatctcagaacacagctccagtgacaaatgggcatgggtgttgtttgaagccaccattgatgtttcacagactacagaggtcatcgccaaagcggttgattcggcggcgaatgttcaaccggaaaatgtggagtcggtgtggaacctaagaggggttctcaacacttcgtggcaccgtgtccttctccggcttggccactctaacttgtag']\n" + ] + } + ], + "source": [ + "file= 'data/test.fasta'\n", + "a = run_genscan(sequence=None, sequence_file=file, organism=\"Arabidopsis\", exon_cutoff=0.5, sequence_name=\"Arabidopsis thaliana chromosome 3\")\n", + "print(a)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 2. `Class OpenFasta`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Выведем все записи из фаста файла _example_fasta.fasta_ и последнюю запись из него" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "All records in fasta:\n", + "FastaRecord(id= >GTD323452, description= 5S_rRNA NODE_272_length_223_cov_0.720238:18-129(+), ACGGCCATAGGACTTTGAAAGCACCGCATCCCGTCCGATCTGCGAAGTTAACCAAGATGCCGCCTGGTTAGTACCATGGTGGGGGACCACATGGGAATCCCTGGTGCTGTG)\n", + "FastaRecord(id= >GTD678345, description= 16S_rRNA NODE_80_length_720_cov_1.094737:313-719(+), TTGGCTTCTTAGAGGGACTTTTGATGTTTAATCAAAGGAAGTTTGAGGCAATAACAGGTCTGTGATGCCCTTAGATGTTCTGGGCCGCACGCGCGCTACACTGAGCCCTTGGGAGTGGTCCATTTGAGCCGGCAACGGCACGTTTGGACTGCAAACTTGGGCAAACTTGGTCATTTAGAGGAAGTAAAAGTCGTAACAAGGT)\n", + "FastaRecord(id= >GTD174893, description= 16S_rRNA NODE_1_length_2558431_cov_75.185164:2153860-2155398(+), TTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAACGGTAACAGGAAACAGCTTGCTGTTTCGCTGACGAGTGGGAAGTAGGTAGCTTAACCTTCGGGAGGGCGCTTACCACTTTGTGATTCATGACTGGGGTGAAGTCGTAACAAGGTAACCGTAGGGGAACCTGCGGTTGGATCACCTCCTT)\n", + "FastaRecord(id= >GTD906783, description= 16S_rRNA NODE_1_length_2558431_cov_75.185164:793941-795479(-), TTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAACGGTAACAGGAAACAGCTTGCTGTTTCGCTGACGAGTGGGAAGTAGGTAGCTTAACCTTCGGGAGGGCGCTTACCACTTTGTGATTCATGACTGGGGTGAAGTCGTAACAAGGTAACCGTAGGGGAACCTGCGGTTGGATCACCTCCTT)\n", + "FastaRecord(id= >GTD129563, description= 16S_rRNA NODE_4_length_428221_cov_75.638017:281055-282593(-), CGGACGGGTGAGTAATGTCTGGGAAACTGCCTGATGGAGGGGGATAACTACTGGAAACGGTAGCTAATACCGCATAACGTCGCAAGACCAAAGAGGGGGACCGAAGTAGGTAGCTTAACCTTCGGGAGGGCGCTTACCACTTTGTGATTCATGACTGGGGTGAAGTCGTAACAAGGTAACCGTAGGGGAACCTGCGGTTGGATCACCTCCTT)\n", + "Last record:\n", + "FastaRecord(id= >GTD129563, description= 16S_rRNA NODE_4_length_428221_cov_75.638017:281055-282593(-), CGGACGGGTGAGTAATGTCTGGGAAACTGCCTGATGGAGGGGGATAACTACTGGAAACGGTAGCTAATACCGCATAACGTCGCAAGACCAAAGAGGGGGACCGAAGTAGGTAGCTTAACCTTCGGGAGGGCGCTTACCACTTTGTGATTCATGACTGGGGTGAAGTCGTAACAAGGTAACCGTAGGGGAACCTGCGGTTGGATCACCTCCTT)\n" + ] + } + ], + "source": [ + "from bio_files_processor import OpenFasta\n", + "\n", + "\n", + "print('All records in fasta:')\n", + "with OpenFasta('data/example_fasta.fasta') as my_fasta:\n", + " for i in my_fasta:\n", + " print(i)\n", + "print('Last record:')\n", + "with OpenFasta('data/example_fasta.fasta') as my_fasta:\n", + " print(my_fasta.read_records()[-1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 3. `BiologicalSequences`" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "from main_script import BiologicalSequence, NucleicAcidSequence, DNASequence, RNASequence, AminoAcidSequence" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Посмотрим на работу с последовательностями нуклеиновых кислот на примере `DNASequences`" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "mRNA = AUCGAUCG\n", + "GC_content = 50.0\n", + "complement_strain = TAGCTAGC\n", + "object type of compliment: \n", + "object type of transcripte: \n" + ] + } + ], + "source": [ + "dna = DNASequence('ATCGATCG')\n", + "print('mRNA = ', dna.transcribe())\n", + "print('GC_content =', dna.gc_content())\n", + "print('complement_strain = ', dna.complement())\n", + "print('object type of compliment:', type(dna.complement()))\n", + "print('object type of transcripte:', type(dna.transcribe()))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "ename": "UnexpectedSymbolInSeqError", + "evalue": "Unexpected symbol 'U' in the sequence", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mUnexpectedSymbolInSeqError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[7], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m bad_dna \u001b[38;5;241m=\u001b[39m \u001b[43mDNASequence\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mAUCGATC\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/HW18/main_script.py:243\u001b[0m, in \u001b[0;36mDNASequence.__init__\u001b[0;34m(self, seq)\u001b[0m\n\u001b[1;32m 242\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__init__\u001b[39m(\u001b[38;5;28mself\u001b[39m, seq):\n\u001b[0;32m--> 243\u001b[0m \u001b[38;5;28;43msuper\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[38;5;21;43m__init__\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mseq\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/HW18/main_script.py:205\u001b[0m, in \u001b[0;36mNucleicAcidSequence.__init__\u001b[0;34m(self, seq)\u001b[0m\n\u001b[1;32m 203\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mseq \u001b[38;5;241m=\u001b[39m seq\n\u001b[1;32m 204\u001b[0m \u001b[38;5;28msuper\u001b[39m()\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__init__\u001b[39m(seq)\n\u001b[0;32m--> 205\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcheck_alphabet\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/HW18/main_script.py:155\u001b[0m, in \u001b[0;36mBiologicalSequence.check_alphabet\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 153\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mseq:\n\u001b[1;32m 154\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m i \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mget_alphabet():\n\u001b[0;32m--> 155\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m UnexpectedSymbolInSeqError(i)\n", + "\u001b[0;31mUnexpectedSymbolInSeqError\u001b[0m: Unexpected symbol 'U' in the sequence" + ] + } + ], + "source": [ + "bad_dna = DNASequence('AUCGATC')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Посмотрим возможности `AminoAcidSequence`" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "M.W. = 639.81\n", + "Peptide length = 6\n" + ] + } + ], + "source": [ + "peptide = AminoAcidSequence('VALINE')\n", + "print('M.W. = ', peptide.molecular_weight())\n", + "print('Peptide length = ', len(peptide))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "ename": "UnexpectedSymbolInSeqError", + "evalue": "Unexpected symbol 'X' in the sequence", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mUnexpectedSymbolInSeqError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[9], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m wrong_peptide \u001b[38;5;241m=\u001b[39m \u001b[43mAminoAcidSequence\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mVALIXX\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/HW18/main_script.py:181\u001b[0m, in \u001b[0;36mAminoAcidSequence.__init__\u001b[0;34m(self, seq)\u001b[0m\n\u001b[1;32m 179\u001b[0m \u001b[38;5;28msuper\u001b[39m()\u001b[38;5;241m.\u001b[39m\u001b[38;5;21m__init__\u001b[39m(seq)\n\u001b[1;32m 180\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mseq \u001b[38;5;241m=\u001b[39m seq\n\u001b[0;32m--> 181\u001b[0m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcheck_alphabet\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/HW18/main_script.py:155\u001b[0m, in \u001b[0;36mBiologicalSequence.check_alphabet\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 153\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mseq:\n\u001b[1;32m 154\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m i \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mget_alphabet():\n\u001b[0;32m--> 155\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m UnexpectedSymbolInSeqError(i)\n", + "\u001b[0;31mUnexpectedSymbolInSeqError\u001b[0m: Unexpected symbol 'X' in the sequence" + ] + } + ], + "source": [ + "wrong_peptide = AminoAcidSequence('VALIXX')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ML", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/bio_files_processor.py b/bio_files_processor.py index b65fe39..4dc12ae 100644 --- a/bio_files_processor.py +++ b/bio_files_processor.py @@ -1,13 +1,73 @@ import os +from dataclasses import dataclass +@dataclass +class FastaRecord: + id: str + seq: str + description: str -def convert_multiline_fasta_to_oneline(input_fasta, output_fasta = ''): + def __repr__(self): + return f'FastaRecord(id= {self.id}, description= {self.description}, {self.seq})' + + +class OpenFasta: + def __init__(self, file_path: str): + self.file_path = file_path + self.file = open(file_path, 'r') + self.title = None + self.stop = False + + def __enter__(self): + return self + + def __exit__(self, exc_type, exc_value, traceback): + self.file.close() + + def __next__(self): + if self.stop: + raise StopIteration + if not self.title: + self.title = self.file.readline().strip().split(' ', 1) + id_cur = self.title[0] + description = self.title[1] + seq_parts = [] + line = self.file.readline().strip() + while line and not line.startswith('>'): + seq_parts.append(line) + line = self.file.readline().strip() + seq = ''.join(seq_parts) + res = FastaRecord(id_cur, seq, description) + self.title = line.split(' ', 1) + id_cur = self.title[0] + if id_cur == '': + self.stop = True + return res + + def read_records(self): + records = [] + while not self.stop: + record = self.read_record() + if record: + records.append(record) + else: + break + return records + + def __iter__(self): + return self + + def read_record(self): + return self.__next__() + + +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 + :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'): + 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') @@ -26,7 +86,7 @@ def convert_multiline_fasta_to_oneline(input_fasta, output_fasta = ''): output_file.write(line) current = [] break - + for line in input_file: if line.startswith('>'): output_file.write(''.join(current) + '\n') @@ -35,9 +95,9 @@ def convert_multiline_fasta_to_oneline(input_fasta, output_fasta = ''): else: current.append(line.strip()) output_file.write(''.join(current) + '\n') - -def change_fasta_start_pos(input_fasta, n: int, output_fasta=''): + +def change_fasta_start_pos(input_fasta, n = 0, output_fasta=''): """ This function moves the position to a new starting nucleotide. :param input_fasta: data for processing (1-string single fasta) @@ -62,7 +122,8 @@ def change_fasta_start_pos(input_fasta, n: int, output_fasta=''): s2 = s2 + string[i] else: s1 = s1 + string[i] - new_fasta = s1 +s2 + new_fasta = s1 + s2 with open(output_fasta, 'w') as output_file: - output_file.write(name +'\n') - output_file.write(new_fasta[0] + '\n') + output_file.write(name + '\n') + output_file.write(new_fasta + '\n') + diff --git a/data/example_fasta.fasta b/data/example_fasta.fasta new file mode 100644 index 0000000..ae5498d --- /dev/null +++ b/data/example_fasta.fasta @@ -0,0 +1,18 @@ +>GTD323452 5S_rRNA NODE_272_length_223_cov_0.720238:18-129(+) +ACGGCCATAGGACTTTGAAAGCACCGCATCCCGTCCGATCTGCGAAGTTAACCAAGATGCCGCCTGGTTAGTACCATGGTGGGGGACCACATGGGAATCCCT +GGTGCTGTG +>GTD678345 16S_rRNA NODE_80_length_720_cov_1.094737:313-719(+) +TTGGCTTCTTAGAGGGACTTTTGATGTTTAATCAAAGGAAGTTTGAGGCAATAACAGGTCTGTGATGCCCTTAGATGTTCTGGGCCGCACGCGCGCTACACT +GAGCCCTTGGGAGTGGTCCATTTGAGCCGGCAACGGCACGTTTGGACTGCAAACTTGGGCAAACTTGGTCATTTAGAGGAAGTAAAAGTCGTAACAAGGT +>GTD174893 16S_rRNA NODE_1_length_2558431_cov_75.185164:2153860-2155398(+) +TTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAACGGTAACAGGAAACAGCTTGCTGTTTCGCTGACGAGTGG +GAAGTAGGTAGCTTAACCTTCGGGAGGGCGCTTACCACTTTGTGATTCATGACTGGGGTGAAGTCGTAACAAGGTAACCGTAGGGGAACCTGCGGTTGGATC +ACCTCCTT +>GTD906783 16S_rRNA NODE_1_length_2558431_cov_75.185164:793941-795479(-) +TTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAACGGTAACAGGAAACAGCTTGCTGTTTCGCTGACGAGTGG +GAAGTAGGTAGCTTAACCTTCGGGAGGGCGCTTACCACTTTGTGATTCATGACTGGGGTGAAGTCGTAACAAGGTAACCGTAGGGGAACCTGCGGTTGGATC +ACCTCCTT +>GTD129563 16S_rRNA NODE_4_length_428221_cov_75.638017:281055-282593(-) +CGGACGGGTGAGTAATGTCTGGGAAACTGCCTGATGGAGGGGGATAACTACTGGAAACGGTAGCTAATACCGCATAACGTCGCAAGACCAAAGAGGGGGACC +GAAGTAGGTAGCTTAACCTTCGGGAGGGCGCTTACCACTTTGTGATTCATGACTGGGGTGAAGTCGTAACAAGGTAACCGTAGGGGAACCTGCGGTTGGATC +ACCTCCTT \ No newline at end of file diff --git a/data/sequence.py b/data/sequence.py new file mode 100644 index 0000000..3942d8e --- /dev/null +++ b/data/sequence.py @@ -0,0 +1 @@ +SEQ = " CTTTTTCGTCTGGGCTGCCAACATGGTAGGTGTTTCGTTTCTTGCCTCCTCTTCCTTGCCGGCGGAGACCCCTAAGCTGTATTCCCATTGCCCCTAGTCATCCACTCCCTACCATGGTCGGGGCTTCCAGGCTGCGCATGGCCGCCTGCGGGGCAGGGTGGCCGGCGCGGGCCCGGGGCGGGGCTCCCGGAGCCGTGTGTTAGGCCCGCGGTTCGGATCTCTAGGACACGCGGGCCCCTGCGCTACCGTGGTGAGACCTCACGGCCCTGAGCGGATCGGTACCCTCAGCTTTCCCAAACGCTCCAGAAGTTAGGTCTTTGACCCACAGGCTTACAGGACCATCTCGGCTGGCGGGCATCGCCCCCTGCCCCTAATTCCTTAGGCCTTACCACCAAGCTTTTTCCACACAGCCATCCAGACTGAGGAAGACCCGGAAACTTAGGGGCCACGTGAGCCACGGCCACGGCCGCATAGGTAAGTGCCGGCTTCCCCTCGGGGTGGGCCTTGGGCTCTCTTCGGGTGCTTAGCTAGTCTGGAGATCGGTAGCCTATAAGTGGGTTAGAATAAGACCTTTTTGTGGTCAAGTTGCACAGCTGTTGATTTTTTTCTGACGATCCTCTAGTATTCCAGTTCTAAGGAATTTCACATCAGTGGGGTAATAGGAATTGAGCAGGCACGGTATTGGGTTAGTTGAAGACATGGAGTACTGTGGGAATGCTGTGATGTGGAACCTGAAAAGATGTTTCACCCGGAATCCTAAAGTAATCGCATTGCTGAAAACCGGCATCGGTAGGGTGGGAACAGCGTAAGCGGGACACAGAAGTCTGGGAAACACTCTGCTTTTGTGCGAGGAAGTATTGAGATGCATGAGAAGGCTGTGTGGTGCATGTAGCTTTTTTGTGTGTGTGTGAGACTGATCACTGTCGCCCAGGCTGGAGTGCAGTGGCGAAATCTCGGCTCACTGCAACCCCTGCCTGCCGGGCTCAAGCGATTCTCCTGCCTCAGCCTCCCGAGTAGCTGATATTACAGGTGTGCGCCACGACGCCCGGCTAATTCTTTTTCTATTTTTAGTAGAGTCGGGGGTTTCTCCGTGTTGGCCAGGCTGGTCTCGAACTCCTGACCTCAGGTGATCTACCCCCCCTCGTCCTCCTAAAATGCTGGGATTACAGGCATGAGCCATCACACCGGCCCCACGTAGCTTTGTATTCCTGCAGGCAAGCACCGGAAGCACCCCGGCGGCCGCGGTAATGCTGGTGGTCTGCATCACCACCGGATCAACTTCGACAAATAGTAAGTGTCCTTGGACTGCTTTTATTGACACAGCTTGGGAGGTAGGGGCAGAGAGAGGGCTGGCTTAAACAAAAAGTTTAGAAGCAAGCCTTGCCTATTGCTGTTTTTTACCAAGTTAACACTTGGTGTGAACTGAGAACCTGTCATCGAGGCTAGAGTCACGCTTGGGTATCGGCTATTGCCTGAGTGTGCTAGAGTCCTCGAAGAGTAACTGCTGACCTTATTCACTGGCTGTGGGCCTTATGGCACAGTCAGTCACCAGGTTAGAGACATGCTTCACATTCACCTACCCACAAACTAGTGGATGATAAATTTTGGCTATTCAGAAGACGTTTATTATAGGAGTATGTAGATTTTCCATAGAGTGCTGTTATGTGACTTGAATTTTAGTCTCGGCCCTGCCTCTGACATTGTCGGTGGTTTATCCTGGTTCCAGGAAATAAGACTAGCCTTTTCCTCATGATAGTCTTTGGTGGTTTTTAAAACAGTTGTTTAAGTCAACAGATGTATCATATGCCTGACACTGCTCTACACCAGTGAATAATTTACACTCTAATAGGGGGTGGTAACTATAAAGATGATAAACATAGCATCTTAATTGGAGTGTGTATGAAGGTGGTTGTTACCTCTTCCTAGCCACCCAGGCTACTTTGGGAAAGTTGGTATGAAGCATTACCACTTAAAGAGGAACCAGAGCTTCTGCCCAACTGTCAACCTTGACAAATTGTGGACTTTGGTCAGTGAACAGACACGGGTGAATGCTGCTAAAAACAAGACTGGGGCTGCTCCCATCATTGATGTGGTGCGATCGGTAAGTTAATTGGATGTTTTTCTGTACTTCCATACCTTCCCTTACAAAACTCTGGCTTAATCTAATCCACTTATATAATCTGTACTTCCCAGTTACCTACCAGACATTGATATTCTTCCTGTGGTAGAATTATCATAGGTAGTTCCCTATCCGTAGCAGTGCCTACTGTCACTGCCCAGGTTGTATCAGGTTTGCATTTCGTGCTTGAACTATAGCTGGTTTTCACTGAGCACAGCTCTTGGCCCTTCATGTTCTCCAGATAATAGAATCCTAATATGTTCCATTGATACTCAGTGCCATGCATTATCTGAAGAGATTTTCCCCCAAAACAGATGTATTATGTCTGTCCTTGCGGGGGTTCTGGTCCCTGTGTCAGTCTTAACTCTCATGAATATAGAGGTAGTGTTAAGAGGCCAGAACCCTAGGGACGCTTTAAATTCACTTCCCAGCCTATTTAATGTCCATTGAGTAGTTCTGGTGGTCAGGAAGGTGGTTGTCTTCTTTTGCTTAGCAGGGGGTATTTGAGCAGGAGGAGGCTTATGCTTTGCCGAGACTAGAGTCACATCCTGACACAACTCTTGTCCTGGTGTGCTAGAGTACTCGAAGAGAATCTACTGGTCTTGATTCACTGGTGGGGGCAGTCGGTGCCCCCGTTAGTGCCCAGATCAGAAACATACATACCCTGCCTAGGGATTTAGAAAGTGGGTTGGCAGTCTTTCCTCACGCCCATCACGCAGTTGGTACCTACTACAGTGTATTGTAAACTTTTTTCTCTGTTCTTCTAGGGCTACTACAAAGTTCTGGGAAAGGGAAAGCTCCCAAAGCAGCCTGTCATCGTGAAGGCCAAATTCTTCAGCAGAAGAGCTGAGGAGAAGATTAAGAGTGTTGGGGGGGCCTGTGTCCTGGTGGCTTGAAGCCACATGGAGGGAGTTTCATTAAATGCTAACTACTTTTTCCTTGTGGTGTGAGTGTAGGTTCTTCAGTGGCACCTCTACATCCTGTGTGCATTGGGAGCCCAGGTTCTAGTACTTAGGGTATGAAGACATGGGGTCCTCTCCTGACTTCCCTCAAATATATGGTAAACGTAAGACCAACACAGACGTTGGCCAGTTAAACATTTCTGTTTATAAAGTCAGAATAATACCTGTTGATCACTGAAAGGCCTGCATGTATTGTACTCTGAATTTTACAGTGAATGAGAGAATGTACCCTAATTGTTCAACAGGGCTCAAAAGGAAAGATTCCATTTTGATGGGTCACATTCTAAAGAGGGGCAGTGTGATAGGAATGAGATGGTCCTTTAGGACTTAAGTTCTCAGCCCAAGGTTTTTCCACGTGGCCCCCTCATCTTTTTTTTTTTTTTAAACGGAGTCTCTCTTGCCAGGCTGGAGTGCAGTGGCACGATCTCGGCTCACTGCAGCCTCCGCCTCCCAGGTTAAGCGATTCTCCTGCCTCAGCTTCCTGACTAACTGGGATTACAGGCGCCCACCACCATGCCCAGCTAATTTTTGTATTTTCAGTAGAGATGGGGTTTCACCATGTTGGCCATGCTGGTCTCTAACTCCTAACCTCAAGTGATCTGCCCACATCGGCCTCCAAAAGTTCTGGGATTATAGTGTGAGCCACTGCGCCCGGCCATGGCTCCTTAATCTTGATCCAAATTATTGTTACATCCAGAATGTGATGAATCAAAATCTCGAGATGGGGGTCCAGCAATCTGAAATTTCAGTATGCCAGGGCTTTTCTGTATGTCAAAGTGGGTTTGAAATAGTTAATTTTTCTTCTAGTCTGAAATGTATCGGGAAAATTTGGAAATCCTGAAGGCTGGAAATTGAAATAAGTTTTTCTAGGATTTGTGTCTCTTGCTATTGGAAAACTGATGGTGACCA" \ No newline at end of file diff --git a/data/test.fasta b/data/test.fasta new file mode 100644 index 0000000..d0a7f96 --- /dev/null +++ b/data/test.fasta @@ -0,0 +1,42 @@ +>NC_003074.8:c317431-314693 Arabidopsis thaliana chromosome 3, partial sequence +TGTTTGGGCTCAATTTGGGCTCTGTCTGTTATCCACAAATTAAGAGACGCAAACCTAAAACGACGCCGTT +CCATAGTTTTCCTTTGAAGAAGAGAGAGACGTCACAAAGGAAAGCACTCAATCGTCATTTCCTCGGTATC +TCCTCCTCGAAAGAAGAATGCCTGGAATTAGAGGTCCTTCGGAATACTCGCAGGAACCACCTCGTCACCC +TTCTCTCAAGGTCAACGCCAAGGTAAGCTTCTTCTTCTTCTTCTTCTCTCTTTCCAATCGAATCAGTCAC +TAGTCTGATATAGCCTTTTTCTTCATGGCATCTAAATTTGGGTGATCCGCAGGAACCTTTCAACGCCGAG +CCTCCCCGCTCGGCCTTAGTCTCATCTTATGTCACTCCCGTCGACCTTTTCTACAAGCGAAATCATGGTC +CCATCCCCATTGTTGATCACCTTCAAAGGTCTCTCTCTCACTCTCATTCATCTTCAATTTTGCCTCCTCC +GCTCATTTCTGCTATTTCTTTTGCCAAACCAGCTACTCCGTCACCCTTACTGGATTGATCCAGAACCCGA +GAAAGCTCTTTATCAAAGACATCAGGTACCTTTATTTATTTGCCTCTCTTTACATTCTCCAACTTTGATC +TTGTATCTCTCCATCAAGCACCACTCTTGTTTCCTCCTCAAAATTGCAGGTCCCTCCCAAAGTACAATGT +TACTGCTACTCTTCAGGTTCATTCTTCCTTTGAGTAAATCTTCTACCTTTAGTTTTTCTCCACATTTAAT +TTAACATCACTTTTGTTTCAATGCAGTGTGCCGGTAACAGAAGGACTGCCATGAGCAAAGTTAGGAATGT +TAGAGGTGTTGGATGGGATGTTTCTGCTATTGGCAACGGTTCGTCTCTCTTTCTCCCATTCCCATGCTTT +CAATACATACTTGTCTTGAAAGAATTTGCCTTTCAGCTGTCTGGGGTGGGGCGAAACTGGCCGATGTTCT +TGAGCTTGTGGGGATACCAAAGCTGACTGCTTCTACCAATTTAGGAGCCAGACATGTTGAGTTCGTTAGT +GTTGATCGCTGTAAGGTAGATATATATTAATAAATTCCACCTTTTTATGCTTTGTTATGAATTCCCATGA +AAAGTCTATGCTCTTTCTATATCTCTAGGAGGAAAATGGGGGCCCTTATAAGGCGTCAATCACTCTAAGT +CAAGCCACAAATCCTGAAGCGGATGTTCTACTCGCTTATGAGATGAATGGAGAGGTTTGATATCGCTCTC +CTAAATTAGTCTCACTTGTGTTTAAATCATCACTTGTTTGACATGTGAGTATGAGTGAGTACTCGGAAAT +GTGCTTATTTGTTGGAAAACGTGTTTCTTCTGTGAAGACCCTGAACAGGGATCACGGATTTCCGTTAAGG +GTGGTTGTCCCTGGTGTGATTGGTGCTCGTTCGGTCAAATGGCTTGATTCCATCAATGTCATCGCTGAAG +AAAGCCAGGTCTCAACTAAAGTGTCCGGCTTTTACCCTCTTTTTAGCAATCAGAAACTGTTCTAACATTG +TAAGGGACTTGGATCAAAATGTGATGTAGGGATTCTTCATGCAAAAAGATTACAAAATGTTTCCACCCTC +TGTCAATTGGGATAATATCAACTGGTCCTCTAGGAGACCGCAAATGGATTTCCCTGTTCAGGTAAATACA +ACTTTCCCTTCTTTCTTACAGAATACTTGACTTTGTAATGACGTTTTTCGTTTATATGCGCTACAAAATT +AAGATTTTCCAATTTTCTTGCGAAGTTCGTAATTCTGTATTCTTCCATGAAAACAGTTGATGATGGTGAT +GAATCGACTGATCATATTTTCTACTGGTATACTGAAAATTTTACTTGTTTGGACCAGTTGCCTTATGAAT +TTGATCAGATAGTAGAATTTGTAAATGAACATTAGTTTGTCTTGAGTCTTGACCCCATACTCGTTTATCA +TTGTTAACAGAGTGCAATCTGCTCTGTGGAGGATGTGCAAATGGTGAAGCCTGGAAAGGTTAGAACGCGG +AATTCACTACCTAATCTGATCTGACCAATGTAGGAAAGCAAGAGCTAATTCAATTAAAAAAAAAAAATTG +CAGGTAAGTATCAAAGGATATGCGGTTTCAGGAGGTGGACGCGGGATAGAAAGAGTGGACATATCCCTGG +ATGGAGGCAAAAACTGGGTGGAAGCTTCTAGAACGCAGGAACCAGGAAAGCAGTACATCTCAGAACACAG +CTCCAGTGACAAATGGGCATGGGTGTTGTTTGAAGCCACCATTGATGTTTCACAGACTACAGAGGTCATC +GCCAAAGCGGTAAAGCCTCAACATCTTGCTCAGCTAATATAAATGGCTTCTGACAAGAGTTAATAGAAAA +CTGTGGAAATGTGCAGGTTGATTCGGCGGCGAATGTTCAACCGGAAAATGTGGAGTCGGTGTGGAACCTA +AGAGGGGTTCTCAACACTTCGTGGCACCGTGTCCTTCTCCGGCTTGGCCACTCTAACTTGTAGACCAAGT +TACCACTATGCGTCTGGTCTCTCTTTATCTTTTCTTTAAGAACTAAAACATCCTTGAATATGCTTGTCCT +TGTGGATCTGTGAGAGCAAAATACTGTTTATAGTTAAATAAAATGGGAAAACCCATCTGTTGTCATGTGG +TGGATTATAAATTATTTATTTATATGATAGAACTGTTATTCATTATTATATATTAATGAAATCATATGGT +TGAGCAAAT + diff --git a/draft_custom_ran_forest.ipynb b/draft_custom_ran_forest.ipynb new file mode 100644 index 0000000..446bca9 --- /dev/null +++ b/draft_custom_ran_forest.ipynb @@ -0,0 +1,1019 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 76, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import multiprocessing as mp\n", + "import seaborn as sns\n", + "import random\n", + "import pandas as pd\n", + "from sklearn import datasets\n", + "from sklearn.base import BaseEstimator\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.ensemble import (RandomForestClassifier,\n", + " ExtraTreesClassifier,\n", + " VotingClassifier)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "metadata": {}, + "outputs": [], + "source": [ + "SEED = 111\n", + "random.seed(SEED)\n", + "np.random.seed(SEED)" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "metadata": {}, + "outputs": [], + "source": [ + "#пример из ML\n", + "X, y = datasets.make_moons(n_samples=500, noise=0.30, random_state=SEED)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "metadata": {}, + "outputs": [], + "source": [ + "class RandomForestClassifierCustom(BaseEstimator):\n", + "\n", + " def __init__(\n", + " self, n_estimators=10, max_depth=None, max_features=None, random_state=SEED, n_jobs = 1\n", + " ):\n", + " self.n_estimators = n_estimators\n", + " self.max_depth = max_depth\n", + " self.max_features = max_features\n", + " self.random_state = random_state\n", + " self.random_state_initial = random_state\n", + " self.n_jobs = n_jobs\n", + " self.trees = []\n", + " self.feat_ids_by_tree = []\n", + "\n", + " def fit_single_tree(self, i):\n", + " self.random_state = np.random.seed(self.random_state_initial + i)\n", + " indices_features = np.random.choice(range(X.shape[1]), size=self.max_features, replace=False)\n", + " pseudo_X_indices = np.random.choice(range(X.shape[0]), size=X.shape[0], replace=True)\n", + " pseudo_X = X[pseudo_X_indices]\n", + " pseudo_y = y[pseudo_X_indices]\n", + " tree = DecisionTreeClassifier(max_depth=self.max_depth, random_state=self.random_state)\n", + " tree.fit(pseudo_X[:, indices_features], pseudo_y)\n", + " return tree, indices_features \n", + "\n", + " def fit(self, X, y, n_jobs = 1):\n", + " self.classes_ = sorted(np.unique(y))\n", + "\n", + " with mp.Pool(processes=self.n_jobs) as pool:\n", + " results = pool.map(self.fit_single_tree, range(self.n_estimators))\n", + "\n", + " self.trees, self.feat_ids_by_tree = zip(*results)\n", + " return self\n", + "\n", + "# TODO: add n_jobs for predict methods\n", + " def predict_proba(self, X):\n", + " predictions = np.zeros((X.shape[0], len(self.classes_))) \n", + "\n", + " for i in range(0, self.n_estimators-1):\n", + " tree = self.trees[i]\n", + " feat_ids = self.feat_ids_by_tree[i]\n", + " prediction = tree.predict(X[:, feat_ids])\n", + "\n", + " for j in range(X.shape[0]):\n", + " pred_class = prediction[j]\n", + " predictions[j, pred_class] += 1 \n", + "\n", + " return predictions / self.n_estimators\n", + " \n", + " def predict(self, X):\n", + " probas = self.predict_proba(X)\n", + " predictions = np.argmax(probas, axis=1)\n", + " \n", + " return predictions" + ] + }, + { + "cell_type": "code", + "execution_count": 80, + "metadata": {}, + "outputs": [], + "source": [ + "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=SEED)" + ] + }, + { + "cell_type": "code", + "execution_count": 107, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 39.8 ms, sys: 41 ms, total: 80.8 ms\n", + "Wall time: 253 ms\n" + ] + }, + { + "data": { + "text/html": [ + "
RandomForestClassifierCustom(max_depth=10, max_features=2, n_estimators=40)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
" + ], + "text/plain": [ + "RandomForestClassifierCustom(max_depth=10, max_features=2, n_estimators=40)" + ] + }, + "execution_count": 107, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%%time\n", + "clf = RandomForestClassifierCustom(n_estimators=40, max_depth=10, max_features=2, random_state=SEED, n_jobs = 1)\n", + "clf.fit(X_train, y_train)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 109, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 29.1 ms, sys: 44.9 ms, total: 74 ms\n", + "Wall time: 153 ms\n" + ] + }, + { + "data": { + "text/html": [ + "
RandomForestClassifierCustom(max_depth=10, max_features=2, n_estimators=40,\n",
+       "                             n_jobs=2)
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
" + ], + "text/plain": [ + "RandomForestClassifierCustom(max_depth=10, max_features=2, n_estimators=40,\n", + " n_jobs=2)" + ] + }, + "execution_count": 109, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%%time\n", + "clf = RandomForestClassifierCustom(n_estimators=40, max_depth=10, max_features=2, random_state=SEED, n_jobs = 2)\n", + "clf.fit(X_train, y_train, n_jobs = 2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Пока все странно, возможно надо работать с данными, на которых fit идет долго, чтобы все было стабильнее" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ML", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/main_script.py b/main_script.py new file mode 100644 index 0000000..fdc394e --- /dev/null +++ b/main_script.py @@ -0,0 +1,348 @@ +from abc import ABC, abstractmethod +from Bio import SeqIO +from Bio import SeqUtils +import os +import sys +import io +import datetime +import requests +import dotenv +import re +from dataclasses import dataclass +from typing import List, Dict +from bs4 import BeautifulSoup + + +def telegram_logger(chat_id: str): + """ + Decorator that logs the execution of the decorated function and sends a Telegram message. + Message content: + Time of function execution, if success + Error type and Error Message in case of error + + Input Args: + chat_id (str): The ID of the Telegram chat where the message will be sent. + """ + + def decorator(func): + def inner_func(*args, **kwargs): + result = None + start_time = datetime.datetime.now() + try: + result = func(*args, **kwargs) + x_time = datetime.datetime.now() - start_time + message = f"🦥 Функция {func.__name__} выполнена успешно.\nВремя работы:\n*{x_time:}*" + except Exception as err: + message = f"👿 Функция {func.__name__} выполнена с ошибкой *{type(err).__name__}*\n{str(err)}" + send_message(chat_id, message) + return result + return inner_func + + return decorator + pass + + +def send_message(chat_id: str, text: str): + ''' + Support function for telegram_logger + Send telegram message for telegram bot with appropriate text + ''' + dotenv.load_dotenv() + token = os.getenv("TG_API_TOKEN") + url = f'https://api.telegram.org/bot{token}/sendMessage' + params = {'chat_id': chat_id, 'text': text, 'parse_mode': 'Markdown'} + response = requests.post(url, params) + return response + + +@dataclass +class GenscanOutput: + """ Represents the output of the GENSCAN analysis. + (http://argonaute.mit.edu/GENSCAN.html) + + Attributes: + - status: str: The status of the analysis, either "Successful" or "Failed". + - cds_list: List[str]: List of predicted peptide sequences. + - intron_list: List[Any]: Dict of found introns, where each Key is intron number, + value is tuple (start position, end position) + - exon_list: List[Any]: Dict of found exons, where each Key is exon description and number, + value is tuple (start position, end position) + """ + status: int + name: str + intron_list: Dict[str, tuple] + exon_list: Dict[str, tuple] + cds_list: List[str] + + def __str__(self): + return f'GenscanOutput🦠 :\nSequence name is: {self.name}\n📬Request status code is <{self.status}>\n📈exon_list= {self.exon_list}\n📉intron_list= {self.intron_list}\n🧬cds_list = {self.cds_list}' + + +def run_genscan(sequence=None, sequence_file=None, organism="Vertebrate", exon_cutoff=1, sequence_name=""): + """ + Executes a GENSCAN analysis by making a request similar to filling out the form on the website. + Accepts all parameters that can be specified on the website (except Print options). + + Args: + - sequence: str, optional: The sequence as a string or any convenient data type. + - sequence_file: str, optional: The path to the file containing the sequence, which can be uploaded and used instead of the sequence. + - organism: str, optional: The organism type for the analysis. Default is "Vertebrate". + - exon_cutoff: float, optional: The exon probability cutoff. Default is 1.00. + - sequence_name: str, optional: The name of the sequence. Default is an empty string. + + Returns: + - GenscanOutput: An object containing the analysis output. + """ + url = "http://hollywood.mit.edu/cgi-bin/genscanw_py.cgi" + data = { + "-o": organism, + "-e": exon_cutoff, + "-n": sequence_name, + "-u": sequence_file, + "-s": sequence, + "-p": "Predicted CDS and peptides" + } + if sequence: + data["-s"] = sequence + elif sequence_file: + with open(sequence_file, 'r') as file: + data["-s"] = file.read() + response = requests.post(url, data=data) + soup = BeautifulSoup(response.content, 'html.parser') + lines = soup.prettify().split('\n') + status = response.status_code + + name = '---' + if sequence_name: + name = sequence_name + cds_list = [] + intron_list = {} + exon_list = {} + exon_lines = [line for line in lines if re.match(r"^ \d", line)] + i = 1 + end = 0 + for l in exon_lines: + parts = l.strip().split() + exon_list[f"{parts[1]}_{i}"] = (int(parts[3]), int(parts[4])) + intron_list[f"{i}"] = (end + 1, int(parts[3]) - 1) + end = int(parts[4]) + i += 1 + current_seq = "" + no_empty_lines = [i for i in lines if i] + for line in no_empty_lines: + if line[0] in {'a', 't', 'c', 'g', 'n'}: + current_seq += line.strip() + elif current_seq != '': + cds_list.append(current_seq) + current_seq = "" + if current_seq: + cds_list.append(current_seq) + return GenscanOutput(status=status, name=name, cds_list=cds_list, intron_list=intron_list, exon_list=exon_list) + + +class BiologicalSequence(ABC, str): + + def __init__(self, seq): + self.seq = seq + + @abstractmethod + def get_alphabet(self): + pass + + def check_alphabet(self): + for i in self.seq: + if i not in self.get_alphabet(): + raise UnexpectedSymbolInSeqError(i) + + + @abstractmethod + def to_print_seq(self): + return self.seq + pass + + +class UnexpectedSymbolInSeqError(ValueError): + def __init__(self, symbol): + self.message = f"Unexpected symbol '{symbol}' in the sequence" + super().__init__(self.message) + + +class AminoAcidSequence(BiologicalSequence): + alphabet = set('ARNDCEQGHILKMFPSTWYVUOarndceqghilkmfpstwyvuo') + 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 get_alphabet(self): + return set('ARNDCEQGHILKMFPSTWYVUOarndceqghilkmfpstwyvuo') + + + 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 round(m, 2) + +class NucleicAcidSequence(BiologicalSequence): + dict_comp = {'A': 'T', 'C': 'G', 'T': 'A', 'G': 'C', 'a': 't', 'c': 'g', 't': 'a', 'g': 'c', 'u' : 'a', 'U': 'A'} + def __init__(self, seq): + self.seq = seq + super().__init__(seq) + self.check_alphabet() + + + def get_alphabet(self): + return set('AUCGTaucgt') + + + 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) + if isinstance(self, DNASequence): + return DNASequence(res) + elif isinstance(self, RNASequence): + return RNASequence(res) + else: + return NucleicAcidSequence(res) + +class DNASequence(NucleicAcidSequence): + dict_trans = {'A': 'A', 'C': 'C', 'T': 'U', 'G': 'G', 'a': 'a', 'c': 'c', 't': 'u', 'g': 'g'} + + def get_alphabet(self): + return set('ATCGatcg') + + + 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(NucleicAcidSequence): + def __init__(self, seq): + super().__init__(seq) + def get_alphabet(self): + return set('AUCGaucg') + + +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_threshold: threshold value for average quality per nucleotide (phred33 scale) + :return:(dict) records 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_threshold=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 sequence length in number of nucleotide + quality_threshold: (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_threshold) + SeqIO.write((record for record in filtered_records_q), output_fastq, "fastq") + return + + + + diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..adae817 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,5 @@ +beautifulsoup4==4.12.3 +Bio==1.7.0 +pytest==8.2.0 +python-dotenv==1.0.1 +Requests==2.31.0 diff --git a/requirments.txt b/requirments.txt deleted file mode 100644 index 04e9ef2..0000000 --- a/requirments.txt +++ /dev/null @@ -1,2 +0,0 @@ -Python 3.12.0 -biopython==1.83 \ No newline at end of file diff --git a/test_module.py b/test_module.py new file mode 100644 index 0000000..9f6659a --- /dev/null +++ b/test_module.py @@ -0,0 +1,107 @@ +import os +import pytest +from bio_files_processor import change_fasta_start_pos +from main_script import (BiologicalSequence, DNASequence, RNASequence, NucleicAcidSequence, AminoAcidSequence, + UnexpectedSymbolInSeqError) + +class TestChangeFastaStartPos: + def test_sequence_length(self): + with open("test.fasta", "w") as f: + f.write(">TestSeq\nACGTACGTGG") + change_fasta_start_pos("test.fasta", n=2, output_fasta="shifted") + with open("test.fasta") as f: + original_seq = f.readlines()[1].strip() + with open("shifted_fasta_results/shifted.fasta") as f: + shifted_seq = f.readlines()[1].strip() + assert len(original_seq) == len(shifted_seq), 'expected the same sequence length ' + os.remove("test.fasta") + os.remove("shifted_fasta_results/shifted.fasta") + + + def test_file_content(self): + with open("test.fasta", "w") as f: + f.write(">TestSeq\nACGTACGTGG") + change_fasta_start_pos("test.fasta", n=2, output_fasta="shifted") + with open("shifted_fasta_results/shifted.fasta") as f: + content = f.readlines() + assert content[0].strip() == '>TestSeq', 'No changes expected in sequence name' + assert content[1].strip() == 'GTACGTGGAC', ' GTACGTGGAC expected ' + os.remove("test.fasta") + os.remove("shifted_fasta_results/shifted.fasta") + + +class TestDNASequences: + def test_dna_alphabet(self): + dna = DNASequence('') + target = set('ATGCatgc') + result = DNASequence.get_alphabet(dna) + assert target == result + + + def test_gc_content(self): + dna = 'ATGC' + target = 50.0 + result = NucleicAcidSequence(dna).gc_content() + assert target == result, '50.0 for ATGC expected' + + + def test_transcribe(self): + dna = 'ATGC' + target = 'AUGC' + result = DNASequence(dna).transcribe() + assert target == result, "expected 'AUGC' for 'ATGC' seq" + + + def test_complemet(self): + dna = 'ATCG' + target = 'TAGC' + result = NucleicAcidSequence(dna).complement() + assert target == result, "'TAGC' for 'ATCG' expected'" + + + def test_bad_dna(self): + bad_dna = 'ATGCU' + with pytest.raises(ValueError): + DNASequence(bad_dna) + + + def test_transcribe_type(self): + dna = 'ATCG' + result = DNASequence(dna).transcribe() + target = RNASequence('AUCG') + assert type(target) == type(result) + + + + def test_not_rna(self): + dna = 'ATGC' + with pytest.raises(ValueError): + RNASequence(dna) + + +class TestAminoAcidSequences(): + def test_aminoacids_alphabet(self): + peptide = AminoAcidSequence('') + target = set('ARNDCEQGHILKMFPSTWYVUOarndceqghilkmfpstwyvuo') + result = AminoAcidSequence.get_alphabet(peptide) + assert target == result + + + def test_zero_peptid_weight(self): + peptide = AminoAcidSequence('') + result = peptide.molecular_weight() + target = 0 + assert target == result, 'expected 0 for empty string' + + + def test_peptid_weight_calculation(self): + peptide = AminoAcidSequence('VALINE') + result = peptide.molecular_weight() + target = 639.81 + assert target == result, 'expected 639.81 for model peptide VALINE' + + def test_error_message_for_bad_seq(self): + bad_peptide = 'XAXAXA' + with pytest.raises(ValueError) as exinfo: + AminoAcidSequence(bad_peptide) + assert exinfo.value.message == "Unexpected symbol 'X' in the sequence" diff --git a/tool_for_bioinformatics.py b/tool_for_bioinformatics.py deleted file mode 100644 index f2acf42..0000000 --- a/tool_for_bioinformatics.py +++ /dev/null @@ -1,190 +0,0 @@ -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='')