使用python计算fasta的GC内容没有错误

Question

人类基因组由 24 条不同的染色体组成（实际上 23 对 = 46 条染色体）。这些染色体被称为1、2、3、...、22、X 和Y。每条染色体都是一个很长的字符串，由'G'、'C'、'A' 和'T' 字符组成（例如，染色体 1 由近 2400 万个字符组成）。我将每条染色体都放在一个文件中（1 号染色体链在1.fa 文件中）。

*.fa 文件称为 fasta 文件，它是 DNA 链信息的标准文件。该文件的结构如下：

>gi|568815591|ref|NC_000007.14| Homo sapiens chromosome 7, GRCh38 Primary Assembly
CATTGCACTCCAGCCTGGGCAAAAACAGCGAAACTCCGTCTCAAAAAAAAAAAAAAGAAAAAAT
TAGCCAGGCATGGTGAAGTTGCAGTGAGCTGAGACTGCACCATTGCACTCCAGCCTGGGTAGCA

如您所见，此类文件的第一行提供了一些有关 GCAT 字符串来源的信息。

我写了这段代码来统计GC内容（G+C字符数占所有字符的比例）：

homo_sapiens_chromosomes_List=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, "x", "y"]
for i in homo_sapiens_chromosomes_List:
    i=str (i)
    file_open= open (i+".fa", "r") #Opening each file

    file_read= file_open.read() #Reading file

    file_read=file_read.upper() #Uppercase characters

    G=float(file_read.count("G"))   #Count G in file

    C=float(file_read.count("C"))   #Count C in file

    A=float(file_read.count("A"))   #Count A in file

    T=float(file_read.count("T"))   #Count A in file

    print "There are %d Gs, %d Cs, %d As and %d Ts, in the DNA strand, of chromosome number %s." % (G, C, A, T, i)

    print "GC content of this chromosome is:", (G+C)*100/(A+T+G+C), "percent"       #Prints GC Content

现在我有一些问题：

1. 如何使这段代码更高效（更快、更短或...）

2. 当我尝试计算 GC 内容时，不属于 DNA 链的 fasta 文件的第一行也会被计算在内。在计算 GC 内容之前，我写了这个函数来删除这一行（这段代码在这行之后：file_read=file_read.upper()）：

代码：

def Fasta_Clean(): #a function to delete the first line of fasta file
    global file_read
    if file_read.isalnum()==False:
        file_read=file_read[1:]
        Fasta_Clean()
    Fasta_Clean()

但是这段代码返回了："RuntimeError: maximum recursion depth exceeded in cmp"，所以我写了这个：

def Fasta_Clean(): #a function to delete the first line of fasta file
    global file_read
    fas=file_read[0:number]
    if fas.isalnum()==False:
        file_read=file_read[1:]
        Fasta_Clean()
    Fasta_Clean()

现在，当变量fas 大于fas=file_read[0:90] 时，我又看到了"RuntimeError: maximum recursion depth exceeded in cmp"。我该如何解决这个问题？

1. 如果 fasta 文件由多个链组成并且具有类似这样的结构（在此示例中，文件由三个不同的链组成）：

例子：

>gi|568815591|ref|NC_000007.14| Homo sapiens chromosome 7, GRCh38 Primary Assembly
GCGAAACTCCGTCTCAAAAAAAAAAAAAAGAAAAAATCATTGCACTCCAGCCTGGGCAAAAACA
CACCATTGCACTCCAGCCTGGGTAGCATAGCCAGGCATGGTGAAGTTGCAGTGAGCTGAGACTG

>gi|568815864|ref|NC_000009.14| Homo sapiens chromosome 8, GRCh38 Primary Assembly
CATTGCACTCCAGCCTGGGCAAAAACAGCGAAACTCCGTCTCAAAAAAAAAAAAAAGAAAAAAT
TAGCCAGGCATGGTGAAGTTGCAGTGAGCTGAGACTGCACCATTGCACTCCAGCCTGGGTAGCA

>gi|568815325|ref|NC_000009.14| Homo sapiens chromosome 9, GRCh38 Primary Assembly
CTGGGCAAAAACAGCGAAACTCCGTCTCAAAAAAAAAAAAAAGAAACATTGCACTCCAGCAAAT
GTGAGCTGAGACTGCACCATTGCTAGCCAGGCATGGTGAAGTTGCAACTCCAGCCTGGGTAGCA

在这种情况下，如何分别计算每条链的 GC 含量？

Answer 1

这应该是您需要的所有代码。

from collections import Counter
chrome_list=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
             12, 13, 14, 15, 16, 17, 18, 19, 20,
             21, 22, "x", "y"]
for i in chrome_list:
    file_ = open('{}.fa'.format(i), 'r')
    broken_file = file_.read().split('\n\n')
    for line in broken_file:
        print Counter(line.split('\n')[1])
    file_.close()

如果您使用的是 windows 或 mac，您可能需要更改 \n\n

Answer 2

对于第 1 点，您可以使用以下方法一次计算多次出现：

from collections import Counter
z = ['G', 'C', 'C', 'A', 'T']
Counter(z)
>>>Counter({'G': 1, 'C': 2, 'A': 1, 'T':1})

第 1 点和第 2 点可以循环吗？：

d = {'G':0, 'C':0, 'A':0, 'T':0}
count = 0
for line in inputfile:
    if count == 0:
        count += 1
        continue
    d[line] += 1