搜索
楼主: 深藏功与名

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

[复制链接]

0

主题

5

帖子

43

积分

新手上路

Rank: 1

积分
43
发表于 2016-10-11 00:32:24 | 显示全部楼层
小饭盆 发表于 2016-10-11 00:28
#发现大家写的都很好,我的好low。本来我看这个题的时候,是没有人回复的。谁知道我做了这么几天的题,就这 ...

截图显示不了,我也不知道为啥。。。
回复 支持 反对

使用道具 举报

0

主题

5

帖子

43

积分

新手上路

Rank: 1

积分
43
发表于 2016-10-11 12:40:54 | 显示全部楼层
挺好的,学习了
回复 支持 反对

使用道具 举报

4

主题

23

帖子

628

积分

高级会员

Rank: 4

积分
628
发表于 2017-12-14 15:37:53 | 显示全部楼层
#!/usr/bin/perl -w
#"GC含量最高的序列"
use strict;
print "ID\tGenelength\tGCcontent\n";
open FA,"$ARGV[0]" or die $!;
$/=">";
<FA>;
my (@GC,$gc_con);
while(<FA>){
    next if(/^\s+$/);
    chomp;
    my ($id,$seq)=(split /\n/,$_,2)[0,1];
    $seq=~ s/\n//g;
    my $len=length($seq);
    my $gc=($seq=~ s/G/G/g+$seq=~ s/C/C/g);
    $gc_con=$gc/$len;
    push @GC,$gc_con;
    print "$id\t$len\t$gc_con\n";
}

my @sort_GC=sort {$a <=> $b} @GC;
my $maxGC=$sort_GC[-1];
print "maxGC_content:", $maxGC,"\n";

close FA;
回复 支持 反对

使用道具 举报

1

主题

55

帖子

832

积分

高级会员

Rank: 4

积分
832
发表于 2017-12-14 17:15:24 | 显示全部楼层
本帖最后由 生信小小菜鸟 于 2017-12-14 17:40 编辑

[Perl] 纯文本查看 复制代码
 use strict;
 open IN,"<","sequence.fa";
 open ONLYB,">","results.fa";
 
 my ($max,$seq_total,@array,%hash);
 local $/=">";
 while (<IN>){
         #chomp;
         my ($ID,$seq1,$seq2)= (split /\n/)[0,1,2];                                
         $seq_total = $seq1.$seq2;                                                        
         my $A = ($seq_total =~tr/A/A/);
         my $T = ($seq_total =~tr/T/T/);
         my $C = ($seq_total =~tr/C/C/);
         my $G = ($seq_total =~tr/G/G/);
         #my $gc = $C+$G;
         unless($gc==0){		$CG = ($C+$G)/($C+$G+$A+$T);
          }
         #print "$gc\n";
         $hash{$gc} = $_;
         }
  foreach (sort {$a<=>$b} keys %hash){    
        push @array,$_;                                
  }
   $max = pop @array;
  printf ONLYB "$hash{$max}\nGC含量:%.2f%",$max;
  close IN;
  close ONLYB;
回复 支持 反对

使用道具 举报

8

主题

26

帖子

157

积分

注册会员

Rank: 2

积分
157
发表于 2017-12-17 22:00:09 | 显示全部楼层
path=r"C:\Users\WIN7\Desktop\将fasta文件中的序列一行转为多行\GC.fasta"
d={}
num=0
with open (path,"r") as f:
    for line in f.readlines():
        if ">" in line:
            name=line.strip()[1:]
            if line not in d:
                d[name]=""
        else:
            d[name]+=line.strip()
for key,value in d.items():
    print("%s GC含量为:%s"%(key,(value.count('C')+value.count('G'))/len(value)))
回复 支持 反对

使用道具 举报

0

主题

2

帖子

85

积分

注册会员

Rank: 2

积分
85
发表于 2018-3-2 09:54:34 | 显示全部楼层
本帖最后由 代码小推车 于 2018-3-2 09:59 编辑

[Perl] 纯文本查看 复制代码
#!perl -w
use strict;
use warnings;

my ($title,$sequence,$GC,$percent);
my ($max,$max_name,$max_gc_sequence) = (0,"","");
while(<DATA>){
    chomp;
    # 假如只是判断行首第一个是不是>,可以用index,使用正则文件数目多了时间明显长于使用index
    if(index($_,">") == 0){
        if($sequence){
            $GC = length((($sequence =~ tr/atAT/ /r) =~ s/ //gr));
            $percent = sprintf "%.2f",($GC * 100 / length($sequence));
            if($max < $percent){$max = $percent;$max_name = $title;$max_gc_sequence = $sequence;}
        }
        $sequence = "";	
        $title = $_;
    }elsif($_){   #防止空行
        $sequence .= $_;	
    }
}
{
# 由于文件末尾的时候,其实最后一条序列没有被考虑,不知道这里有没有好的办法
    $GC = length((($sequence =~ tr/atAT/ /r) =~ s/ //gr));
    $percent = sprintf "%.2f",($GC * 100 / length($sequence));
    if($max < $percent){$max = $percent;$max_name = $title;$max_gc_sequence = $sequence;}
}
print "title:\t${max_name}\nGC:\t${max}%\nsequence:\t$max_gc_sequence\n";

__DATA__;

>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT
回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2018-4-25 15:02:41 | 显示全部楼层
小白第一次评论,在自己写的过程中遇到了一些问题,学习了楼上几人的代码,也算磕磕绊绊完成了任务。
from collections import OrderedDict
with open(r"C:\Users\Lizhiwei\Desktop\新建文本文档.txt") as f :
    d = OrderedDict()
    for line in f :
        if line.startswith(">"):
            seqname = line.strip()
            d[seqname] = ""
            continue
        else:
            d[seqname] += line.upper().strip()
    print d      
##sort函数不能给d.items()排序,sorted函数可以            
    d_count_gc = {}
    for key,value in d.items():
        d_count_gc[key] = ( value.count("G") + value.count("C") ) *1.0 / len(value)
    ##d_count_gc.sort(reverse=True ,key = lambda d : d[1])  错误是字典没有sort方法
    ##d_count_gc.items().sort(reverse=False ,key = lambda d : d[1])  也出不来结果
    print d_count_gc
    d_sort = sorted(d_count_gc.items(),reverse=True,key = lambda d_count_gc : d_count_gc[1])
    print d_sort[0][0]
    print "GC含量最高的序列是", d[d_sort[0][0]]
回复 支持 反对

使用道具 举报

9

主题

15

帖子

244

积分

中级会员

Rank: 3Rank: 3

积分
244
发表于 2018-6-1 10:20:27 | 显示全部楼层
本帖最后由 disil 于 2018-6-1 11:27 编辑

不会调包,方法比较糙
结果输出  >Rosalind_0808 0.6091954022988506

[Python] 纯文本查看 复制代码
d={}#将每个reads转换成字典的一对键、值。参考生信宝典的思路
for i in open('d:/fasta.txt'):
    if i[0] == '>':
        key = i.strip()
        d[key] = []
    else:
        d[key].append(i.strip())
for i ,j in d.items():
    d[i]=''.join(j)
e={}#将d中的值替换为CG含量
for i,j in d.items():
    e[(j.count('C')+j.count('G'))/len(j)]=i
print(e[max(e.keys())],max(e.keys()))
f.close()
回复 支持 反对

使用道具 举报

0

主题

4

帖子

67

积分

注册会员

Rank: 2

积分
67
发表于 2018-6-11 14:58:40 | 显示全部楼层
learnyoung 发表于 2016-10-4 15:40
看了0号菜鸟的逻辑,也写了个python版本(python代码可读性真的好[mw_shl_code=python,true]import os
...

你好,我用的是Python3版本,照你写的程序运行后d[name]=d[name]+line.rstrip()这一行报错,显示KeyErroe:'',请问您知道这个问题该怎么解决吗,初学python,还请多多指教,谢谢!
回复 支持 反对

使用道具 举报

0

主题

4

帖子

67

积分

注册会员

Rank: 2

积分
67
发表于 2018-6-11 16:33:09 | 显示全部楼层
使用Python中的biopython,终于写出可以完美运行的,而且代码比较短。
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC
from Bio.SeqUtils import GC
from Bio import SeqIO
d={}
t=[]
for seqrecord in SeqIO.parse("F:/GC.txt","fasta"):
    d[seqrecord.id]=seqrecord.seq
print(d)
for k,v in d.items():
    d[k]=GC(v)
d_sort=sorted(d.items())
print(d,d_sort)
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-10-23 06:45 , Processed in 0.097623 second(s), 22 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.