fasta：删除 n 长度后的序列答案

【问题标题】：fasta: delete sequences after n lengthfasta：删除 n 长度后的序列
【发布时间】：2015-07-03 13:04:35
【问题描述】：

我有多个 fasta 文件，每个文件中包含 1000 个不同长度的 seq。我想只保留每个序列的前 200 个（n）个碱基。如何在 Perl 中做到这一点？

【问题讨论】：

【解决方案1】：

如果序列打印在多个物理行上，则只打印到第 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

【讨论】：

感谢 Tripleee 的回答。我试过你的 perl 脚本，它成功了。

【解决方案2】：

如果序列太长，只保留有趣的部分：

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

【讨论】：

+1 $/ 技巧！但这破坏了换行符。根据en.wikipedia.org/wiki/FASTA_format，行数应保持在 80 个字符以下。
@triplee：这只是一个建议 :-) 您可以随时在 print 行之前添加 $seq =~ s/(.{80})(?=.)/$1\n/g;。

【解决方案3】：

我建议您考虑将 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);
}

【讨论】：

【解决方案4】：

如果不看示例，很难准确理解您的意思，但如果您只需要每行的前 200 个字符，请使用 cut：

cut -c1-200 file

【讨论】：

只为我打印标题行。我在pastebin.com/51nVG5nD 中创建了一个示例输入文件
我无法运行这个脚本，而是使用了下面的脚本，运行流畅：cut -c -200 文件。感谢您的帮助
@Ronn 你是说我的原始答案cut -c1-200 解决了你的问题吗？

【解决方案5】：

如果有人有兴趣尝试另一种方法，这就是我的解决方法我使用 biolinux 中包含的一个名为 Fasta_formatter 的工具将实际序列放在一行中（-w 0），然后按照@sudo_O 所说的进行修剪，最后回到 80 个字母的宽度。

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

【讨论】：