1

所以这个问题基本上给了我 19 个 DNA 序列,并要我制作一个基本的文本表。第一列必须是序列 ID,第二列是序列长度,第三列是“A”的数量,第四列是“G”,第五列是“C”,第六列是“T” ,第7个是%GC,第8个是序列中是否有“TGA”。然后我得到所有这些值并将表写入“dna_stats.txt”

这是我的代码:

fh = open("dna.fasta","r")
Acount = 0
Ccount = 0
Gcount = 0
Tcount = 0
seq=0
alllines = fh.readlines()
for line in alllines:
    if line.startswith(">"):
        seq+=1
        continue
    Acount+=line.count("A")
    Ccount+=line.count("C")
    Gcount+=line.count("G")
    Tcount+=line.count("T")
    genomeSize=Acount+Gcount+Ccount+Tcount
    percentGC=(Gcount+Ccount)*100.00/genomeSize
    print "sequence", seq
    print "Length of Sequence",len(line)
    print Acount,Ccount,Gcount,Tcount
    print "Percent of GC","%.2f"%(percentGC)
    if "TGA" in line:
        print "Yes"
    else:
        print "No"
    fh2 = open("dna_stats.txt","w")
    for line in alllines:
        splitlines = line.split()
        lenstr=str(len(line))
        seqstr = str(seq)
        fh2.write(seqstr+"\t"+lenstr+"\n")

我发现您必须将变量转换为字符串。当我在终端中打印出来时,我已经正确计算了所有值。但是,我一直只得到第一列的 19,而它应该是 1、2、3、4、5 等。来表示所有的序列。我尝试了其他变量,它只得到了整个文件的总量。我开始尝试制作桌子,但还没有完成。

所以我最大的问题是我不知道如何获取每个特定行的变量值。

我是 python 和一般编程的新手,所以任何提示或技巧或任何东西都会真正有帮助。

我正在使用python 2.7版

4

2 回答 2

1

好吧,你最大的问题:

for line in alllines: #1
    ...
    fh2 = open("dna_stats.txt","w")
    for line in alllines: #2
        ....

缩进很重要。这说“对于每一行(#1),打开一个文件,然后再次遍历每一行(#2)......”

取消缩进这些东西。

于 2013-10-17T05:44:41.457 回答
0

这会将信息随时放入字典中,并允许 DNA 序列遍历多行

from __future__ import division # ensure things like 1/2 is 0.5 rather than 0
from collections import defaultdict

fh = open("dna.fasta","r")
alllines = fh.readlines()


fh2 = open("dna_stats.txt","w")
seq=0
data = dict()
for line in alllines:
    if line.startswith(">"):
        seq+=1
        data[seq]=defaultdict(int) #default value will be zero if key is not present hence we can do +=1 without originally initializing to zero        
        data[seq]['seq']=seq
        previous_line_end = "" #TGA might be split accross line
        continue

    data[seq]['Acount']+=line.count("A")
    data[seq]['Ccount']+=line.count("C")
    data[seq]['Gcount']+=line.count("G")
    data[seq]['Tcount']+=line.count("T")
    data[seq]['genomeSize']+=data[seq]['Acount']+data[seq]['Gcount']+data[seq]['Ccount']+data[seq]['Tcount']

    line_over = previous_line_end + line[:3]        
    data[seq]['hasTGA']= data[seq]['hasTGA'] or ("TGA" in line) or (TGA in line_over)
    previous_line_end = str.strip(line[-4:]) #save previous_line_end for next line removing new line character.


for seq in data.keys():    
    data[seq]['percentGC']=(data[seq]['Gcount']+data[seq]['Ccount'])*100.00/data[seq]['genomeSize']    
    s = '%(seq)d, %(genomeSize)d, %(Acount)d, %(Ccount)d, %(Tcount)d, %(Tcount)d, %(percentGC).2f, %(hasTGA)s'
    fh2.write(s % data[seq])

fh.close()    
fh2.close() 
于 2013-10-17T06:44:02.467 回答