从fasta文件中找出CG含量最大的DNA片段

Jimmy

这个很简单啦，首先明白什么是fasta，然后写一个简单脚本即可，正好论坛有人问过这个问题，我转述一下！
FASTA序列如下：

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

脚本很简单，就统计每个序列的GC含量即可， 然后输出最大的，用perl一句话命令即可！
说真的，要学好 perl 还真的不简单，因为 perl 的程序代码比C来得精简一半1，靠得就是在撰写时的大脑运作。程设师得花更多时间写出精简的 code，同时也要将「语意上的错误」减少到最低，这就是要靠经验的累积。

所以不建议初学者用perl单行命令，还是好好写吧！

disheng

我来一个菜鸟水准的
[user132@mu01 perl]$ cat test_3
>salind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT
[user132@mu01 perl]$ cat test_3.pl

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

#!/usr/bin/perl 
open IN,"<$ARGV[0]" or die $!;

$/ = ">";<IN>;
my $max_gc;
my $max_id;
while(<IN>){
        s/\r?\n>//;
        my($id,$seq) =split /\n/,$_,2;
        my $length = length $seq;
        my $gc =($seq =~ s/C/C/g) + ($seq =~ s/G/G/g);
        my $GC =$gc/$length*100;
        if($max_gc < $GC){
                $max_gc = $GC;
                $max_id = $id;
        }
#        print "id$id\tlength$length\tgc$gc\tmax_gc$max_gc\n";
}
print "$max_id\t$max_gc\n";
close IN;

[user132@mu01 perl]$ perl test_3.pl test_3
Rosalind_0808        59.5505617977528

Jimmy

disheng 发表于 2016-9-26 20:27
我来一个菜鸟水准的
$ cat test_3
>salind_6404

帮你加上了一个高亮~

Panda姐

这好像是一个网站上的练习题

happylittlepig

1.用length统计碱基个数，出现了错误。
cat tmp.txt |perl -alne '$key=$_ if /\>/;if(/^(T|C|G|A)/){$le{$key}+=length{$_};$GC{$key}+=($_=~tr/GC//);$per{$key}=$GC{$key}/$le{$key}}END{foreach $k(sort keys %per){print"##$k $per{$k} $le{$k}##"}}'

file:///C:/Users/yuli/AppData/Local/YNote/data/happylittlepig@126.com/975b346cd2a84722b5e0df373c1d76fc/clipboard.png
2.改用tr统计碱基个数
cat tmp.txt |perl -alne '$key=$_ if /\>/;if(/^(T|C|G|A)/){$le{$key}+=($_=~tr/TCGA//);$GC{$key}+=($_=~tr/GC//);$per{$key}=$GC{$key}/$le{$key}}END{foreach $k(sort keys %per){print"##$k $per{$k} $le{$k}##"}}'

file:///C:/Users/yuli/AppData/Local/YNote/data/happylittlepig@126.com/cccb3bcec2ac42ce844466bc6a94c9ce/clipboard.png
3.找最大的GC%，用到了一个模块中的max函数。
cat tmp.txt |perl -MList::Util=max -alne '$key=$_ if /\>/;if(/^(T|C|G|A)/){$le{$key}+=($_=~tr/TCGA//);$GC{$key}+=($_=~tr/GC//);$per{$key}=$GC{$key}/$le{$key}}END{foreach $k(sort keys %per){print"##$k $per{$k} $le{$k}##"} print "##max GC% is ##\n"; print max values %per}'

file:///C:/Users/yuli/AppData/Local/YNote/data/happylittlepig@126.com/fdc5ced4edbe4b05b4d661eb65ce3ac8/clipboard.png

如果大家有更好的方法请告诉我哦

Jimmy

happylittlepig 发表于 2016-9-23 17:39
1.用length统计碱基个数，出现了错误。
cat tmp.txt |perl -alne '$key=$_ if /\>/;if(/^(T|C|G|A)/){$le{$ ...

panda的那个帖子，的确已经讲得很全面了，我并没有更好的办法，代码量精简不是我的长处

小饭盆

disheng 发表于 2016-9-26 20:27
我来一个菜鸟水准的
$ cat test_3
>salind_6404

佩服您，这么几行代码就解决问题了，不像我写了那么长的臭臭的代码。。。

小饭盆

disheng 发表于 2016-9-26 20:27
我来一个菜鸟水准的
$ cat test_3
>salind_6404

我好像看懂了，谢谢您，学习啦，大神。大神咱俩认识不，能不能加一下您qq呀！我对您的代码真的是特别崇拜

经典蓝色

Done.
perl的代码好精简啊！

gevin_li

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]]