【菜鸟Python练习5】[ROSALIND-GC] 计算GC含量最大的DNA序列

bioinfo.dong

【5】 Computing GC Content
为方便无法翻墙的同学，我把题目也放在帖子里。本题之前有很多人做过了，思路不难，在这里教大家两个小技巧。
* OrderedDict in collections。创建有序的dictionary。两年前我之所以弃perl从Python很重要的一个原因就是Python里有很多实用的包。比如这个。当我们希望dictionary里的元素顺序为元素加入该dictionary的顺序时，就可以创建一个OrderedDict()，这样我们就可以控制元素的输入和输出顺序.
* itemgetter in operator。使用itemgetter可以指定dictionary按key <itemgetter(0) >排序或者按value <itemgetter(1) >排序，从而找出最大GC含量的DNA序列。

[Python] 纯文本查看 复制代码

### 5. Computing GC Content ###
from operator import itemgetter
from collections import OrderedDict
 
seqTest = OrderedDict()
gcContent = OrderedDict()
 
with open('/Users/DONG/Downloads/rosalind_gc.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))

颛顼

程序的18到21行，两个整数做除法，结果很容易变成零，将其中一个数转为浮点数就好了

biostudy

根据楼主的代码，依葫芦画瓢写了一段，也是第一次知道字典原来也能排序。然后也附上自己思路写的一段代码，仅供参考，因为不太熟悉一些python包，所以用嵌套列表解决的排序问题。
楼主代码练习：

[Python] 纯文本查看 复制代码

#!/usr/bin/env python
#_*_ coding: utf-8 _*
from operator import itemgetter
from collections import OrderedDict
#创建一个有序字典集合，根据输入的先后顺序排序
seqTest = OrderedDict()
gcCountent = OrderedDict()
with open('E:\\bioinfo\study\data\\test5.txt','r') as input_file:
    #逐行读取输入文件
    for line in input_file:
        line = line.rstrip()    #删除字符串末尾的指定字符，默认为空格
        #将DNA序列信息保存在有序字典seqTest中
        if line.startswith('>'):
            seqName = line.strip('> ')
        else:
            seqTest[seqName] = line

    for key, value in seqTest.items():
        #生成一个包含GC含量的有序字典
        seq_length = len(value)
        GC_ratio = (float(value.count('C') + value.count('G'))) / seq_length
        gcCountent[key] = GC_ratio
#将字典排序，根据每个元祖中的第二个元素排序
#将字典排序完之后，其形式也不为字典格式
gcCountent_sort = sorted(gcCountent.items(), key=itemgetter(1))
#取最大GC比
large_Name = gcCountent_sort[-1][0]
large_GCRation = gcCountent_sort[-1][1]
print 'GC比最大的DNA序列为：\n%s\n%.6f' % (large_Name, large_GCRation)

我的代码：

[Python] 纯文本查看 复制代码

#!/usr/bin/env python
#_*_ coding: utf-8 _*
gcCount = []
def seq_GCRatio(sequence):
    #输入含有序列信息的字符串，输出该序列中的GC比
    GC_count = float(sequence.count('C') + sequence.count('G'))
    seq_length = len(sequence)
    GC_ratio = GC_count / seq_length * 100
    return GC_ratio

with open('E:\\bioinfo\study\data\\test5.txt', 'r') as input_file:
    for line in input_file:
        seq_list = []
        if line.startswith('>'):
            seq_name = line.strip('[> |  \n]')
        else:
            sequence = line.strip()
            GC_ratio = seq_GCRatio(sequence)
            seq_list.append(seq_name)
            seq_list.append(GC_ratio)
            gcCount.append(seq_list)
    GC_Ratio_sort = sorted(gcCount, key=lambda x:x[1])
print 'GC比最高的序列为:\n%s\n%.6f' % (GC_Ratio_sort[-1][0], GC_Ratio_sort[-1][1])

zckoo007

大神 我用你的代码没跑成功，原来fasta 格式是
>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT

读出来的是
[('Rosalind_6404', 'CCTGCGGAAGATCGGCACTAGAATAGC CAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCCTCCCACTAATAATTCTGAGG'), ('Rosalind_5959', 'CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCTATATCCATTTGTCAGCAGACACGC'), ('Rosalind_0808', '')]
最后 Rosalind_0808,后面的序列读不出来

hxoxh

看代码感觉OrderedDict() 跟dict()区别不大啊，就是最后排序的时候简单了一点点

BACK

请问这里 Orderdict  和dict 有什么区别

liukili

因为直接输出最大值和序列名称，所以没有排序直接用max()
但是前面的代码不够精简，还是楼主的好一些。

[Python] 纯文本查看 复制代码

with open('test.txt') as Dnaseq:
    Dnadic = {}
    while Dnaseq:
        line_1 = Dnaseq.readline().strip()
        line_1 = line_1.strip('>')
        if not (Dnaseq and line_1):
            break
        line_2 = Dnaseq.readline().strip()
        line_3 = Dnaseq.readline().strip()
        key = line_1
        length = float(len(line_2) + len(line_3))
        gcnumber = (line_2.count('G') + line_2.count('C')
                    + line_3.count('G') + line_3.count('C'))
        value = gcnumber / length * 100
        Dnadic[key] = value

print(max(Dnadic, key=Dnadic.get))
print(Dnadic[max(Dnadic, key=Dnadic.get)])

546397641

dna_1 = "CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCCTCCCACTAATAATTCTGAGG"
dna_2 = "CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCTATATCCATTTGTCAGCAGACACGC"
dna_3 = "CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGACTGGGAACCTGCGGGCAGTAGGTGGAAT"
a = b = c = 0
for i in dna_1:
    if i == "C" or i == "G":
        a += 1
GC_1=a/len(dna_1)
for j in dna_2:
    if j == "C" or j == "G":
        b += 1
GC_2=b/len(dna_2)
for k in dna_3:
    if k == "C" or k == "G":
        c += 1
GC_3=c/len(dna_3)
if max(GC_1,GC_2,GC_3) == GC_1:
    print("dna_1:",GC_1)
elif max(GC_1,GC_2,GC_3) == GC_2:
    print("dna_2:",GC_2)
else:
    print("dna_3:",GC_3)

lin_

from Bio import SeqIO
from Bio.SeqUtils import GC
seq_dict = {}
for seq_record in SeqIO.parse(r"D:\bioinformatics\test_seq.txt", "fasta"):
    seq_dict[seq_record.id] = GC(seq_record.seq)
seq_max_id = max(seq_dict, key=seq_dict.get)
print(seq_max_id, seq_dict[seq_max_id])

藏青衣

from Bio import SeqIO
from Bio.SeqUtils import GC


def main():

    gc_values = sorted(GC(rec.seq) for rec in SeqIO.parse("HY12.fasta", "fasta"))
    print(gc_values)


if __name__ == '__main__':
    main()