搜索
查看: 5649|回复: 12

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

[复制链接]

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-9-23 10:54:09 | 显示全部楼层 |阅读模式
这个很简单啦,首先明白什么是fasta,然后写一个简单脚本即可,正好论坛有人问过这个问题,我转述一下!
FASTA序列如下:
>Rosalind_6404
CCTGCGGAAGATCGGCACTAGAATAGCCAGAACCGTTTCTCTGAGGCTTCCGGCCTTCCC
TCCCACTAATAATTCTGAGG
>Rosalind_5959
CCATCGGTAGCGCATCCTTAGTCCAATTAAGTCCCTATCCAGGCGCTCCGCCGAAGGTCT
ATATCCATTTGTCAGCAGACACGC
>Rosalind_0808
CCACCCTCGTGGTATGGCTAGGCATTCAGGAACCGGAGAACGCTTCAGACCAGCCCGGAC
TGGGAACCTGCGGGCAGTAGGTGGAAT

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

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






上一篇:免疫相关的CD分子有多少个?分布如何?
下一篇:基因表达数据库
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

1

主题

15

帖子

180

积分

注册会员

Rank: 2

积分
180
发表于 2016-9-26 20:27:40 | 显示全部楼层
我来一个菜鸟水准的
[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
回复 支持 1 反对 0

使用道具 举报

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2016-9-28 09:51:02 | 显示全部楼层
disheng 发表于 2016-9-26 20:27
我来一个菜鸟水准的
$ cat test_3
>salind_6404

帮你加上了一个高亮~
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 1 反对 0

使用道具 举报

58

主题

103

帖子

756

积分

版主

Rank: 7Rank: 7Rank: 7

积分
756
QQ
发表于 2016-9-23 16:02:18 | 显示全部楼层
这好像是一个网站上的练习题
回复 支持 反对

使用道具 举报

21

主题

51

帖子

176

积分

信息监察员

Rank: 9Rank: 9Rank: 9

积分
176
发表于 2016-9-23 17:39:09 | 显示全部楼层
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

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

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
学习的痛苦是一时的,学到的快乐是长久的。
回复 支持 反对

使用道具 举报

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
 楼主| 发表于 2016-9-23 20:51:41 | 显示全部楼层
happylittlepig 发表于 2016-9-23 17:39
1.用length统计碱基个数,出现了错误。
cat tmp.txt |perl -alne '$key=$_ if /\>/;if(/^(T|C|G|A)/){$le{$ ...

panda的那个帖子,的确已经讲得很全面了,我并没有更好的办法,代码量精简不是我的长处
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

0

主题

5

帖子

43

积分

新手上路

Rank: 1

积分
43
发表于 2016-10-12 18:00:42 | 显示全部楼层
disheng 发表于 2016-9-26 20:27
我来一个菜鸟水准的
$ cat test_3
>salind_6404

佩服您,这么几行代码就解决问题了,不像我写了那么长的臭臭的代码。。。
回复 支持 反对

使用道具 举报

0

主题

5

帖子

43

积分

新手上路

Rank: 1

积分
43
发表于 2016-10-12 18:02:32 | 显示全部楼层
disheng 发表于 2016-9-26 20:27
我来一个菜鸟水准的
$ cat test_3
>salind_6404

我好像看懂了,谢谢您,学习啦,大神。大神咱俩认识不,能不能加一下您qq呀!我对您的代码真的是特别崇拜
回复 支持 反对

使用道具 举报

0

主题

5

帖子

135

积分

注册会员

Rank: 2

积分
135
发表于 2016-10-22 21:50:46 | 显示全部楼层
Done.
perl的代码好精简啊!
回复 支持 反对

使用道具 举报

0

主题

6

帖子

131

积分

注册会员

Rank: 2

积分
131
发表于 2018-4-25 15:09:00 | 显示全部楼层
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]]
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-13 03:55 , Processed in 0.042620 second(s), 43 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.