Problem
The GC-content of a DNA string is given by the percentage of symbols in the string that are 'C' or 'G'. For example, the GC-content of "AGCTATAG" is 37.5%. Note that the reverse complement of any DNA string has the same GC-content.
DNA strings must be labeled when they are consolidated into a database. A commonly used method of string labeling is called FASTA format. In this format, the string is introduced by a line that begins with '>', followed by some labeling information. Subsequent lines contain the string itself; the first line to begin with '>' indicates the label of the next string.
In Rosalind's implementation, a string in FASTA format will be labeled by the ID "Rosalind_xxxx", where "xxxx" denotes a four-digit code between 0000 and 9999.
Given: At most 10 DNA strings in FASTA format (of length at most 1 kbp each).
Return: The ID of the string having the highest GC-content, followed by the GC-content of that string. Rosalind allows for a default error of 0.001 in all decimal answers unless otherwise stated; please see the note on absolute error below.
Sample Dataset
>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT
Sample Output
Rosalind_0808
60.919540 方法一:
# -*- coding: utf-8 -*- # to open FASTA format sequence file:
s=open('Computing_GC_Content.txt','r').readlines() # to create two lists, one for names, one for sequences
name_list=[]
seq_list=[] data='' # to put the sequence from several lines together for line in s:
line=line.strip()
for i in line:
if i == '>':
name_list.append(line[1:])
if data:
seq_list.append(data) #将每一行的的核苷酸字符串连接起来
data='' # 合完后data 清零
break
else:
line=line.upper()
if all([k==k.upper() for k in line]): #验证是不是所有的都是大写
data=data+line
seq_list.append(data) # is there a way to include the last sequence in the for loop?
GC_list=[]
for seq in seq_list:
i=0
for k in seq:
if k=="G" or k=='C':
i+=1
GC_cont=float(i)/len(seq)*100.0
GC_list.append(GC_cont) m=max(GC_list)
print name_list[GC_list.index(m)] # to find the index of max GC
print "{:0.6f}".format(m) # 保留6位小数
方法二:
# -*- coding: utf-8 -*- def parse_fasta(s):
results = {}
strings = s.strip().split('>')
# Python split()通过指定分隔符对字符串进行切片,如果参数num 有指定值,则仅分隔 num 个子字符串 for s in strings:
if len(s) == 0:
continue
# 如果字符串长度为0,就跳出循环。 parts = s.split()
label = parts[0]
bases = ''.join(parts[1:]) results[label] = bases return results def gc_content(s):
n = len(s)
m = 0 for c in s:
if c == 'G' or c == 'C':
m += 1 return 100 * (float(m) / n) if __name__ == "__main__": small_dataset = """
>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT
""" #large_dataset = open('datasets/rosalind_gc.txt').read() results = parse_fasta(small_dataset)
results = dict([(k, gc_content(v)) for k, v in results.iteritems()])
# 这里iteritem()和item()功能是一样的
# 前一个results输出,名称+序列,后一个results输出,名称+百分比 highest_k = None
highest_v = 0 for k, v in results.iteritems():
if v > highest_v:
highest_k = k
highest_v = v
# 输出GC含量高的
print highest_k
print '%f%%' % highest_v
方法三:
# -*- coding: utf-8 -*- ### 5. Computing GC Content ###
from operator import itemgetter
from collections import OrderedDict seqTest = OrderedDict()
gcContent = OrderedDict() with open('Computing_GC_Content.txt', 'rt') as f:
for line in f:
line = line.rstrip()
if line.startswith('>'):
seqName = line[1:]
seqTest[seqName] = ''
continue
seqTest[seqName] += line.upper() for ke, val in seqTest.items():
totalLength = len(val)
gcNumber = val.count('G') + val.count('C')
gcContent[ke] = (float(gcNumber) / totalLength)*100 sortedGCContent = sorted(gcContent.items(), key=itemgetter(1))
largeName = sortedGCContent[-1][0]
largeGCContent = sortedGCContent[-1][1] print ('most GC ratio gene is %s and it is %s ' % (largeName, largeGCContent))