搜索
查看: 1903|回复: 1

<PERL> 如何对给定的motif序列匹配 fasta文件,并给出定位信息

[复制链接]
回帖奖励 25 金钱 回复本帖可获得 5 金钱奖励! 每人限 1 次

1

主题

4

帖子

140

积分

注册会员

Rank: 2

积分
140
发表于 2017-9-28 16:04:51 | 显示全部楼层 |阅读模式
最近刚刚开始学Perl语言,向大家请教道题:

给定的fasta文件:分别为染色体名称,及序列信息,如
>chr_1
ATGCGAGAGAGAGAGAAGTGCTGTGTAGCTGATGCGCTAGTTTCGCGCTAGAGA....................

现在想对给定的 短序列motif进行定位如对 GTGTA 定位,能分别对正反向链 所有的motif匹配定位,要求给出具体起始和终止的位置信息
输出结果例如为:
chr_1    GTGTA     +     23     27
chr_1    GTGTA     +     23550     23554
chr_1    GTGTA     -      33333    33337
chr_1    GTGTA     -      323333    323337。。

实在缺乏思路,想向大家请教下思路跟具体代码。谢谢!!



上一篇:比对工具代码大全
下一篇:16S扩增子数据拆分问题
回复

使用道具 举报

1

主题

4

帖子

140

积分

注册会员

Rank: 2

积分
140
 楼主| 发表于 2017-10-1 23:50:16 | 显示全部楼层
纠结了几天终于在翻小骆驼时找到了解决答案,利用index函数。
初学者代码如下,如若有更好的或简化的方法,万分感谢告知!

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
open QUERY, $ARGV[0] or die "$!";
open IN, $ARGV[1] or die "$!";
chomp($query=<QUERY>);
$len=length$query;
while (<IN>){
        chomp;
        if(/^>/){
                $id=$_;
        }else{
                ($pos,$now)=(0,-1);
                until($pos == -1){
                        $pos=index($_,$query,$now+1);
                        $now=$pos;
                        print "$id\t$query\t+\t$pos\t",($pos+$len-1),"\n" unless $pos<0;
                }
                $reverse_seq=reverse($_);
                $reverse_seq=~tr/ATGCatgc/TACGtacg/;
                ($pos_re,$now_re)=(0,-1);
                until($pos_re== -1){
                        $pos_re=index($reverse_seq,$query,$now_re+1);
                        $now_re=$pos_re;
                        print"$id\t$query\t-\t$pos_re\t",($pos_re+$len-1),"\n" unless $pos_re<0;
                }
        }
}
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-5-21 15:43 , Processed in 0.117980 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.