如果第一次閱讀本系列文檔請(qǐng)先移步閱讀【ROSALIND】【練Python，學(xué)生信】00 寫在前面 ?謝謝配合~

題目：

計(jì)算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。

?

測(cè)試數(shù)據(jù)

>Rosalind_6404

CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC

TCCCACTAATAATTCTGAGG

>Rosalind_5959

CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT

ATATCCATTTGTCAGCAGACACGC

>Rosalind_0808

CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC

TGGGAACCTGCGGGCAGTAGGTGGAAT

測(cè)試輸出

Rosalind_0808

60.919540

?

背景

在分析語言時(shí)，首先要做的往往是計(jì)算各個(gè)字母出現(xiàn)的頻率。這種思想可以推廣到生物學(xué)的序列分析中。在雙鏈DNA中，由于堿基互補(bǔ)配對(duì)，G與C出現(xiàn)的頻率相同，A與T出現(xiàn)的頻率相同。同一物種內(nèi)的不同個(gè)體GC含量基本相同，物種間則存在差異，尤其是真核和原核生物之間差異較大，真核生物GC含量約為50%，原核生物則顯著高于50%。因此面對(duì)一段信息未知的DNA序列，可以利用GC含量初步判斷其種屬。

FASTA格式是一種基于文本用于表示核酸序列或多肽序列的格式，其第一行是由大于號(hào)“>”（較常用）或分號(hào)“;”打頭的任意文字說明，用于序列標(biāo)記；從第二行開始為序列本身，只允許使用既定的核苷酸或氨基酸編碼符號(hào)，通常核苷酸符號(hào)大小寫均可，而氨基酸常用大寫字母。一般每行60～80個(gè)字母。FASTA格式已成為生物信息學(xué)領(lǐng)域的一項(xiàng)標(biāo)準(zhǔn)。

?

思路

本題出現(xiàn)了fasta格式文件，為了后續(xù)操作方便，做好的方法是定義一個(gè)專門讀入fasta文件的函數(shù)，可以直接輸出兩個(gè)列表，分別包含序列說明和序列本身，且索引相同。隨后就可以用字符串本身帶的方法實(shí)現(xiàn)GC含量的計(jì)算。

?

Python知識(shí)點(diǎn)

Python使用def關(guān)鍵字定義函數(shù)，將常用的操作封裝成函數(shù)可以實(shí)現(xiàn)代碼的復(fù)用，還可以將所用的工具函數(shù)寫在專門的文件里，使用時(shí)用import關(guān)鍵字導(dǎo)入，為閱讀方便，本系列不采用這種方法，所有的用到的函數(shù)都臨時(shí)定義。

?

代碼

def readfasta(lines):
??? """讀入并處理FASTA文件的函數(shù)"""

??? 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", "")? # 把分行的序列拼接成一個(gè)字符串

??????????? numlines += 1

??????? if numlines == len(lines):

??????????? seq.append(seqplast.replace("\n", ""))

??? seq = seq[1:]

return index, seq

?

?

def countGC(list):
??? """計(jì)算GC含量的函數(shù)"""

??? 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))