2
\$\begingroup\$

I have updated my my gene sequencing program from my previous post. That post explains what each functions accomplish. If you need clarifications feel free to ask.

Any tips to make the code more concise or efficient are much appreciated. Are there any functions that you think would be useful?

Here are example FASTA and GenBank files, for reference.

import random
import re
from collections import defaultdict

nucleotides = ("A", "C", "G", "T")
rnanucleotides = ("A", "C", "G", "U")
rev_compliment = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G'}
RNA_codon_table = {"UUU": "F", "CUU": "L", "AUU": "I", "GUU": "V",
    "UUC": "F", "CUC": "L", "AUC": "I", "GUC": "V",
    "UUA": "L", "CUA": "L", "AUA": "I", "GUA": "V",
    "UUG": "L", "CUG": "L", "AUG": "M", "UCU": "S",
    "GUG": "V", "CCU": "P", "ACU": "T", "GCU": "A",
    "UCC": "S", "CCC": "P", "ACC": "T", "GCC": "A",
    "UCA": "S", "CCA": "P", "ACA": "T", "GCA": "A",
    "UCG": "S", "CCG": "P", "ACG": "T", "GCG": "A",
    "UAU": "Y", "CAU": "H", "AAU": "N", "GAU": "D",
    "UAC": "Y", "CAC": "H", "AAC": "N", "GAC": "D",
    "UAA": "_", "CAA": "Q", "AAA": "K", "GAA": "E",
    "UAG": "_", "CAG": "Q", "AAG": "K", "GAG": "E",
    "UGU": "C", "CGU": "R", "AGU": "S", "GGU": "G",
    "UGC": "C", "CGC": "R", "AGC": "S", "GGC": "G",
    "UGA": "_", "CGA": "R", "AGA": "R", "GGA": "G",
    "UGG": "W", "CGG": "R", "AGG": "R", "GGG": "G"
    }
amino_acid_weights = {
"A": 71.03711,
"C": 103.00919,
"D": 115.02694,
"E": 129.04259,
"F": 147.06841,
"G": 57.02146,
"H": 137.05891,
"I": 113.08406,
"K": 128.09496,
"L": 113.08406,
"M": 131.04049,
"N": 114.04293,
"P": 97.05276,
"Q": 128.05858,
"R": 156.10111,
"S": 87.03203,
"T": 101.04768,
"V": 99.06841,
"W": 186.07931,
"Y": 163.06333,
"_": 0,
}

def title_screen():
    print("""-. .-.   .-. .-.   .-. .-.   .  
||\|||\ /|||\|||\ /|||\|||\ /|
|/ \|||\|||/ \|||\|||/ \|||\||
~   `-~ `-`   `-~ `-`   `-~ `-""")
    print("DNA SEQUENCE ANALYZER by Ethan Hetrick\n")

def user_selection():
    response = input("What would you like to do?\n"
                     "1. Input your own DNA sequence.\n"
                     "2. Input your own RNA sequence.\n"
                     "3. Generate random sequence data.\n"
                     "4. Import a file in FASTA format.\n"
                     "5. Import a file in GenBank format.\n")
    if response == '1':
        seq = input("Input your DNA sequence here: \n").upper().strip()
        seq = validate_seq(seq)
        tag = "Your DNA sequence"
        return seq, tag
    elif response == '3':
        try:
            x = int(input("How many nucleotides long do you want your sequence? Enter an integer.\n"))
            seq = ''.join([random.choice(nucleotides) for nuc in range(x)])
            tag = "Random DNA sequence"
            return seq, tag
        except ValueError:
            print("Invalid response.\n")
            run('seq', 'tag')
    elif response == '2':
        rna = input("Input your RNA sequence here: \n").upper().strip()
        rna = validate_seq(rna)
        seq = rev_transcription(rna)
        tag = "Your RNA sequence"
        return seq, tag
    elif response == '4':
        try:
            for i in fasta():
                tag = i[0]
                seq = i[1]
                if "U" in seq:
                    seq = rev_transcription(seq)
                return seq, tag
        except FileNotFoundError:
            print("\nFile not found. Please input a valid file path.\n")
    elif response == '5':
        try:
            seq, tag = genbank_to_fasta()
            return seq, tag
        except FileNotFoundError:
            print("\nFile not found. Please input a valid file path.\n")
    else:
        print("Invalid response.\n")

def fasta():
    sequences = defaultdict(str)
    file = input(r'Input the path to your file: ')
    with open(f'{file}') as f:
        lines = f.readlines()
    current_tag = None
    list = []
    for line in lines:
        m = re.match('^>(.+)', line)
        if m:
            current_tag = m.group(1)
        else:
            sequences[current_tag] += line.strip()
    for tag, seq in sequences.items():
        list.append([tag, seq])
    return list

def genbank_to_fasta():
    file = input(r'Input the path to your file: ')
    with open(f'{file}') as f:
        gb = f.readlines()
        locus = re.search('NC_\d+\.\d+', gb[3]).group()
        region = re.search('(\d+)?\.+(\d+)', gb[2])
        definition = re.search('\w.+', gb[1][10:]).group()
        definition = definition.replace(definition[-1], "")
        tag = locus + ":" + region.group(1) + "-" + region.group(2) + " " + definition
        sequence = ""
        for line in (gb):
            pattern = re.compile('[atgc]{10}')
            matches = pattern.finditer(line)
            for match in matches:
                sequence += match.group().upper()
        end_pattern = re.search('[atgc]{1,9}', gb[-3])
        sequence += end_pattern.group().upper()
        return sequence, tag

def validate_seq(dnaseq):
    if len(dnaseq) == 0:
        print("Empty sequence.\n")
        user_selection()
    else:
        for nuc in dnaseq:
            if nuc not in nucleotides or rnanucleotides:
                print("Invalid sequence.\n")
                user_selection()
            else:
                return dnaseq

def nuc_count(dnaseq, tag):
    print("\n================================================")
    print(f'THE ANALYSIS: >{tag}\n')
    if len(dnaseq) >= 1000:
        print(f'Total: {len(dnaseq)/1000} kbp\n')
    else:
        print(f'Total: {len(dnaseq)} bp\n')
    print("Nucleotide frequency:")
    for letter in nucleotides:
        letter_total = dnaseq.count(letter)
        letter_per = (letter_total / len(dnaseq))
        print(f'{letter}: {letter_total} : {letter_per:%}')
    GC = ((dnaseq.count("G") + dnaseq.count("C"))/(len(dnaseq)))
    print(f'GC content: {GC:%}')

def rev_transcription(rnaseq):
    dnaseq = rnaseq.replace("U", "T")
    return dnaseq

def DNA_to_cDNA(dnaseq):
    cdna = "".join([rev_compliment[nuc] for nuc in dnaseq])
    return cdna

def transcription(dnaseq):
    rna = dnaseq.replace("T", "U")
    print(f"RNA:  5' {rna} 3'")
    return rna

def translation(dnaseq, init_pos=0):
    dnaseq = dnaseq.replace("T", "U")
    return [RNA_codon_table[dnaseq[pos:pos + 3]] for pos in range(init_pos, len(dnaseq) - 2, 3)]

def gen_reading_frames(dnaseq):
    frames = []
    frames.append(translation(dnaseq, 0))
    frames.append(translation(dnaseq, 1))
    frames.append(translation(dnaseq, 2))
    frames.append(translation(DNA_to_cDNA(dnaseq)[::-1], 0))
    frames.append(translation(DNA_to_cDNA(dnaseq)[::-1], 1))
    frames.append(translation(DNA_to_cDNA(dnaseq)[::-1], 2))
    return frames

def prot_from_rf(aa_seq):
    prot1 = []
    proteins = []
    for aa in aa_seq:
        if aa == "_":
            proteins.extend(prot1)
            prot1 = []
        else:
            if aa == "M":
                prot1.append("")
            for i in range(len(prot1)):
                prot1[i] += aa
    return proteins

def all_proteins_from_rfs(dnaseq, startReadPos=0, endReadPos=0):
    if endReadPos > startReadPos:
        rfs = gen_reading_frames(dnaseq[startReadPos: endReadPos])
    else:
        rfs = gen_reading_frames(dnaseq)
    all_proteins = []
    for rf in rfs:
        prots = prot_from_rf(rf)
        for p in prots:
            all_proteins.append(p)
    return all_proteins

def protein_weight(protein):
    weights = (([amino_acid_weights[protein[pos: pos + 1]] for pos in range(0, len(protein))]))
    weight = round(sum(weights), 3)
    return weight

def printing(seq, frames):
    print("\nThe 6 possible reading frames:")
    for i, frame in enumerate(frames):
        print(f'{i + 1}. {"".join(frame)}')
    print("\nAll possible proteins:")
    list = []
    for prot in all_proteins_from_rfs(seq):
        if prot not in list:
            list.append(prot)
    list.sort(key=len, reverse=True)
    for i, prot in enumerate(list):
        print(f'{i + 1}. {prot}: {protein_weight(prot)} Da')
    print("================================================\n")

def run(seq, tag):
    if seq == 'seq' and tag == 'tag':
        seq, tag = user_selection()
    nuc_count(seq, tag)
    cseq = DNA_to_cDNA(seq)
    print(f"\nDNA:  5' {seq} 3'")
    print(f"cDNA: 3' {cseq} 5'")
    transcription(seq)
    translation(seq)
    all_proteins_from_rfs(seq, startReadPos=0, endReadPos=0)
    frames = gen_reading_frames(seq)
    printing(seq, frames)

title_screen()
run('seq', 'tag')

Example output

================================================
THE ANALYSIS: >NC_000006.12:26156329-26157115 Homo sapiens chromosome 6, GRCh38.p13 Primary Assembly

Total: 787 bp

Nucleotide frequency:
A: 221 : 28.081321%
C: 246 : 31.257942%
G: 230 : 29.224905%
T: 90 : 11.435832%
GC content: 60.482846%

DNA:  5' CGAGTCCCGGCCAGTGCCTCTGCTTCCGGCTCGAATTGCTCTCGCTCACGCTTGCCTTCAACATGTCCGAGACTGCGCCTGCCGCGCCCGCTGCTCCGGCCCCTGCCGAGAAGACTCCCGTGAAGAAGAAGGCCCGCAAGTCTGCAGGTGCGGCCAAGCGCAAAGCGTCTGGGCCCCCGGTGTCCGAGCTCATTACTAAAGCTGTTGCCGCCTCCAAGGAGCGCAGCGGCGTATCTTTGGCCGCTCTCAAGAAAGCGCTGGCAGCCGCTGGCTATGACGTGGAGAAGAACAACAGCCGCATCAAGCTGGGTCTCAAGAGCCTGGTGAGCAAGGGCACCCTGGTGCAGACCAAGGGCACCGGCGCGTCGGGTTCCTTCAAACTCAACAAGAAGGCGGCCTCTGGGGAAGCCAAGCCTAAGGCTAAAAAGGCAGGCGCGGCCAAGGCCAAGAAGCCAGCAGGAGCGGCGAAGAAGCCCAAGAAGGCGACGGGGGCGGCCACCCCCAAGAAGAGCGCCAAGAAGACCCCAAAGAAGGCGAAGAAGCCGGCTGCAGCTGCTGGAGCCAAAAAAGCGAAAAGCCCGAAAAAGGCGAAAGCAGCCAAGCCAAAAAAGGCGCCCAAGAGCCCAGCGAAGGCCAAAGCAGTTAAACCCAAGGCGGCTAAACCAAAGACCGCCAAGCCCAAGGCAGCCAAGCCAAAGAAGGCGGCAGCCAAGAAAAAGTAGAAAGTTCCTTTGGCCAACTGCTTAGAAGCCCAACACAACCCAAAGGCTCTTTTCAGAGCCACCCA 3'
cDNA: 3' GCTCAGGGCCGGTCACGGAGACGAAGGCCGAGCTTAACGAGAGCGAGTGCGAACGGAAGTTGTACAGGCTCTGACGCGGACGGCGCGGGCGACGAGGCCGGGGACGGCTCTTCTGAGGGCACTTCTTCTTCCGGGCGTTCAGACGTCCACGCCGGTTCGCGTTTCGCAGACCCGGGGGCCACAGGCTCGAGTAATGATTTCGACAACGGCGGAGGTTCCTCGCGTCGCCGCATAGAAACCGGCGAGAGTTCTTTCGCGACCGTCGGCGACCGATACTGCACCTCTTCTTGTTGTCGGCGTAGTTCGACCCAGAGTTCTCGGACCACTCGTTCCCGTGGGACCACGTCTGGTTCCCGTGGCCGCGCAGCCCAAGGAAGTTTGAGTTGTTCTTCCGCCGGAGACCCCTTCGGTTCGGATTCCGATTTTTCCGTCCGCGCCGGTTCCGGTTCTTCGGTCGTCCTCGCCGCTTCTTCGGGTTCTTCCGCTGCCCCCGCCGGTGGGGGTTCTTCTCGCGGTTCTTCTGGGGTTTCTTCCGCTTCTTCGGCCGACGTCGACGACCTCGGTTTTTTCGCTTTTCGGGCTTTTTCCGCTTTCGTCGGTTCGGTTTTTTCCGCGGGTTCTCGGGTCGCTTCCGGTTTCGTCAATTTGGGTTCCGCCGATTTGGTTTCTGGCGGTTCGGGTTCCGTCGGTTCGGTTTCTTCCGCCGTCGGTTCTTTTTCATCTTTCAAGGAAACCGGTTGACGAATCTTCGGGTTGTGTTGGGTTTCCGAGAAAAGTCTCGGTGGGT 5'
RNA:  5' CGAGUCCCGGCCAGUGCCUCUGCUUCCGGCUCGAAUUGCUCUCGCUCACGCUUGCCUUCAACAUGUCCGAGACUGCGCCUGCCGCGCCCGCUGCUCCGGCCCCUGCCGAGAAGACUCCCGUGAAGAAGAAGGCCCGCAAGUCUGCAGGUGCGGCCAAGCGCAAAGCGUCUGGGCCCCCGGUGUCCGAGCUCAUUACUAAAGCUGUUGCCGCCUCCAAGGAGCGCAGCGGCGUAUCUUUGGCCGCUCUCAAGAAAGCGCUGGCAGCCGCUGGCUAUGACGUGGAGAAGAACAACAGCCGCAUCAAGCUGGGUCUCAAGAGCCUGGUGAGCAAGGGCACCCUGGUGCAGACCAAGGGCACCGGCGCGUCGGGUUCCUUCAAACUCAACAAGAAGGCGGCCUCUGGGGAAGCCAAGCCUAAGGCUAAAAAGGCAGGCGCGGCCAAGGCCAAGAAGCCAGCAGGAGCGGCGAAGAAGCCCAAGAAGGCGACGGGGGCGGCCACCCCCAAGAAGAGCGCCAAGAAGACCCCAAAGAAGGCGAAGAAGCCGGCUGCAGCUGCUGGAGCCAAAAAAGCGAAAAGCCCGAAAAAGGCGAAAGCAGCCAAGCCAAAAAAGGCGCCCAAGAGCCCAGCGAAGGCCAAAGCAGUUAAACCCAAGGCGGCUAAACCAAAGACCGCCAAGCCCAAGGCAGCCAAGCCAAAGAAGGCGGCAGCCAAGAAAAAGUAGAAAGUUCCUUUGGCCAACUGCUUAGAAGCCCAACACAACCCAAAGGCUCUUUUCAGAGCCACCCA 3'

The 6 possible reading frames:
1. RVPASASASGSNCSRSRLPSTCPRLRLPRPLLRPLPRRLP_RRRPASLQVRPSAKRLGPRCPSSLLKLLPPPRSAAAYLWPLSRKRWQPLAMTWRRTTAASSWVSRAW_ARAPWCRPRAPARRVPSNSTRRRPLGKPSLRLKRQARPRPRSQQERRRSPRRRRGRPPPRRAPRRPQRRRRSRLQLLEPKKRKARKRRKQPSQKRRPRAQRRPKQLNPRRLNQRPPSPRQPSQRRRQPRKSRKFLWPTA_KPNTTQRLFSEPP
2. ESRPVPLLPARIALAHACLQHVRDCACRARCSGPCREDSREEEGPQVCRCGQAQSVWAPGVRAHY_SCCRLQGAQRRIFGRSQESAGSRWL_RGEEQQPHQAGSQEPGEQGHPGADQGHRRVGFLQTQQEGGLWGSQA_G_KGRRGQGQEASRSGEEAQEGDGGGHPQEERQEDPKEGEEAGCSCWSQKSEKPEKGESSQAKKGAQEPSEGQSS_TQGG_TKDRQAQGSQAKEGGSQEKVESSFGQLLRSPTQPKGSFQSHP
3. SPGQCLCFRLELLSLTLAFNMSETAPAAPAAPAPAEKTPVKKKARKSAGAAKRKASGPPVSELITKAVAASKERSGVSLAALKKALAAAGYDVEKNNSRIKLGLKSLVSKGTLVQTKGTGASGSFKLNKKAASGEAKPKAKKAGAAKAKKPAGAAKKPKKATGAATPKKSAKKTPKKAKKPAAAAGAKKAKSPKKAKAAKPKKAPKSPAKAKAVKPKAAKPKTAKPKAAKPKKAAAKKK_KVPLANCLEAQHNPKALFRAT
4. WVALKRAFGLCWASKQLAKGTFYFFLAAAFFGLAALGLAVFGLAALGLTALAFAGLLGAFFGLAAFAFFGLFAFLAPAAAAGFFAFFGVFLALFLGVAAPVAFLGFFAAPAGFLALAAPAFLALGLASPEAAFLLSLKEPDAPVPLVCTRVPLLTRLLRPSLMRLLFFSTS_PAAASAFLRAAKDTPLRSLEAATALVMSSDTGGPDALRLAAPADLRAFFFTGVFSAGAGAAGAAGAVSDMLKASVSESNSSRKQRHWPGL
5. GWL_KEPLGCVGLLSSWPKELSTFSWLPPSLAWLPWAWRSLV_PPWV_LLWPSLGSWAPFLAWLLSPFSGFSLFWLQQLQPASSPSLGSSWRSSWGWPPPSPSWASSPLLLASWPWPRLPF_P_AWLPQRPPSC_V_RNPTRRCPWSAPGCPCSPGS_DPA_CGCCSSPRHSQRLPALS_ERPKIRRCAPWRRQQL___ARTPGAQTLCAWPHLQTCGPSSSRESSRQGPEQRARQAQSRTC_RQA_ARAIRAGSRGTGRDS
6. GGSEKSLWVVLGF_AVGQRNFLLFLGCRLLWLGCLGLGGLWFSRLGFNCFGLRWALGRLFWLGCFRLFRAFRFFGSSSCSRLLRLLWGLLGALLGGGRPRRLLGLLRRSCWLLGLGRACLFSLRLGFPRGRLLVEFEGTRRAGALGLHQGALAHQALETQLDAAVVLLHVIASGCQRFLESGQRYAAALLGGGNSFSNELGHRGPRRFALGRTCRLAGLLLHGSLLGRGRSSGRGRRSLGHVEGKREREQFEPEAEALAGT

All possible proteins:
1. MSETAPAAPAAPAPAEKTPVKKKARKSAGAAKRKASGPPVSELITKAVAASKERSGVSLAALKKALAAAGYDVEKNNSRIKLGLKSLVSKGTLVQTKGTGASGSFKLNKKAASGEAKPKAKKAGAAKAKKPAGAAKKPKKATGAATPKKSAKKTPKKAKKPAAAAGAKKAKSPKKAKAAKPKKAPKSPAKAKAVKPKAAKPKTAKPKAAKPKKAAAKKK: 21833.993 Da
2. MTWRRTTAASSWVSRAW: 2034.001 Da
3. MRLLFFSTS: 1082.558 Da
================================================
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.