搜索
查看: 3579|回复: 20

[新人求教]找出CG含量最大的DNA片段

[复制链接]

4

主题

29

帖子

119

积分

注册会员

Rank: 2

积分
119
发表于 2016-9-19 14:05:18 | 显示全部楼层 |阅读模式
原问题如下:
In Rosalind's implementation, a string in FASTA format will be labeled by the ID "Rosalind_xxxx", where "xxxx" denotes a four-digit code between 0000 and 9999.

Given: At most 10 DNA strings in FASTA format (of length at most 1 kbp each).

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.
Sample Dataset

>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT

Sample Output

Rosalind_0808
60.919540

这里主要麻烦的一是处理首行,这个我是用正则表达式做的。其次是找出GC含量最大的片段,因为我把数据放在了hash里,感觉找起来比较麻烦,希望大神能优化一下这部分算法。

这是我的解答:
[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w

use strict;

my %seq;
my ($n, $out, $rato, @m);
$rato = 0;

while(<>) {
	chomp;
	if(/^>(\w+)/) {
		$n = $1;
		next;
	}
	my @list = split(//, $_);
	for my $word (@list) {
		$seq{$n}{$word}++;
		$seq{$n}{'total'}++;
	}
}

for my $lable ( keys %seq ) {
	my $half = $seq{$lable}{'G'} + $seq{$lable}{'C'};
	$seq{$lable}{'rato'} = $half / $seq{$lable}{'total'}; 
	if($seq{$lable}{'rato'} > $rato) {
		$rato = $seq{$lable}{'rato'};
		$out = $lable
	}
}
my $per = 100 * $seq{$out}{'rato'};
print "$out\n$per\n";





上一篇:把hg19基因组按照染色体切割
下一篇:如何提取染色体一段区域内的所有基因的注释信息
回复

使用道具 举报

0

主题

14

帖子

225

积分

中级会员

Rank: 3Rank: 3

积分
225
发表于 2016-9-28 11:07:41 | 显示全部楼层
我也是新手~  python写的,欢迎相互交流~
[Python] 纯文本查看 复制代码
# coding=utf-8
import os

os.chdir('C:/workspace')

fh = open('test.txt')
read_fh = fh.readlines()
number = 0
dic = {}

#整理数据格式
for line in read_fh:
    number += 1
    if line[0] == '>':
        name = line[:-1]
        dic[name] = ''
        continue
    else:
        dic[name] = dic[name]+line[:-1]

#计算GC含量
for x in dic.keys():
    count_all = len(dic[x])
    count_Gg = dic[x].count('G')+dic[x].count('g')
    count_Cc = dic[x].count('C')+dic[x].count('c')
    dic[x] = (count_Cc+count_Gg)*1.0/count_all
dic_sort = sorted(dic.iteritems(),key=lambda d:d[1],reverse=True)
print 'most GC ratio gene is %s and it is %s ' %(dic_sort[0][0],dic_sort[0][1])

fh.close()


回复 支持 2 反对 0

使用道具 举报

1

主题

41

帖子

282

积分

中级会员

Rank: 3Rank: 3

积分
282
发表于 2016-10-4 15:40:04 | 显示全部楼层
看了0号菜鸟的逻辑,也写了个python版本(python代码可读性真的好
[Python] 纯文本查看 复制代码
import os
fh = open('fasta.txt')
read_fh = fh.readlines()
d= {}
t=[]
# 整理数据格(数据名称:数据(ATCG...))
for line in read_fh:
    if line.startswith('>'):
        name = line.rstrip()
        d[name] = ''
        continue
    else:
        d[name] = d[name] + line.rstrip()
for k,v in d.items():
    GC_content=(v.count('C')+v.count('G'))*1.0/len(v)
    t.append((k,GC_content))
t.sort(reverse=False)
print 'max GC_content squence is %s:%s'%(t[0][0],t[0][1])
fh.close()
回复 支持 1 反对 0

使用道具 举报

20

主题

68

帖子

862

积分

版主

Rank: 7Rank: 7Rank: 7

积分
862
QQ
发表于 2016-10-2 14:29:14 | 显示全部楼层
本帖最后由 bioinfo.dong 于 2016-10-2 14:32 编辑

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

使用道具 举报

20

主题

68

帖子

862

积分

版主

Rank: 7Rank: 7Rank: 7

积分
862
QQ
发表于 2016-10-2 14:33:17 | 显示全部楼层
0号菜鸟 发表于 2016-9-28 11:07
我也是新手~  python写的,欢迎相互交流~[mw_shl_code=python,true]# coding=utf-8
import os

Your logic is pretty good but you can polish the code more. There are many built-in functions in Python like startswith(), upper() which could make your life easier. Here I also paste my script based on yours. Let's crack Python! Fighting!

You really shouldn't spend your time reinventing the wheel
回复 支持 反对

使用道具 举报

20

主题

68

帖子

862

积分

版主

Rank: 7Rank: 7Rank: 7

积分
862
QQ
发表于 2016-10-2 14:34:52 | 显示全部楼层
本帖最后由 bioinfo.dong 于 2016-10-2 14:36 编辑
0号菜鸟 发表于 2016-9-28 11:07
我也是新手~  python写的,欢迎相互交流~

[Python] 纯文本查看 复制代码
import os
from operator import itemgetter
from collections import OrderedDict

seqTest = OrderedDict()
gcContent = OrderedDict()

os.chdir('C:/workspace')

with open('test.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)
    
sortedGCContent = sorted(gcContent.items(), key = itemgetter(1))

largeName = sortedGCContent[-1][0]
largeGCContent = round(sortedGCContent[-1][1],2)

print ('most GC ratio gene is %s and it is %s ' %(largeName,largeGCContent))
You really shouldn't spend your time reinventing the wheel
回复 支持 反对

使用道具 举报

0

主题

14

帖子

225

积分

中级会员

Rank: 3Rank: 3

积分
225
发表于 2016-10-4 08:32:57 | 显示全部楼层
bioinfo.dong 发表于 2016-10-2 14:34
[mw_shl_code=python,true]import os
from operator import itemgetter
from collections import OrderedD ...

Thank you for your suggestion, your little tips make me feel refreshed. Let's crack Python! Fighting!
回复 支持 反对

使用道具 举报

1

主题

41

帖子

282

积分

中级会员

Rank: 3Rank: 3

积分
282
发表于 2016-10-4 15:37:01 | 显示全部楼层
0号菜鸟 发表于 2016-9-28 11:07
我也是新手~  python写的,欢迎相互交流~[mw_shl_code=python,true]# coding=utf-8
import os

数据整理的逻辑好棒
回复 支持 反对

使用道具 举报

4

主题

51

帖子

327

积分

中级会员

Rank: 3Rank: 3

积分
327
发表于 2016-10-8 23:36:02 | 显示全部楼层
试问一下楼主rosalind刷到多少了?我最近也在刷。可以交流交流。
回复 支持 反对

使用道具 举报

0

主题

5

帖子

43

积分

新手上路

Rank: 1

积分
43
发表于 2016-10-11 00:28:39 | 显示全部楼层
#发现大家写的都很好,我的好low。本来我看这个题的时候,是没有人回复的。谁知道我做了这么几天的题,就这么多回复了。自己写的太low。但是写了,我就拿出来。
#!perl -w
open FH,"<","fasta.txt"or die $!;
@files=<FH>;#把文件全部读到数组里面
$i=0;
$name=0;
$max=0;
while($i<@files){#只要数组没有越界,就是说文件还没读完,就一直接着读
$q=0;
$w=0;
if($files[$i]=~m/^>/){#如果以大于号开头说明是基因名字,
  $temp=$files[$i];#先把名字暂时存到temp变量中
  #print $temp;
  $i++;
  if($i>=@files){#还是检查越界问题
   last;
  }
}
while(!($files[$i]=~m/^>/)){#如果以不是大于行开头,意思是是碱基序列的话,
  @nn=split//,$files[$i];#就切割,然后放到数组里面
  #$w+=@nn;#这一行不能这样写,这样得出的结果是不对的。个人感觉是切割的问题,比如存在空格,他会当成一个元素,以至于总量w的计算出现问#题,为保险起见,我就有了下面的if语句
  for $p(@nn){
   if($p eq C || $p eq G){#如果是C和G,那么
    $q++;#变量q是C和G的含量,
    $w++;#变量w是总的碱基含量。
   }
   elsif($p eq A || $p eq T){
    $w++;#匹配到ATCG都加
   }
  }
  $i++;
  if($i>=@files){#检查数组是否越界
   last;
  }
}
#print $w;
#print $q," 呵呵 ";
$ratio=$q/$w;#计算CG百分含量

#print $ratio,"        ";
if($max<$ratio){#找CG含量最大的
  $max=$ratio;
  #print $max,"!!";
  $name=$temp;
  #print $name,"....";
}
if($i>=@files){#考虑数组越界
  last;
}

}
print $name,"\n";#输出最终的基因名字
print $max;
正确答案截图:

还有一个错误与大家共享,自己看(代码还是原来的代码没有任何变动,害惨我了)

其实文件里面不匹配行的绝对开头就更好了。完美,睡觉
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-9-22 04:30 , Processed in 0.138098 second(s), 32 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.