Output when compute GC content is different between using user input and using string

Question

I am learning bioinformatics and i was solving a problem on rosalind.info. Input is a string consist of multiple DNA strand names and the DNA sequences, and Output is the name of the sequence that has the highest GC percentage and the GC percentage. this is the input:

>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT

this is the correct output:

Rosalind_0808
60.919540

This is my code (sorry for this filthy mess):

DNA = input()
wholeRosalind_DNA = DNA.splitlines()
Rosalind = []
Rosalind_DNA = []
for i in wholeRosalind_DNA:
    if "Rosalind" in i:
        Rosalind.append(i)
Countstring = []
maxGCcontent = []
GCcontent = []
def eachRosalind_DNA(DNA):
    for i in range(0,len(Rosalind)):
        if Rosalind[i] in wholeRosalind_DNA:
            a = wholeRosalind_DNA.index(Rosalind[i])
            if a == wholeRosalind_DNA.index(Rosalind[-1]):  
                b = len(wholeRosalind_DNA)+1
            else:
                b = wholeRosalind_DNA.index(Rosalind[i+1])
            Rosalind_DNA.append(wholeRosalind_DNA[a:b])
    for i in Rosalind_DNA:
        Countstring.append(''.join(i))
    for i in Countstring:
        bases = len(i)-14
        CG = i.count("C") + i.count("G")
        GCcontent.append([i[1:14],(100*(CG/bases))])
    for i in GCcontent:
        maxGCcontent.append(i[1])
    result = max(maxGCcontent)
    for i in GCcontent:
        if i[1] == result:
            return (i[0]),(round(i[1],6))
print(eachRosalind_DNA(DNA)[0])
print(eachRosalind_DNA(DNA)[1])

but i get this instead:

  File "GC_Content.py", line 26, in eachRosalind_DNA
    GCcontent.append([i[1:14],(100*(CG/bases))])
ZeroDivisionError: division by zero

i will get the desired output if i just paste the input into the variabe "DNA" without using user input. Please help me fix the code and distinguish those two ways of getting input

Answer 1

Is it cheating to use Biopython(?), will a solution might be helpful anyway. This assumes you have your sequences stored in the seq.fa FASTA file:

from Bio import SeqIO
from Bio.SeqUtils import GC

print(max(((r.id, GC(r.seq)) for r in SeqIO.parse("seq.fa", "fasta")),
          key=lambda x: x[1]))

Output:

('Rosalind_0808', 60.91954022988506)

Answer 2

The error is that bases = 0. Briefly looking at the code carriage returns in your original code will trigger this, as will confusing a faster header with sequence. The latter is because the fasta headers are 14 characters long and subtracting 14 from bases which is of course zero. The former will likely result in 'nan' because it would be 14 subtracted from ''.

I've a simpler solution, download Biopython (using conda or last resort pip), you might want to use SeqIO (Biopython), which would be cooler because it allows you to call Rosalind_0808 directly. Then:

from Bio.SeqUtils import GC
seq = 'CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGACTGGGAACCTGCGGGCAGTAGGTGGAAT'
print (GC(seq))

Done. If you want to splice off the last 14 characters its

GC(seq[:-14])

Answer 3

I think the hard thing in this problem is not the calculation of GC content, it is the processing of DNA strings because they are separated on different lines. I would suggest reading the file line by line to generate complete DNA strings for each Rosalind entry. (I am not able to spot what is causing the error in the code).

seqDic = {}
with open("YouInFile", "r") as inFile:
    for ln in inFile:
        ln = ln.rstrip()
        if ln.startswith(">"):
            ln.name = ln[1:]
            seq = inFile.readline().rstrip() + inFile.readline().rstrip()
            seqDic[ln.name] = seq
        else
            continue

Basically, whenever you encounter a ">" in your file, you write down the Rosalind name, and then you read the next two lines and combine them together. You put the name and the sequence in a dictionary, and iterate over the dictionary to calculate GC content.

Note that I did not debug my code. Let me know if it is not working.