搜索
查看: 2610|回复: 2

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

[复制链接]
回帖奖励 20 金钱 回复本帖可获得 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;
                }
        }
}
回复 支持 反对

使用道具 举报

2

主题

23

帖子

169

积分

注册会员

Rank: 2

积分
169
发表于 2019-6-17 00:40:26 | 显示全部楼层

回帖奖励 +5 金钱

python代码
from Bio.Seq import Seq
f=open('对给定的motif序列匹配 fasta文件,并给出定位信息.txt')
dict1={}
for line in f:
    if line.startswith('>'):
        key=line[1:-1]
        dict1[key]=''
    else:
        line.strip()
        dict1[key]+=line
with open('对给定的motif序列匹配 fasta文件,并给出定位信息output.txt','w') as f1:
    for k,v in dict1.items():
        if 'GTGTA' in v:
            r=k+'\t'+'GTGTA'+'\t'+'+'+'\t'+str(v.index('GTGTA')+1)+'\t'+str(v.index('GTGTA')+len('GTGTA'))
            f1.write(r)
    for k,v in dict1.items():
        seq=Seq(v)
        reversed_complement_seq=seq.reverse_complement()
        str_reversed_complement_seq=str(reversed_complement_seq)
        if 'GTGTA' in reversed_complement_seq:
            r = k + '\t' + 'GTGTA' + '\t' + '+' + '\t' + str(str_reversed_complement_seq.index('GTGTA') + 1) + '\t' + str(str_reversed_complement_seq.index('GTGTA') + len('GTGTA'))
            f1.write(r)
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-23 08:24 , Processed in 0.032570 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.