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
================================================
