再求个perl程序，本人分不多，但是需要紧急，望各位帮助小弟！

aids260

这个问题好像前边问过但是现在格式不一样了啊

有这样的些位置信息

1. BGIBMGA000062-TA         nscaf1071        1239531
   
2. BGIBMGA000064-TA         nscaf1071        1313194
   
3. BGIBMGA000129-TA         nscaf1108        1283079
   
4. BGIBMGA000140-TA         nscaf1108        1773552
   
5. BGIBMGA000154-TA         nscaf1108        2355076
   
6. BGIBMGA000155-TA         nscaf1108        2429759
   
7. BGIBMGA000201-TA         nscaf1299        349886
   
8. BGIBMGA000202-TA         nscaf1375        7946
   
9. BGIBMGA000360-TA         nscaf1681        45390
   
10. BGIBMGA000362-TA         nscaf1681        90916
    
11. BGIBMGA000367-TA         nscaf1681        262818
    
12. BGIBMGA000372-TA         nscaf1681        801566
    
13. BGIBMGA000381-TA         nscaf1681        1125109
    
14. BGIBMGA000384-TA         nscaf1681        1165828
    
15. BGIBMGA000394-TA         nscaf1681        1732614

复制代码

有这样的一些  位置信息，第一列是基因的名字上，第二列是nscaf的名称，第三列是起点,是这个基因比对对到ncaf上的位点
现在我的需要是把 起点上游2000bp  到下游500bp的碱基序列(-2000~+500)     从ncaf上提取出来

nscaf格式如下：

1. >nscaf100 /length=4083 /lengthwogaps=4073
   
2. TCTCTTAAAGGCTGCAAGTCCTAAGTCATCTTTTAAAATACGCGACATGGTGCTATCTTCATCTCCCGAG
   
3. ATAAAATCTTTTGCTTTCGGACAGGATTTCTTCGAATTCTTTTCCTTACTGCTTTGACCACCTTTTTCGT
   
4. ACGAACACTACGTGAACGGCCAGATCTCATTCCCGATCCTGCGGACAGAATGGAAAGCAGTCGACGTCGC
   
5. CCAAAACACGTCTCTCGGATCCTCCCGATCCACTAACGGTGCTTTTAGGTACTCCAAGCACCGGTCACCG
   
6. TTCTCGTCGAACCCGTCGCTTGCGACGAAGGGCTCGACGAGTAAATTAACTCTCAGACACA...........................
   
7. .................
   
8. 
   
9. >nscaf1000 /length=20507 /lengthwogaps=15885
   
10. CTCCAAAGATTTATTGGATAATGGTTTTTTTGTTGTTGCTTCCGAAGATGTACTAAAACTTACATAAGTC
    
11. CTTATTTTTGACAAGAAAGACACTTATGCACCTCTTTTAAACTTTATTTGGAGAGGTGCATGTCATGACG
    
12. TGAAAACTCAGAGTAGGTGAGTTCGTGCCATAAGAGATTGGTCGCCAATGATTTTGGTATGCGCATTCTA
    
13. TAAAAGCAGAGATCCGAGGTGGTAGATTTTAGCTCTTGTTCAAGGTATGGGGTAAAAAGTAGACTAAAAT
    
14. ATCACAGGTCTGAAGAACAATTTGAAGAACAATGCCCGTTTGAATTTT.......................................
    
15. .....................................
    

复制代码

提取出的格式最好是这样的：
>pro1/BGIBMGA000062-TA
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
>pro2/BGIBMGA000064-TA
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
>pro3/BGIBMGA000129-TA
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
..........................................................
...........................................................

我把一些格式传为附件了，因为分数有限，就只给最快，最准确的了！

longbow0

输出文件是 'out.fa'.
在 Debian 下已经调试通过。

1. 
   
2. #!/usr/bin/perl
   
3. 
   
4. use strict;
   
5. use warnings;
   
6. 
   
7. use Bio::SeqIO;
   
8. 
   
9. # Read all seqs into a hash
   
10. my $seqi = Bio::SeqIO->new(
    
11.         -file => 'nscaf格式.fa',
    
12.         -format => 'fasta',
    
13. );
    
14. 
    
15. my $rh_seqs;
    
16. 
    
17. while (my $seq = $seqi->next_seq) {
    
18.         $rh_seqs->{$seq->id} = $seq;
    
19. }
    
20. 
    
21. # Output seqs
    
22. my $seqo = Bio::SeqIO->new(
    
23.         -file => ">out.fa",
    
24.         -format => 'fasta',
    
25. );
    
26. 
    
27. my $i = 1;
    
28. 
    
29. open(FH, "位置信息1.txt") or die "Error: '位置信息1.txt' $!\n";
    
30. 
    
31. while (<FH>) {
    
32.         next if /^\s*$/;
    
33.        
    
34.         my ($geneid, $seqid, $pos) = split /\s+/;
    
35. 
    
36.         print 'Gene: ', $geneid, 'Seq: ' . $seqid, "\n";
    
37.        
    
38.         # 提示不存在的 seqid
    
39.         unless ( defined $rh_seqs->{$seqid} ) {
    
40.             warn 'Seq ', $seqid, "does NOT exost!\n";
    
41.             next;
    
42.         }
    
43.         my $gene_str =  $rh_seqs->{$seqid}->subseq($pos-2000, $pos+500);
    
44.        
    
45.         my $sseq = Bio::Seq->new(
    
46.         -id => 'pro' . $i . '/' . $geneid,
    
47.         -seq => $gene_str,
    
48.         );
    
49.        
    
50.         $seqo->write_seq($sseq);
    
51. 
    
52.         $i++;
    
53. }
    
54. 
    
55. close FH;
    

复制代码

aids260

回复 2# longbow0


    哥们高手！

aids260

回复 2# longbow0


     学习了！

leigh111

open(FH,'E:\task\nscaf格式.fa') || die "$!";
open(IN,'E:\task\位置信息.txt') || die "$!";
open(OUT,'>E:\task\nscaf.fasta') || die "$!";

my ($title,@array,%seq);
while (<FH>) {
         chomp;
      if (/^>(nscaf\d+)/) {
         $title = $1;
      }
      else {
         $seq{$title} .= $_;
      }
}
close FH;
my ($up,$down,$num);
while (<IN>) {
   chomp;
   @array = split/\s+/,$_;

   if ($seq{$array[1]}) {
       $num++;
       $up = substr ($seq{$array[1]},$array[2]-2000,2000);
       $down = substr ($seq{$array[1]},$array[2],500);
        print ">pro.$num\/$array[0]\n";
        print $up.$down,"\n";
   }
}

aids260

回复 5# leigh111


    谢谢了啊，哥们！

wuhan_2112

学习了、