2013-05-02 48 views
1

我有不同长度的每个文件中具有1000个seq的多个fasta文件。我想只保留每个序列的前200(n)个碱基。我如何在Perl中做到这一点?fasta:在n长度后删除序列

回答

0

很难准确地理解你的意思,没有看到一个例子,但如果你只需要每行的前200个字符只使用cut

cut -c1-200 file 
+0

只为我打印标题行。我在http://pastebin.com/51nVG5nD创建了一个示例输入文件 – tripleee 2013-05-02 11:06:34

+0

我无法运行这个脚本,而是使用了下面的脚本,它运行流畅:cut -c -200文件。感谢您的帮助 – Ronn 2013-05-02 11:15:12

+0

@Ronn您是否在说我的原创答案'cut -c1-200'解决了您的问题? – 2013-05-02 11:56:57

1

如果序列太长,只保留最有趣的部分:

$/ = '>'; 
<>; 
while (my $seq = <>) { 
    $seq =~ s/>$//; 
    $seq =~ s/^(.*)//; 
    my $id = $1; 
    $seq =~ s/\n//g; 
    $seq = substr $seq, 0, 200; 
    print ">$id\n$seq\n"; 
} 
+0

+1为'$ /'技巧!但是这破坏了换行符。根据http://en.wikipedia.org/wiki/FASTA_format,行应保持在80个字符以下。 – tripleee 2013-05-02 10:31:36

+0

@triplee:这只是一个建议:-)您可以在'print'行之前添加'$ seq =〜s /(。{80})(?=。)/ $ 1 \ n/g;'。 – choroba 2013-05-02 11:00:19

+0

感谢Choroba,它完美的工作 – Ronn 2013-05-02 11:13:51

2

如果序列打印在几个物理线路,只能通过第200个字符打印了。以楔形开始的行是标题行,表示新序列的开始。

awk '/^>/{ seqlen=0; print; next; } 
    seqlen < 200 { if (seqlen + length($0) > 200) 
      $0 = substr($0, 1, 200-seqlen); 
     seqlen += length($0); print }' file.fasta >newfile.fasta 

呵呵,在Perl中?

perl -nle 'if (/^>/) { $seqlen = 0; print; next } 
    next if ($seqlen >= 200); 
    $_ = substr($_, 0, 200-$seqlen) if ($seqlen + length($_) > 200); 
    $seqlen += length($_); 
    print;' file.fasta >newfile.fasta 
+0

感谢Tripleee为您的答案。我试过你的Perl脚本,它的工作。 – Ronn 2013-05-02 11:14:34

1

我建议你考虑使用的BioPerl此排序的事情,因为它很容易完成这些任务,你不必担心格式化等事情。在下面的代码中,脚本的第一个参数是你的fasta,第二个参数是一个仅保存每个序列的前200个基地的文件。

#!/usr/bin/env perl 

use strict; 
use warnings; 
use Bio::Seq; 
use Bio::SeqIO; 

my $usage = "$0 infile outfile\n"; 
my $infile = shift or die $usage; 
my $outfile = shift or die $usage; 

my $seqin = Bio::SeqIO->new(-file => $infile, -format => 'fasta'); 
my $seqout = Bio::SeqIO->new(-file => ">$outfile", -format => 'fasta'); 

while (my $seq = $seqin->next_seq) { 
    my $first200 = $seq->subseq(1,200); # 1-based 
    my $subseq = Bio::Seq->new(-seq => $first200, -id => $seq->id); 
    $seqout->write_seq($subseq); 
} 
0

下面是我如何解决这个问题,如果有兴趣的人尝试了另一种方式来做到这一点 我用包含在biolinux称为Fasta_formatter把实际的序列中的一个线(-w 0)的工具,然后修剪为@sudo_O说,然后最后回到80个字母的宽度。

fasta_formatter -w 0 < FILE | cut -c1-LENGTH | fasta_formatter -w 80 > TRIMMED_FILE