免费注册 查看新帖 |

Chinaunix

  平台 论坛 博客 文库
最近访问板块 发新帖
查看: 2942 | 回复: 5
打印 上一主题 下一主题

perl 比对,出不来结果,不知如何修改 [复制链接]

论坛徽章:
0
跳转到指定楼层
1 [收藏(0)] [报告]
发表于 2015-01-26 10:42 |只看该作者 |倒序浏览
<IN1>:>AT1G50920.1 | [1-104] [403-416] | chr1:18870139-18870554 FORWARD LENGTH=118
GGTTGGGTCAGAAAACAGTCGACCGTCATTTGGTTGGGTCGGCTCGCTTCTCTTATCAGCTAACCCTAAGTACGCCGACTTTCAAATTCGC
CACTCCCGTTCATCTTTTCTCTGCATCGATAAGCGAAG
>AT1G15970.1 | [1-374] | chr1:5488495-5488868 REVERSE LENGTH=374
TCGTTTCGTCGTCGATCAGAAAAAGCTTCTTTCTTTCTTTTTTTTTTCTGGGTTTTTCTGAAGTCAAAAGCTTTTTTTTT
CTTCATTGTTCAAAGATTTTATTTTGTTGCTCGCTTGGGTCAACCTTGCTTTTCTCGTGAACACAGCAAGCAACTTGGTT
TAGTTTTACTAGAAGATTCCTAAAAATGAGATATTTAGCGGTGACAAAAAAGTTTCATTTTTGGTGCTAATGTAGCTGTA
AAGTCTAGGCCTTCCGATATATAATCTCTTGAAATAAGTTTCGAGACTGTTGGAAATAGCAATTTTGTTGTTGTTGTGTG
ATCCAGTGATAGATCTGTGAAATCAAAAAGCAGTTTTTTGTGAGAATCGGGTCA
>AT1G73440.1 | [1-55] | chr1:27611363-27611417 FORWARD LENGTH=55
ATTGAAAAGAAAACACATCCCAACTCTGGAGAAACAAATCCCTAATTAGCTAACA
>AT1G75120.1 | [1-146] | chr1:28198657-28198802 REVERSE LENGTH=146
CCGTGTTAAAAAACACGTTGGGTCAAAATCAATGGTGGGTGGGTCTTAACTTGCACTACTAGGATTTGATCACAGAAGAAA
GAAAGGGGTTTATTAATTTGGTTACTGTATCTTACTCTTTTTACTTTCGGATAACTTCAATTGGCA
<IN>:GGTTGGGTC
[ATC]GTTGGGTC
G[ATC]TTGGGTC
GG[ACG]TGGGTC
GGT[ACG]GGGTC
GGTT[ATC]GGTC
GGTTG[ATC]GTC
GGTTGG[ATC]TC
GGTTGGG[ACG]C
GGTTGGGT[AGT]
目的:将<IN>中的每行序列各自比对到<IN1>中所有基因序列上,将比对到的输出,输出位置信息及相应的比对序列
我的程序:open IN1,"<seq_list.txt" ;
while (<IN1>
{
chomp;
my @aa1=split /\s+/,$_;
my $a=$aa1[0];
nn($a);
}
close IN1;
sub nn{
        my @a=$_;
open IN,"<test.fa" ;
open OUT,">>test_out.txt";
my %hash;
my $seq;
$/=">";
while (<IN>
{
chomp;
my @aa=split /\n/;
my $gene=(split /\s+/,shift @aa)[0];
$hash{$gene}=$_;
$seq=join('',@aa[0..$#aa]);
while ($seq =~ m/$a[0]/g)
{
my $pos = pos $seq;
push (@pos,$pos-;
pos($seq)=$pos-8;
if (pos($seq)){
print OUT "$gene\t", pos $seq, "\n";}
}
}
close IN;
}
出来的结果和我想要的不太一样:AT1G50920.1        1        GGTTGGGTC
                                          AT1G50920.1        32        GGTTGGGTC
实际上还有结果没有输出

论坛徽章:
0
2 [报告]
发表于 2015-01-26 12:13 |只看该作者
额,没有人会吗

论坛徽章:
8
技术图书徽章
日期:2013-08-22 11:21:28未羊
日期:2015-01-19 22:22:25巳蛇
日期:2014-08-11 16:53:08子鼠
日期:2014-05-29 09:04:44摩羯座
日期:2014-04-11 14:15:07丑牛
日期:2014-01-24 12:41:28金牛座
日期:2013-11-21 17:38:28射手座
日期:2015-01-21 08:50:32
3 [报告]
发表于 2015-01-26 13:14 |只看该作者
本帖最后由 huang6894 于 2015-01-26 13:20 编辑

回复 2# 1002279625

  1. my($fa, $seq) = @ARGV;

  2. my @all_seq;
  3. open IN,"< $seq" or die "$!";
  4. while(<IN>){
  5.     chomp;
  6.     if(/^\[([ATCG]+)\]([ATCG]+)$/){
  7.         map{push @all_seq,$_.$2}(split //,$1);
  8.     }elsif(/^([ATCG]+)\[([ATCG]+)\]([ATCG]+)$/){
  9.         map{push @all_seq,$1.$_.$3}(split //,$2);
  10.     }elsif(/^([ATCG]+)\[([ATCG]+)\]$/){
  11.         map{push @all_seq,$1.$_}(split //,$2);
  12.     }else{
  13.         push @all_seq,$_;
  14.     }
  15. }
  16. close IN;
  17. $/ = '>';
  18. open IN1,"< $fa" or die "$!";
  19. <IN1>;
  20. while(<IN1>){
  21.     chomp;
  22.     my ($key, @seq) = split /\n/;
  23.     my ($gene) = $key =~/^([^\|]+)(\s+)?|/;
  24.     my $s =  join("",@seq);
  25.     foreach my $key2(@all_seq){
  26.         my ($pos, $now) = (0, -1);
  27.         until ( $pos eq -1 ) {
  28.             $pos = index( $s, $key2, $now + 1 );
  29.             $now = $pos;
  30.             print $gene."\t".$pos."\t".$key2."\n" unless $pos < 0;
  31.         }
  32.     }
  33. }
  34. close IN2;
  35. $/ = "\n";
复制代码
结果:
  1. AT1G50920.1     0       GGTTGGGTC
  2. AT1G50920.1     31      GGTTGGGTC
  3. AT1G15970.1     112     GCTTGGGTC
  4. AT1G75120.1     15      CGTTGGGTC
  5. AT1G75120.1     36      GGGTGGGTC
复制代码
如果数据量小也就罢了,太大,可以考虑用kmp算法,或者直接blat

论坛徽章:
0
4 [报告]
发表于 2015-01-26 13:29 |只看该作者
回复 3# huang6894


    我是seq_list这个比较长,有116行,考虑到9bp中最多有两个错配。用过blat,但因为9bp太短,出不来结果,所以才自己写程序比对。

论坛徽章:
0
5 [报告]
发表于 2015-01-26 16:02 |只看该作者
回复 3# huang6894
我发现,这样写,有漏掉的


   

论坛徽章:
0
6 [报告]
发表于 2015-01-26 16:16 |只看该作者
回复 3# huang6894
留着以后慢慢研究大神程序~

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

本版积分规则 发表回复

  

北京盛拓优讯信息技术有限公司. 版权所有 京ICP备16024965号-6 北京市公安局海淀分局网监中心备案编号:11010802020122 niuxiaotong@pcpop.com 17352615567
未成年举报专区
中国互联网协会会员  联系我们:huangweiwei@itpub.net
感谢所有关心和支持过ChinaUnix的朋友们 转载本站内容请注明原作者名及出处

清除 Cookies - ChinaUnix - Archiver - WAP - TOP