【ROSALIND】【练Python，学生信】05 计算DNA序列GC含量
如果第一次阅读本系列文档请先移步阅读【ROSALIND】【练Python，学生信】00 写在前面 谢谢配合~
题目：
计算DNA序列GC含量（Computing GC content）
Given: At most 10 DNA strings in FASTA format (of length at most 1 kbp each).
所给：不超过10条DNA序列，每条最少1kbp，以FASTA格式提供。
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.
需得：GC含量最高序列的ID，以及其GC含量。允许误差0.001。
测试数据
>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT
测试输出
Rosalind_0808
60.919540
背景
在分析语言时，首先要做的往往是计算各个字母出现的频率。这种思想可以推广到生物学的序列分析中。在双链DNA中，由于碱基互补配对，G与C出现的频率相同，A与T出现的频率相同。同一物种内的不同个体GC含量基本相同，物种间则存在差异，尤其是真核和原核生物之间差异较大，真核生物GC含量约为50%，原核生物则显著高于50%。因此面对一段信息未知的DNA序列，可以利用GC含量初步判断其种属。
FASTA格式是一种基于文本用于表示核酸序列或多肽序列的格式，其第一行是由大于号“>”（较常用）或分号“;”打头的任意文字说明，用于序列标记；从第二行开始为序列本身，只允许使用既定的核苷酸或氨基酸编码符号，通常核苷酸符号大小写均可，而氨基酸常用大写字母。一般每行60～80个字母。FASTA格式已成为生物信息学领域的一项标准。
思路
本题出现了fasta格式文件，为了后续操作方便，做好的方法是定义一个专门读入fasta文件的函数，可以直接输出两个列表，分别包含序列说明和序列本身，且索引相同。随后就可以用字符串本身带的方法实现GC含量的计算。
Python知识点
Python使用def关键字定义函数，将常用的操作封装成函数可以实现代码的复用，还可以将所用的工具函数写在专门的文件里，使用时用import关键字导入，为阅读方便，本系列不采用这种方法，所有的用到的函数都临时定义。
代码
def readfasta(lines):
"""读入并处理FASTA文件的函数"""
seq = []
index = []
seqplast = ""
numlines = 0
for i in lines:
if '>' in i: # 判断是序列行还是说明行
index.append(i.replace("\n", "").replace(">", ""))
seq.append(seqplast.replace("\n", ""))
seqplast = ""
numlines += 1
else:
seqplast = seqplast + i.replace("\n", "") # 把分行的序列拼接成一个字符串
numlines += 1
if numlines == len(lines):
seq.append(seqplast.replace("\n", ""))
seq = seq[1:]
return index, seq
def countGC(list):
"""计算GC含量的函数"""
numGC = list.count('G') + list.count('C')
perGC = float(numGC) / len(list)
return perGC * 100
f = open('rosalind_gc.txt', 'r')
lines = f.readlines()
f.close()
(index, seq) = readfasta(lines) # 接收序列名和序列
result = []
for i in seq:
result.append(countGC(i))
maxID = index[result.index(max(result))].replace('>', "").replace("\n", "")
seqGC = max(result)
print(maxID)
print(round(seqGC,6))
