Looking for feedback on my DNA analyzer program that reads FASTA and generates cDNA, RNA, reading frames, and codons

Question

import random

# Variables, lists, tuples, and dictionaries.
nucleotides = ("A", "C", "G", "T")
rev_compliment = {'A': 'T', 'T': 'A', 'G': 'C', 'C': 'G'}
seq = ""
description = ""
frames = []
proteins = []
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,
}

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

# Lets user choose a custom sequence, generate a random one, or import a FASTA file.
def randomvschoice(retry = 1):
    global seq
    global description
    while retry:
        response = input("What would you like to do?\n"
                         "1. Input your own DNA sequence.\n"
                         "2. Generate a random DNA sequence.\n"
                         "3. Import a FASTA file.\n")
        if response == '1':
            seq = input("Paste your sequence here: \n").upper().replace(' ', '').replace('\n', '')
            return seq
        elif response == '2':
            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)])
                return seq
            except ValueError:
                print("\nInvalid response.\n")
                retry += 1
        elif response == '3':
            try:
                seq, description = open_fasta()
                return seq, description
            except FileNotFoundError:
                print("\nFile not found. Please input a valid file path.\n")
                retry += 1
        else:
            print("You must answer 1, 2, or 3.\n")
            retry += 1

# Converts the FASTA file into a readable sequence.
def open_fasta():
    path = input(r'Input path to the fasta file: ')
    with open(path, 'r') as file:
        fasta = file.readlines()
        seq = fasta[1:]
        description = fasta[0]
    seq = "".join(seq).replace("\n", "")
    return seq, description

# Validates the sequence to make sure it only contains A, C, G, or T.
def validate_seq(dnaseq):
    try:
        1/len(dnaseq)
        for nuc in dnaseq:
            if nuc not in nucleotides:
                print("Invalid sequence. Only A, C, T, and G are accepted characters.\n")
                randomvschoice()
    except ZeroDivisionError:
        print("Invalid sequence. Only A, C, T, and G are accepted characters.\n")
        randomvschoice()

# Counts nucleotides (A, T, G, and C) and gives percentages + total nucleotides.
def nuc_count(dnaseq):
    print("\n================================================")
    print(f'THE ANALYSIS: {description}\n')
    print(f'Total: {len(dnaseq)} nucleotides\n')
    print("Nucleotide frequency:")
    for letter in nucleotides:
        letter_total = dnaseq.count(letter)
        letter_per = round(letter_total / len(dnaseq) * 100, 1)
        print(letter + ": " + str(letter_total) + "  :  " + str(letter_per) + "%")

# Converts DNA to cDNA by using the rev_compliment dictionary.
def DNA_to_cDNA1(dnaseq):
    dnaseq = "".join([rev_compliment[nuc] for nuc in dnaseq])[::1]
    print("\nDNA:  " + "5' " + seq + " 3'")
    # This is to show the base pairing for shorter sequences as if you run it in a .exe, over 150 nucleotides
    # and it looks like a mess if the sequence is more than one line.
    if len(seq) <= 150:
        print("         " + "|" * len(dnaseq))
    print("cDNA: " + "3' " + dnaseq + " 5'")
    return ''.join([rev_compliment[nuc] for nuc in dnaseq])[::-1]

# I would like both of these to be one function since this code is redundant. Was not sure how to only access the
# return and not the print statement of DNA_to_cDNA1.
def DNA_to_cDNA(dnaseq):
    dnaseq = "".join([rev_compliment[nuc] for nuc in dnaseq])[::1]
    return ''.join([rev_compliment[nuc] for nuc in dnaseq])[::-1]

# Converts cDNA to RNA by replacing T with U.
def transcription(dnaseq):
    dnaseq = dnaseq.replace("T", "U")
    print("RNA:  5' " + str(dnaseq) + " 3'")
    return dnaseq

# GC content calculator. Calculates percentages of G+C in the sequence.
def GC_content(dnaseq):
    content = round(((dnaseq.count("C") + dnaseq.count("G")) / len(dnaseq)) * 100, 3)
    print("GC content: " + str(content) + "%")

# Converts mRNA to protein using the RNA_codon_table.
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)]

# Generates all 6 reading frames from nucleotide positions 0, 1, and 2 from both the DNA and cDNA sequences.
def gen_reading_frames(dnaseq):
    global frames
    frames.append(translation(dnaseq, 0))
    frames.append(translation(dnaseq, 1))
    frames.append(translation(dnaseq, 2))
    frames.append(translation(DNA_to_cDNA(dnaseq), 0))
    frames.append(translation(DNA_to_cDNA(dnaseq), 1))
    frames.append(translation(DNA_to_cDNA(dnaseq), 2))
    return frames

# Generates a protein from the reading frames.
def prot_from_rf(aa_seq):
    prot1 = []
    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

# This code I took inspiration from a video but seems like too much. It runs the gen_reading_frames function
# and generates all proteins from the 6 reading frames generated. The first statement is to make sure the
# dnaseq is not 0.
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

# Protein weight calculator. Uses the amino_acid_weights dictionary to calculate the mass of each protein in Daltons.
def protein_weight(protein):
    for aa in protein:
        weights = (([amino_acid_weights[protein[pos: pos + 1]] for pos in range(0, len(protein))]))
        weight = round(sum(weights), 3)
        return weight

# This was all the extra code I had to write to make the program present the data in a user-friendly fashion.
# This is the main function I would like to do away with or change.
def stupid():
    print("\nThe 6 possible reading frames:")
    x = 0
    # Numbers and prints the reading frames in a string.
    for frame in frames:
        x += 1
        print(f'{x}. {"".join(frame)}')
    print("\nAll possible proteins:")
    list = []
    # Even more code to get the proteins to print numbered and sorted by length.
    for prot in all_proteins_from_rfs(seq):
        if prot not in list:
            list.append(prot)
    list.sort(key=len, reverse=True)
    y = 0
    for prot in list:
        y += 1
        print(f'{y}. {prot}: {protein_weight(prot)} Da')
    print("================================================\n")

# This is how I run the program. Also probably not optimal.
def run():
    randomvschoice()
    validate_seq(seq)
    nuc_count(seq)
    DNA_to_cDNA1(seq)
    transcription(seq)
    translation(seq)
    GC_content(seq)
    all_proteins_from_rfs(seq, startReadPos=0, endReadPos=0)
    stupid()
    run()

title_screen()
run()

I am a biologist that decided to learn coding and apply it to my field. This is my first project and I would love some feedback. I noticed by my last post that I have some bad coding habits. I was told the mutable global variables was not the most appropriate in this context and I was able to remove some. I have a lot of redundant code and readability issues. Even when I go back through it is hard for me to understand exactly what I wrote. I would like to condense it and clarify it as much as possible as I am planning to add more functions to it and possibly turning it into an application in the future. I simply hate the current structure of my code and it makes it hard for me to keep adding to it in its current state. Any advice will be much appreciated!

Answer 1

Data domain restriction

Your nucleotides should be represented as an Enum, given that there are only four possible values. The values of the enum can be strings, but usage in the rest of the program should be through the enum to get more confidence in correctness.

External data

Your codon table should be externalized, to a CSV or JSON file.

Globals

the mutable global variables was not the most appropriate in this context

Indeed. Constants (nucleotides, rev_compliment, the codon table and the amino acids) are fine as global constants, but the others are problematic. Pull them out of global scope and pass them around as parameters and return values for your various methods.

Docstrings

# Validates the sequence to make sure it only contains A, C, G, or T.
def validate_seq(dnaseq):

should actually look more like

def validate_seq(dnaseq):
    """
    Validates the sequence to make sure it only contains A, C, G, or T.
    """

for a few reasons, including that Python can reflect on that string.

Length checking

This:

1/len(dnaseq)
...
except ZeroDivisionError:

is a nasty way to do a length check. Your error message string is also wrong in that case, because it isn't a matter of invalid characters - the sequence is empty. Instead, spell it out:

if len(dnaseq) == 0:
    raise ValueError('Empty DNA sequence')

Further, attempt to keep console printing and scanning at the top, so that methods like this can use exceptions, which should consistently be the internal error-signalling mechanism of the application.

String interpolation

    print(letter + ": " + str(letter_total) + "  :  " + str(letter_per) + "%")

can be

    print(f'{letter}: {letter_total}  :  {letter_per:%}')

Note that the use of % in this way obviates a multiplication by 100.

String formation

This:

    dnaseq = "".join([rev_compliment[nuc] for nuc in dnaseq])[::1]

has a few issues:

• Don't make an inner [] list; pass the generator directly to join
• The step size of 1 for a slice is the default, so you can omit it and write [:]
• But why do a slice at all? Just omit the slice.

Console lines

it looks like a mess if the sequence is more than one line.

Sure, but 150 will be far above the width of some consoles and far below the width of others. Getting the exact width is platform-dependent and can be tricky. Most console applications simply assume that the user is an adult and will be able to pipe through a scrolling frame buffer if they want, and do not do anything special with long lines.

Counting

This:

dnaseq.count("C") + dnaseq.count("G")

is not as efficient as making a Counter instance. It's a built-in and will do the counting for both symbols in one pass.

Answer 2

The answer from @Reinderien covers a lot of things well. I want to particularly pick up on the Global Variables aspect as I think it's the single part that will help you to improve your code the most.

You seem to have down the idea that functions take an input and that you can get an output from them using return, but you don't seem to use that output in the run function. Giving the intermediate variables name would also make it easier to work out what stage you were at. e.g.

def run(): 
    ...
    c_dna_seq = DNA_to_cDNA(seq)
    rna_seq = transcription(c_dna_seq)
    ...

Another thing that I would recommend for readability is using verbs for function names rather than nouns. e.g. convert_dna_to_cdna() or transcribe(seq)