搜索
查看: 14620|回复: 9

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

[复制链接]

20

主题

68

帖子

870

积分

版主

Rank: 7Rank: 7Rank: 7

积分
870
QQ
发表于 2016-10-11 12:53:23 | 显示全部楼层 |阅读模式
本帖最后由 bioinfo.dong 于 2016-10-13 03:29 编辑

【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序列。

C9130D89-F117-4B45-A5C1-36237F031983.png


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






You really shouldn't spend your time reinventing the wheel
回复

使用道具 举报

0

主题

1

帖子

37

积分

新手上路

Rank: 1

积分
37
发表于 2017-8-14 09:33:12 | 显示全部楼层
程序的18到21行,两个整数做除法,结果很容易变成零,将其中一个数转为浮点数就好了
回复 支持 3 反对 0

使用道具 举报

0

主题

7

帖子

85

积分

注册会员

Rank: 2

积分
85
发表于 2017-11-12 15:23:21 | 显示全部楼层
根据楼主的代码,依葫芦画瓢写了一段,也是第一次知道字典原来也能排序。然后也附上自己思路写的一段代码,仅供参考,因为不太熟悉一些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])

回复 支持 1 反对 1

使用道具 举报

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-7-30 12:11:55 | 显示全部楼层
大神 我用你的代码没跑成功,原来fasta 格式是
>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT

读出来的是
[('Rosalind_6404', 'CCTGCGGAAGATCGGCACTAGAATAGC CAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCCTCCCACTAATAATTCTGAGG'), ('Rosalind_5959', 'CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCTATATCCATTTGTCAGCAGACACGC'), ('Rosalind_0808', '')]
最后 Rosalind_0808,后面的序列读不出来
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复 支持 反对

使用道具 举报

3

主题

34

帖子

318

积分

中级会员

Rank: 3Rank: 3

积分
318
发表于 2017-8-15 10:57:38 | 显示全部楼层
看代码感觉OrderedDict() 跟dict()区别不大啊,就是最后排序的时候简单了一点点
回复 支持 反对

使用道具 举报

0

主题

1

帖子

71

积分

注册会员

Rank: 2

积分
71
发表于 2017-9-9 11:50:19 | 显示全部楼层
请问这里 Orderdict  和dict 有什么区别
回复 支持 反对

使用道具 举报

0

主题

7

帖子

93

积分

注册会员

Rank: 2

积分
93
发表于 2017-10-10 15:45:53 | 显示全部楼层
因为直接输出最大值和序列名称,所以没有排序直接用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)])
回复 支持 反对

使用道具 举报

8

主题

26

帖子

157

积分

注册会员

Rank: 2

积分
157
发表于 2017-11-4 00:16:27 | 显示全部楼层
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)
回复 支持 反对

使用道具 举报

0

主题

19

帖子

68

积分

注册会员

Rank: 2

积分
68
发表于 2020-6-3 11:34:16 | 显示全部楼层
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])
回复 支持 反对

使用道具 举报

0

主题

10

帖子

52

积分

注册会员

Rank: 2

积分
52
发表于 2021-2-6 17:10:01 | 显示全部楼层
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()
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2021-3-5 19:10 , Processed in 0.035545 second(s), 33 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.