diff --git a/README.md b/README.md index c176f0a..253535e 100644 --- a/README.md +++ b/README.md @@ -1 +1,177 @@ # HW5. Modules + +> *This is repository for the fifth homework assignment of the BI Python 2023 course* +## Table of contents + + * [Project description](#project-description) + * [Main part](#Main-part) + * [Usage](#Usage) + * [Parametrs](#Parameters) + * [Examples](#examples) + * [Requirements](#Requirements) + * [Contact](#contact) + * [Educational result](#Учебный-результат) + * [See also](#see-also) + +## Project Description +This project includes several projects, and the first thought that came to my mind when viewing the task was ~~это starдец~~ it would be difficult. By that time I had neither 3 nor 4 homework assignments. But I didn't have time to cry. To reach the goal, a person needs only one thing – to go. Slowly I did 3 homework, then I got to 4 and now I have reached the writing of the README file of the fifth homework + +![go to goal](https://st3.depositphotos.com/1001201/12919/i/450/depositphotos_129195498-stock-photo-businessman-running-fast.jpg) + +## Main part +This project consists of 2 points; it combines 3 mini-projects from previous tasks [HERE](https://github.com/1799mrd/HW4_Functions2/tree/HW4_Mukhametshina), [HERE](https://github.com/Python-BI-2023/hw3-functions-1799mrd/tree/develop) and [HERE](https://github.com/1799mrd/Bioinformatic_tools/blob/develop/all_functions/fastq_filtration_tools.py). +In this work need to get a real package with bioinformatics utilities.
+ +There are 2 points in this project: + +### Writing a utility for working with fastq sequences. +The main function takes 4 arguments as input: `seqs`, `gc_bounds`, `length_bounds`, `quality_threshold`: + +- `seqs` - a dictionary consisting of fastq sequences. The structure is as follows. Key - string, sequence name. The value is a tuple of two strings: sequence and quality. Since I haven’t read files yet, we use a python dictionary. +The `example_data.py` script contains an example for debugging. +- `gc_bounds` - GC composition interval (in percentage) for filtering (by default it is equal to (0, 100), i.e. all reads are saved). If you pass one number as an argument, it is considered to be the upper limit. Examples: `gc_bounds = (20, 80)` - we save only reads with GC content from 20 to 80%, `gc_bounds = 44.4` - we save reads with GC content less than 44.4%. +- `length_bounds` - length interval for filtering, everything is similar to `gc_bounds`, but by default it is (0, 2**32). +- `quality_threshold` - the threshold value of the average quality of the read for filtering, by default it is 0 (phred33 scale). Reads with average quality across all nucleotides below the threshold are discarded.

+ + +**Based on the results of the work** the main function returns a similar dictionary, consisting only of those sequences that pass all the conditions. All described intervals include both upper and lower boundaries.
+ + +### Assembling projects for working with fastq sequences, amino acid sequences, DNA and RNA sequences. +HERE HERE AND HERE into a single package for a novice bioinformatician + + ```python + -/ + |- bioinf_script.py # (импорты и 3 функции) + |- README.md + |- all_functions/ + |- amino_acid_tools.py # (ДНК-функции) + |- dna_rna_tools.py # (белковые функции) + |- fastq_filtration_tools.py # (FASTQ-функции) + ``` + +# Bioinformatic tool - DNA/RNA sequence processing/amino acid sequence processing/fastq sequences filtering +This tool is designed to work with fastq sequences, amino acid sequences, DNA and RNA sequences. + +## Usage +There three functions: `filter_fastq_sequences`, `amino_acid_tools`, `run_dna_rna_tools`. + +If you are working with fastq sequences, run `filter_fastq_sequences` function
+It takes as input a dictionary containing an unlimited number of fastq sequences, filtering bounds on gc composition, length, and average read quality. The function then filters all sequences in the dictionary. Sequences can contain both uppercase and lowercase letters. The result of the function is a dictionary with filtered sequences. + +If you are working with amino acid sequences, run `amino_acid_tools` function
+It takes as input an arbitrary number of arguments with amino acid sequences (str) as well as the name of the procedure to be performed. + +If you are working with dna or rna sequences, run `run_dna_rna_tools` function
+It takes as input an arbitrary number of arguments with DNA or RNA sequences (str) and the name of the procedure to be performed. + +For `amino_acid_tools` and `run_dna_rna_tools` input sequences can contain both uppercase and lowercase letters, but the last argument with the function name must match the possible functions. 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. + + +## Parameters + +#### For `filter_fastq_sequences` function + +Function works with fastq file:
+##### Parameters:
`seqs`, `gc_bounds`, `length_bounds`, `quality_threshold`:
+- **seqs** is a dictionary consisting of fastq sequences. The structure is as follows. The key is a string, the name of the sequence. The value is a tuple of two strings: sequence and quality. In fact, this is the content of a fastq file, but you and I have not yet been through reading files. So for now we are using a python dictionary. Then it will be enough to add reading files and writing them to a dictionary of this kind so that everything works from beginning to end. In the script `example_data.py ` an example has been made for you to debug.
+- **gc_bounds** - GC composition interval (in percent) for filtering (by default is (0, 100), i.e. all reads are saved). If a single number is passed to the argument, it is assumed that this is the upper bound. Examples: `gc_bounds = (20, 80)` - save only reads with GC composition from 20 to 80%, `gc_bounds = 44.4` - save reads with GC composition less than 44.4%.
+- **length_bounds** - the length interval for filtering, everything is similar to `gc_bounds`, but by default it is `(0, 2**32)`.
+- **quality_threshold** is the threshold value of the average quality of the read for filtering, by default it is 0 (phred33 scale). Reads with average quality for all nucleotides below the threshold are discarded.

+##### Return:
+`similar dictionary` consisting only of those sequences that have passed all the conditions. All the described intervals include both upper and lower bounds. + +#### For `amino_acid_tools` function +The following options for amino acid sequence processing are available at the moment:
+##### Parametrs: +- `molecular weight`: calculates the molecular weight of an amino acid chain +- `long_amino_code`: converts translations from one letter to translations from three letters +- `amino_to_rna` : translates a sequence of amino acids into a possible sequence of nucleic acids +- `amino_seq_charge`: estimates the total charge of the amino acid chain in a neutral aqueous solution (pH = 7) +##### Return: +If one sequence is supplied, a string with the result is returned.
+If several are submitted, a list of strings is returned.
+Depending on the function performed, the following returns will occur: +- `long_amino_code` (str) or (list): translated sequence from one-letter in three-letter code +- `molecular_weight` (int) or (list): amino acid sequence molecular weight number or list of numbers +- `amino_to_rna` (str) or (list): possible RNA sequence +- `amino_seq_charge` (str) or (list): "positive", "negative" or "neutral" + +#### For `run_dna_rna_tools` function +- **transcribe**: print transcribed sequence +- **reverse**: print reversed sequence +- **complement**: print complementary sequence +- **reverse_complement**: print reverse complementary sequence + +## Examples + +### For `filter_fastq_sequences` function +Below is an example of processing a fastq sequences. +```python +EXAMPLE_FASTQ = { + '@SRX079804:1:SRR292678:1:1101:21885:21885': ('ACAGCAACATAAACATGATGGGATGGCGTAAGCCCCCGAGATATCAGTTTACCCAGGATAAGAGATTAAATTATGAGCAACATTATTAA', 'FGGGFGGGFGGGFGDFGCEBB@CCDFDDFFFFBFFGFGEFDFFFF;D@DD>C@DDGGGDFGDGG?GFGFEGFGGEF@FDGGGFGFBGGD'), + '@SRX079804:1:SRR292678:1:1101:24563:24563': ('ATTAGCGAGGAGGAGTGCTGAGAAGATGTCGCCTACGCCGTTGAAATTCCCTTCAATCAGGGGGTACTGGAGGATACGAGTTTGTGTG', 'BFFFFFFFB@B@A<@D>BDDACDDDEBEDEFFFBFFFEFFDFFF=CC@DDFD8FFFFFFF8/+.2,@7<<:?B/:<><-><@.A*C>D'), + '@SRX079804:1:SRR292678:1:1101:30161:30161': ('GAACGACAGCAGCTCCTGCATAACCGCGTCCTTCTTCTTTAGCGTTGTGCAAAGCATGTTTTGTATTACGGGCATCTCGAGCGAATC', 'DFFFEGDGGGGFGGEDCCDCEFFFFCCCCCB>CEBFGFBGGG?DE=:6@=>AD?D8DCEE:>EEABE5D@5:DDCA;EEE-DCD'), + '@SRX079804:1:SRR292678:1:1101:47176:47176': ('TGAAGCGTCGATAGAAGTTAGCAAACCCGCGGAACTTCCGTACATCAGACACATTCCGGGGGGTGGGCCAATCCATGATGCCTTTG', 'FF@FFBEEEEFFEFFD@EDEFFB=DFEEFFFE8FFE8EEDBFDFEEBE+EFBEEB=EFA>EEBD=DAEEEEB9)99>B99BC)@,@<9CDD=C,5;B::?@;A'), + '@SRX079804:1:SRR292678:1:1101:170868:170868': ('CTGCCGAGACTGTTCTCAGACATGGAAAGCTCGATTCGCATACACTCGCTGAGTAAGAGAGTCACACCAAATCACAGATT', 'E;FFFEGFGIGGFBG;C6D<@C7CDGFEFGFHDFEHHHBBHHFDFEFBAEEEEDE@A2=DA:??C3:@>EEBEEHEFEHHFFHH?FGBGFBBD77B;;C?FFFFGGFED.BBABBG@DBBE'), + '@SRX079804:1:SRR292678:1:1101:190845:190845': ('CCTCAGCGTGGATTGCCGCTCATGCAGGAGCAGATAATCCCTTCGCCATCCCATTAAGCGCCGTTGTCGGTATTCC', 'FF@FFCFEECEBEC@@BBBBDFBBFFDFFEFFEB8FFFFFFFFEFCEB/>BBA@AFFFEEEEECE;ACD@DBBEEE'), + '@SRX079804:1:SRR292678:1:1101:198993:198993': ('AGTTATTTATGCATCATTCTCATGTATGAGCCAACAAGATAGTACAAGTTTTATTGCTATGAGTTCAGTACAACA', '<<<=;@B??@<>@><48876EADEG6B.BB@.?+98204<:<>@?A=@EFEFFFEEFB'), + '@SRX079804:1:SRR292678:1:1101:204480:204480': ('AGTGAGACACCCCTGAACATTCCTAGTAAGACATCTTTGAATATTACTAGTTAGCCACACTTTAAAATGACCCG', '<98;<@@@:@CD@BCCDD=DBBCEBBAAA@9???@BCDBCGF=GEGDFGDBEEEEEFFFF=EDEE=DCD@@BBC') + } +filter_fastq = filter_fastq_sequences(EXAMPLE_FASTQ, (20, 80), (0, 80), 30) +print(filter_fastq) + +#{'@SRX079804:1:SRR292678:1:1101:170868:170868': ('CTGCCGAGACTGTTCTCAGACATGGAAAGCTCGATTCGCATACACTCGCTGAGTAAGAGAGTCACACCAAATCACAGATT', #'E;FFFEGFGIGGFBG;C6D<@C7CDGFEFGFHDFEHHHBBHHFDFEFBAEEEEDE@A2=DA:??C3:@>EEBEEHEFEHHFFHH?FGBGFBBD77B;;C?FFFFGGFED.BBABBG@DBBE'), +#'@SRX079804:1:SRR292678:1:1101:190845:190845': ('CCTCAGCGTGGATTGCCGCTCATGCAGGAGCAGATAATCCCTTCGCCATCCCATTAAGCGCCGTTGTCGGTATTCC', #'FF@FFCFEECEBEC@@BBBBDFBBFFDFFEFFEB8FFFFFFFFEFCEB/>BBA@AFFFEEEEECE;ACD@DBBEEE'), +#'@SRX079804:1:SRR292678:1:1101:204480:204480': ('AGTGAGACACCCCTGAACATTCCTAGTAAGACATCTTTGAATATTACTAGTTAGCCACACTTTAAAATGACCCG','<98;<@@@:@CD@BCCDD=DBBCEBBAAA@9???@BCDBCGF=GEGDFGDBEEEEEFFFF=EDEE=DCD@@BBC')} +``` + +### For `amino_acid_tools` function +Below is an example of processing an amino acid sequence. + +```python +Using the function for translated sequence from one-letter in three-letter code: +amino_acid_tools('PLfHnfPdD','long_amino_code') # 'ProLeuPheHisAsnPheProAspAsp' +amino_acid_tools('random', 'CpUPQWhmrY','molecular_weight') # [('рандомная последовательнсть', 'FySiDfGym', 1124.43), 1367.39] +amino_acid_tools('cvdrLepaW', 'amino_to_rna') # 'ugugucgaccggCUAgagccggcgUGG' +amino_acid_tools('cvdrLepaW', 'VurgdOhio', 'amino_seq_charge') # ['negativ', 'positiv'] +``` + +#### For `run_dna_rna_tools` function +Below is an example of processing DNA and RNA sequences. + +```python +run_dna_rna_tools('ATG', 'transcribe') # 'AUG' +run_dna_rna_tools('ATG', 'reverse') # 'GTA' +run_dna_rna_tools('AtG', 'complement') # 'TaC' +run_dna_rna_tools('ATg', 'reverse_complement') # 'cAT' +run_dna_rna_tools('ATG', 'aT', 'reverse') # ['GTA', 'Ta'] +``` + +## Requirements +You need to import modules: +```python +from all_functions.amino_acid_tools import amino_acid_tools +from all_functions.dna_rna_tools import run_dna_rna_tools +from all_functions.fastq_filtration_tools import filter_fastq_sequences +``` +## Учебный-результат +> * В данном проекте я научилась пользоваться pyCharm в котором довольно удобно писать Readme file. А также работу с модулями +> * К сожалению, я не успела получить комментарии к 4 дз и не смогла исправить код к данному проекту, но перечитать все написанное мной пришлось. +> * Вспомнила fastq + +## Contacts +- [Mukhametshina Regina] 1709mrd@gmail.com + +![health](https://u.9111s.ru/uploads/202209/12/27ec4ad913eaa7cf6277658d30dbb3cf.jpg) + +## See also + +Про fastq-файлы можно посмотреть, например, [здесь](https://stepik.org/lesson/32398/step/1?unit=12379), про определение качества на основе ASCII-кода можно почитать [здесь](https://support.illumina.com/help/BaseSpace_OLH_009008/Content/Source/Informatics/BS/QualityScoreEncoding_swBS.htm), также вам может пригодиться таблица [кодировки ASCII](https://www.asciitable.com/) и функция ord (документация [en](https://docs.python.org/3/library/functions.html#ord), [ru](https://docs-python.ru/tutorial/vstroennye-funktsii-interpretatora-python/funktsija-ord/)). diff --git a/all_functions/amino_acid_tools.py b/all_functions/amino_acid_tools.py new file mode 100644 index 0000000..d1f5e00 --- /dev/null +++ b/all_functions/amino_acid_tools.py @@ -0,0 +1,195 @@ +amino_acid = 'ARNDCEQGHILKMFPSTWYVUOarndceqghilkmfpstwyvuo' +short_code = list(amino_acid) +long_code = ['Ala', 'Arg', 'Asn', 'Asp', 'Cys', 'Glu', 'Gln', 'Gly', 'His', 'Ile', 'Leu', 'Lys', 'Met', 'Phe', 'Pro', + 'Ser', 'Thr', 'Trp', 'Tyr', 'Val', 'Sec', 'Pyl', + 'Ala', 'Arg', 'Asn', 'Asp', 'Cys', 'Glu', 'Gln', 'Gly', 'His', 'Ile', 'Leu', 'Lys', 'Met', 'Phe', 'Pro', + 'Ser', 'Thr', 'Trp', 'Tyr', 'Val', 'Sec', 'Pyl'] +codon_table = { + 'A': ['GCU', 'GCC', 'GCA', 'GCG'], + 'R': ['CGU', 'CGC', 'CGA', 'CGG', 'AGA', 'AGG'], + 'N': ['AAU', 'AAC'], + 'D': ['GAU', 'GAC'], + 'C': ['UGU', 'UGC'], + 'Q': ['CAA', 'CAG'], + 'E': ['GAA', 'GAG'], + 'G': ['GGU', 'GGC', 'GGA', 'GGG'], + 'H': ['CAU', 'CAC'], + 'I': ['AUU', 'AUC', 'AUA'], + 'L': ['UUA', 'UUG', 'CUU', 'CUC', 'CUA', 'CUG'], + 'K': ['AAA', 'AAG'], + 'M': ['AUG'], + 'F': ['UUU', 'UUC'], + 'P': ['CCU', 'CCC', 'CCA', 'CCG'], + 'S': ['UCU', 'UCC', 'UCA', 'UCG', 'AGU', 'AGC'], + 'T': ['ACU', 'ACC', 'ACA', 'ACG'], + 'W': ['UGG'], + 'Y': ['UAU', 'UAC'], + 'V': ['GUU', 'GUC', 'GUA', 'GUG'], + 'STOP': ['UAA', 'UAG', 'UGA'], + 'f': ['uuu', 'uuc'], + 'l': ['uua', 'uug', 'cuu', 'cuc', 'cua', 'cug'], + 's': ['ucu', 'ucc', 'uca', 'ucg', 'agu', 'agc'], + 'y': ['uau', 'uac'], + 'c': ['ugu', 'ugc'], + 'w': ['ugg'], + 'p': ['ccu', 'ccc', 'cca', 'ccg'], + 'h': ['cau', 'cac'], + 'q': ['caa', 'cag'], + 'r': ['cgu', 'cgc', 'cga', 'cgg', 'aga', 'agg'], + 'i': ['auu', 'auc', 'aua'], + 'm': ['aug'], + 't': ['acu', 'acc', 'aca', 'acg'], + 'n': ['aau', 'aac'], + 'k': ['aaa', 'aag'], + 'v': ['guu', 'guc', 'gua', 'gug'], + 'a': ['gcu', 'gcc', 'gca', 'gcg'], + 'd': ['gau', 'gac'], + 'e': ['gaa', 'gag'], + 'g': ['ggu', 'ggc', 'gga', 'ggg'], + 'stop': ['uaa', 'uag', 'uga'] +} +weight_amino = [71.08, 156.2, 114.1, 115.1, 103.1, 129.1, 128.1, 57.05, 137.1, 113.2, 113.2, 128.2, 131.2, 147.2, 97.12, 87.08, + 101.1, 186.2, 163.2, 99.13, 168.05, 255.3, + 71.08, 156.2, 114.1, 115.1, 103.1, 129.1, 128.1, 57.05, 137.1, 113.2, 113.2, 128.2, 131.2, 147.2, 97.12, 87.08, + 101.1, 186.2, 163.2, 99.13, 168.05, 255.3] + +import random + +def long_amino_code(sequence): + """ + Function translates a given sequence of one-letter amino acids + into a more understandable sequence of amino acids consisting of three letters + + Parameters: + sequence (str): each letter refers to one-letter coded proteinogenic amino acids or "random" + Returns: + (str) translated in three-letter code + """ + if sequence != 'random': + d_names = dict(zip(short_code, long_code)) + recording = sequence.maketrans(d_names) + return sequence.translate(recording) + else: + len = int(input("введите желаемую длину: ")) + bases = list(amino_acid) + amino_sequencqe = ''.join(random.choice(bases) for i in range(len)) + d_names = dict(zip(short_code, long_code)) + recording = amino_sequencqe.maketrans(d_names) + return "рандомная последовательнсть", amino_sequencqe, amino_sequencqe.translate(recording) + +def molecular_weight(sequence): + """ + Function calculates molecular weight of the amino acid chain + Parameters: + sequence (str): each letter refers to one-letter coded proteinogenic amino acids or "random" + Returns: + weight (float) Molecular weight of tge given amino acid chain in Da + """ + if sequence != 'random': + molecular_weights = dict(zip(short_code, weight_amino)) + weight = sum(molecular_weights.get(aa, 0) for aa in sequence) + return weight + else: + len = int(input("введите желаемую длину: ")) + bases = list(amino_acid) + amino_sequencqe = ''.join(random.choice(bases) for i in range(len)) + molecular_weights = dict(zip(short_code, weight_amino)) + weight = sum(molecular_weights.get(aa, 0) for aa in amino_sequencqe) + return "рандомная последовательнсть", amino_sequencqe, weight + +def amino_to_rna(amino_sequence): + """ + Function translates an amino acid sequence into a possible RNA sequence + Parameters: + amino_sequence (str) or "random" + Returns: + (str) possible RNA sequence + """ + if amino_sequence != 'random': + rna_sequence = "" + + for aminoacid in amino_sequence: + if aminoacid in codon_table: + codons = codon_table[aminoacid] + # Selecting one random codon + codon = random.choice(codons) + rna_sequence += codon + else: + print("Unknown amino acid code: ", aminoacid) + + return rna_sequence + else: + len = int(input("введите желаемую длину: ")) + bases = list(amino_acid) + amino_sequencqe = ''.join(random.choice(bases) for i in range(len)) + rna_sequence = "" + + for aminoacid in amino_sequencqe: + if aminoacid in codon_table: + codons = codon_table[aminoacid] + # Selecting one random codon + codon = random.choice(codons) + rna_sequence += codon + else: + print("Unknown amino acid code: ", aminoacid) + return "рандомная последовательнсть", amino_sequencqe, rna_sequence + + +def amino_seq_charge(amino_sequence): + """ + Function evaluates the overall charge of the aminoacid chain in neutral aqueous solution (pH = 7) + Parameters: + amino_sequence (str): amino acid sequence of proteinogenic amino acids or "random" + Returns: + (str): "positive", "negative" or "neutral" + """ + if amino_sequence != 'random': + aminoacid_charge = {'R': 1, 'D': -1, 'E': -1, 'K': 1, 'O': 1} + charge = 0 + for aminoacid in amino_sequence.upper(): + if aminoacid in 'RDEKO': + charge += aminoacid_charge[aminoacid] + if charge > 0: + return 'positiv' + elif charge < 0: + return 'negativ' + else: + return 'neutral' + else: + len = int(input("введите желаемую длину: ")) + bases = list(amino_acid) + amino_sequencqe = ''.join(random.choice(bases) for i in range(len)) + aminoacid_charge = {'R': 1, 'D': -1, 'E': -1, 'K': 1, 'O': 1} + charge = 0 + for aminoacid in amino_sequencqe.upper(): + if aminoacid in 'RDEKO': + charge += aminoacid_charge[aminoacid] + if charge > 0: + return "рандомная последовательнсть", amino_sequencqe, 'positiv' + elif charge < 0: + return "рандомная последовательнсть", amino_sequencqe, 'negativ' + else: + return "рандомная последовательнсть", amino_sequencqe, 'neutral' + +def amino_seqs(amino_sequence): + """ + Leaves only the amino acid sequences from the fed into the function. + Parameters: + amino_sequence (list): amino acid sequence list or "random" + Returns: + amino_seqs (list): amino acid sequence list without non amino acid sequence + """ + if amino_sequence != 'random': + aminoac_seqs = [] + for seq in amino_sequence: + unique_chars = set(seq) + amino_acids = set(amino_acid) + if unique_chars <= amino_acids: + aminoac_seqs.append(seq) + return aminoac_seqs + else: + len = int(input("введите желаемую длину: ")) + bases = list(amino_acid) + amino_sequencqe = ''.join(random.choice(bases) for i in range(len)) + aminoac_seqs = list(amino_sequencqe) + return "рандомная последовательнсть", amino_sequencqe, aminoac_seqs diff --git a/all_functions/dna_rna_tools.py b/all_functions/dna_rna_tools.py new file mode 100644 index 0000000..27e925e --- /dev/null +++ b/all_functions/dna_rna_tools.py @@ -0,0 +1,39 @@ +def transcribe(seq): + transcribed_seq = "" + for base in seq: + if base == 'T': + transcribed_seq += 'U' + elif base == 't': + transcribed_seq += 'u' + else: + transcribed_seq += base + return transcribed_seq + +def reverse(seq): + return seq[::-1] + +def complement(seq): + complement_seq = "" + for base in seq: + if base == 'A': + complement_seq += 'T' + elif base == 'T': + complement_seq += 'A' + elif base == 'C': + complement_seq += 'G' + elif base == 'G': + complement_seq += 'C' + elif base == 'a': + complement_seq += 't' + elif base == 't': + complement_seq += 'a' + elif base == 'c': + complement_seq += 'g' + elif base == 'g': + complement_seq += 'c' + else: + complement_seq += base + return complement_seq + +def reverse_complement(seq): + return complement(reverse(seq)) diff --git a/all_functions/fastq_filtration_tools.py b/all_functions/fastq_filtration_tools.py new file mode 100644 index 0000000..bbaf04a --- /dev/null +++ b/all_functions/fastq_filtration_tools.py @@ -0,0 +1,46 @@ +EXAMPLE_FASTQ = { + '@SRX079804:1:SRR292678:1:1101:21885:21885': ('ACAGCAACATAAACATGATGGGATGGCGTAAGCCCCCGAGATATCAGTTTACCCAGGATAAGAGATTAAATTATGAGCAACATTATTAA', 'FGGGFGGGFGGGFGDFGCEBB@CCDFDDFFFFBFFGFGEFDFFFF;D@DD>C@DDGGGDFGDGG?GFGFEGFGGEF@FDGGGFGFBGGD'), + '@SRX079804:1:SRR292678:1:1101:24563:24563': ('ATTAGCGAGGAGGAGTGCTGAGAAGATGTCGCCTACGCCGTTGAAATTCCCTTCAATCAGGGGGTACTGGAGGATACGAGTTTGTGTG', 'BFFFFFFFB@B@A<@D>BDDACDDDEBEDEFFFBFFFEFFDFFF=CC@DDFD8FFFFFFF8/+.2,@7<<:?B/:<><-><@.A*C>D'), + '@SRX079804:1:SRR292678:1:1101:30161:30161': ('GAACGACAGCAGCTCCTGCATAACCGCGTCCTTCTTCTTTAGCGTTGTGCAAAGCATGTTTTGTATTACGGGCATCTCGAGCGAATC', 'DFFFEGDGGGGFGGEDCCDCEFFFFCCCCCB>CEBFGFBGGG?DE=:6@=>AD?D8DCEE:>EEABE5D@5:DDCA;EEE-DCD'), + '@SRX079804:1:SRR292678:1:1101:47176:47176': ('TGAAGCGTCGATAGAAGTTAGCAAACCCGCGGAACTTCCGTACATCAGACACATTCCGGGGGGTGGGCCAATCCATGATGCCTTTG', 'FF@FFBEEEEFFEFFD@EDEFFB=DFEEFFFE8FFE8EEDBFDFEEBE+EFBEEB=EFA>EEBD=DAEEEEB9)99>B99BC)@,@<9CDD=C,5;B::?@;A'), + '@SRX079804:1:SRR292678:1:1101:170868:170868': ('CTGCCGAGACTGTTCTCAGACATGGAAAGCTCGATTCGCATACACTCGCTGAGTAAGAGAGTCACACCAAATCACAGATT', 'E;FFFEGFGIGGFBG;C6D<@C7CDGFEFGFHDFEHHHBBHHFDFEFBAEEEEDE@A2=DA:??C3:@>EEBEEHEFEHHFFHH?FGBGFBBD77B;;C?FFFFGGFED.BBABBG@DBBE'), + '@SRX079804:1:SRR292678:1:1101:190845:190845': ('CCTCAGCGTGGATTGCCGCTCATGCAGGAGCAGATAATCCCTTCGCCATCCCATTAAGCGCCGTTGTCGGTATTCC', 'FF@FFCFEECEBEC@@BBBBDFBBFFDFFEFFEB8FFFFFFFFEFCEB/>BBA@AFFFEEEEECE;ACD@DBBEEE'), + '@SRX079804:1:SRR292678:1:1101:198993:198993': ('AGTTATTTATGCATCATTCTCATGTATGAGCCAACAAGATAGTACAAGTTTTATTGCTATGAGTTCAGTACAACA', '<<<=;@B??@<>@><48876EADEG6B.BB@.?+98204<:<>@?A=@EFEFFFEEFB'), + '@SRX079804:1:SRR292678:1:1101:204480:204480': ('AGTGAGACACCCCTGAACATTCCTAGTAAGACATCTTTGAATATTACTAGTTAGCCACACTTTAAAATGACCCG', '<98;<@@@:@CD@BCCDD=DBBCEBBAAA@9???@BCDBCGF=GEGDFGDBEEEEEFFFF=EDEE=DCD@@BBC') + } + +def filter_fastq_sequences(seqs, gc_bounds=(0, 100), length_bounds=(0, 2**32), quality_threshold=0): + filtered_seqs = {} + + for name, (sequence, quality) in seqs.items(): + gc_content = calculate_gc_content(sequence) + seq_length = len(sequence) + avg_quality = calculate_avg_quality(quality) + + if is_within_bounds(gc_content, gc_bounds) and is_within_bounds(seq_length, length_bounds) and avg_quality >= quality_threshold: + filtered_seqs[name] = (sequence, quality) + + return filtered_seqs + +def calculate_gc_content(sequence): + gc_count = sequence.count('G') + sequence.count('C') + total_count = len(sequence) + + return (gc_count / total_count) * 100 + +def calculate_avg_quality(quality): + total_quality = sum(ord(q)-33 for q in quality) + seq_length = len(quality) + + return total_quality / seq_length + +def is_within_bounds(value, bounds): + if type(bounds) in (int, float): + return value < bounds + else: + lower_bound, upper_bound = bounds + return lower_bound <= value <= upper_bound diff --git a/bioinf_script.py b/bioinf_script.py new file mode 100644 index 0000000..0975709 --- /dev/null +++ b/bioinf_script.py @@ -0,0 +1,105 @@ +from all_functions.amino_acid_tools import amino_acid_tools +from all_functions.dna_rna_tools import run_dna_rna_tools +from all_functions.fastq_filtration_tools import filter_fastq_sequences + +def filter_fastq_sequences(seqs, gc_bounds=(0, 100), length_bounds=(0, 2**32), quality_threshold=0): + filtered_seqs = {} + + for seq_name, (sequence, quality) in seqs.items(): + # Фильтрация по GC составу + gc_content = calculate_gc_content(sequence) + if isinstance(gc_bounds, tuple): + if not gc_bounds[0] <= gc_content <= gc_bounds[1]: + continue + elif isinstance(gc_bounds, (int, float)): + if gc_content > gc_bounds: + continue + + # Фильтрация по длине + sequence_length = len(sequence) + if isinstance(length_bounds, tuple): + if not length_bounds[0] <= sequence_length <= length_bounds[1]: + continue + elif isinstance(length_bounds, (int, float)): + if sequence_length > length_bounds: + continue + + # Фильтрация по качеству + average_quality = calculate_average_quality(quality) + if average_quality < quality_threshold: + continue + + # Добавление отфильтрованной последовательности в результат + filtered_seqs[seq_name] = (sequence, quality) + + return filtered_seqs + +def run_dna_rna_tools(*args): + sequences = args[:-1] + procedure = args[-1] + + results = [] + for seq in sequences: + if procedure == 'transcribe': + result = transcribe(seq) + elif procedure == 'reverse': + result = reverse(seq) + elif procedure == 'complement': + result = complement(seq) + elif procedure == 'reverse_complement': + result = reverse_complement(seq) + else: + result = "Unknown procedure" + + results.append(result) + + if len(results) == 1: + return results[0] + else: + return results + +def amino_acid_tools(*args: str): + """ + Performs functions for working with amino acid sequences. + + Parameters: + The function should accept an unlimited number of protein sequences (str) as input, + the last variable should be the function (str) that you want to execute. + The amino acid sequence can consist of both uppercase and lowercase letters. + Input example: + amino_acid_tools('LVElkPL','CpUPQWhmrY','McgMmLcTTG','molecular_weight') + + or + + amino_acid_tools('LVElkPL','CpUPQWhmrY','random','molecular_weight') + + + Function: + molecular weight: calculates the molecular weight of an amino acid chain + long_amino_code: converts translations from one letter to translations + from three letters + amino_to_rna translates a sequence of amino acids into a possible sequence of nucleic acids + amino_seq_charge: estimates the total charge of the amino acid chain in a neutral aqueous solution (pH = 7) + + Returns: + If one sequence is supplied, a string with the result is returned. + If several are submitted, a list of strings is returned. + Depending on the function performed, the following returns will occur: + long_amino_code (str) or (list): translated sequence from one-letter in three-letter code + molecular_weight (int) or (list): amino acid sequence molecular weight number or list of numbers + amino_to_rna (str) or (list): possible RNA sequence + amino_seq_charge (str) or (list): "positive", "negative" or "neutral" + """ + *seqs, function = args + d_of_functions = {'long_amino_code' : long_amino_code, + 'molecular_weight': molecular_weight, + 'amino_to_rna' : amino_to_rna, + 'amino_seq_charge' : amino_seq_charge} + answer = [] + aminoacid_seqs = amino_seqs(seqs) + for sequence in aminoacid_seqs: + answer.append(d_of_functions[function](sequence)) + if len(answer) == 1: + return answer[0] + else: + return answer