编写一个接受 fasta 并反转所有序列的 Perl 脚本(没有 BioPerl)?
Write a Perl script that takes in a fasta and reverses all the sequences (without BioPerl)?
我不知道这是否只是 Stawberry Perl 的一个怪癖,但我似乎无法做到 运行。我只需要做一个 fasta 并反转其中的每个序列。
-问题-
我有一个 multifasta 文件:
>seq1
ABCDEFG
>seq2
HIJKLMN
预期输出为:
>REVseq1
GFEDCBA
>REVseq2
NMLKJIH
脚本在这里:
$NUM_COL = 80; ## set the column width of output file
$infile = shift; ## grab input sequence file name from command line
$outfile = "test1.txt"; ## name output file, prepend with “REV”
open (my $IN, $infile);
open (my $OUT, '>', $outfile);
$/ = undef; ## allow entire input sequence file to be read into memory
my $text = <$IN>; ## read input sequence file into memory
print $text; ## output sequence file into new decoy sequence file
my @proteins = split (/>/, $text); ## put all input sequences into an array
for my $protein (@proteins) { ## evaluate each input sequence individually
$protein =~ s/(^.*)\n//m; ## match and remove the first descriptive line of
## the FATA-formatted protein
my $name = ; ## remember the name of the input sequence
print $OUT ">REV$name\n"; ## prepend with #REV#; a # will help make the
## protein stand out in a list
$protein =~ s/\n//gm; ## remove newline characters from sequence
$protein = reverse($protein); ## reverse the sequence
while (length ($protein) > $NUM_C0L) { ## loop to print sequence with set number of cols
$protein =~ s/(.{$NUM_C0L})//;
my $line = ;
print $OUT "$line\n";
}
print $OUT "$protein\n"; ## print last portion of reversed protein
}
close ($IN);
close ($OUT);
print "done\n";
这会按照你的要求去做
它从 FASTA 文件中构建一个散列 %fasta
,保持数组 @keys
以保持序列的顺序,然后打印出散列的每个元素
序列的每一行在添加到哈希之前使用reverse
反转,使用unshift
以相反的顺序添加序列的行
程序期望输入文件作为命令行参数,并将结果打印到STDOUT,可在命令行重定向
use strict;
use warnings 'all';
my (%fasta, @keys);
{
my $key;
while ( <> ) {
chomp;
if ( s/^>\K/REV/ ) {
$key = $_;
push @keys, $key;
}
elsif ( $key ) {
unshift @{ $fasta{$key} }, scalar reverse;
}
}
}
for my $key ( @keys ) {
print $key, "\n";
print "$_\n" for @{ $fasta{$key} };
}
输出
>REVseq1
GFEDCBA
>REVseq2
NMLKJIH
更新
如果您更喜欢重新包装序列以便短行在末尾,那么您只需要重写转储哈希的代码
该方案以原文件中最长行的长度为限,将反转后的序列重新包装成相同的长度。更清楚的是,指定一个显式长度而不是计算它会很简单
您需要在程序顶部添加 use List::Util 'max'
my $len = max map length, map @$_, values %fasta;
for my $key ( @keys ) {
print $key, "\n";
my $seq = join '', @{ $fasta{$key} };
print "$_\n" for $seq =~ /.{1,$len}/g;
}
给定原始数据,输出与上述解决方案相同。我用这个作为输入
>seq1
ABCDEFGHI
JKLMNOPQRST
UVWXYZ
>seq2
HIJKLMN
OPQRSTU
VWXY
有了这个结果。所有行都换行到 11 个字符 - 原始数据中最长 JKLMNOPQRST
行的长度
>REVseq1
ZYXWVUTSRQP
ONMLKJIHGFE
DCBA
>REVseq2
YXWVUTSRQPO
NMLKJIH
我不知道这是否仅适用于使用玩具数据集的 class 或大小可能达到千兆字节的实际研究 FASTA。如果是后者,那么不像您的程序和 Borodin 那样将整个数据集保存在内存中是有意义的,而是一次读取一个序列,将其打印出来并忘记它。以下代码执行此操作,还处理 may have asterisks as sequence-end markers 的 FASTA 文件,只要它们以 >
开头,而不是 ;
.
#!/usr/bin/perl
use strict;
use warnings;
my $COL_WIDTH = 80;
my $sequence = '';
my $seq_label;
sub print_reverse {
my $seq_label = shift;
my $sequence = reverse shift;
return unless $sequence;
print "$seq_label\n";
for(my $i=0; $i<length($sequence); $i += $COL_WIDTH) {
print substr($sequence, $i, $COL_WIDTH), "\n";
}
}
while(my $line = <>) {
chomp $line;
if($line =~ s/^>/>REV/) {
print_reverse($seq_label, $sequence);
$seq_label = $line;
$sequence = '';
next;
}
$line = substr($line, 0, -1) if substr($line, -1) eq '*';
$sequence .= $line;
}
print_reverse($seq_label, $sequence);
我不知道这是否只是 Stawberry Perl 的一个怪癖,但我似乎无法做到 运行。我只需要做一个 fasta 并反转其中的每个序列。
-问题-
我有一个 multifasta 文件:
>seq1
ABCDEFG
>seq2
HIJKLMN
预期输出为:
>REVseq1
GFEDCBA
>REVseq2
NMLKJIH
脚本在这里:
$NUM_COL = 80; ## set the column width of output file
$infile = shift; ## grab input sequence file name from command line
$outfile = "test1.txt"; ## name output file, prepend with “REV”
open (my $IN, $infile);
open (my $OUT, '>', $outfile);
$/ = undef; ## allow entire input sequence file to be read into memory
my $text = <$IN>; ## read input sequence file into memory
print $text; ## output sequence file into new decoy sequence file
my @proteins = split (/>/, $text); ## put all input sequences into an array
for my $protein (@proteins) { ## evaluate each input sequence individually
$protein =~ s/(^.*)\n//m; ## match and remove the first descriptive line of
## the FATA-formatted protein
my $name = ; ## remember the name of the input sequence
print $OUT ">REV$name\n"; ## prepend with #REV#; a # will help make the
## protein stand out in a list
$protein =~ s/\n//gm; ## remove newline characters from sequence
$protein = reverse($protein); ## reverse the sequence
while (length ($protein) > $NUM_C0L) { ## loop to print sequence with set number of cols
$protein =~ s/(.{$NUM_C0L})//;
my $line = ;
print $OUT "$line\n";
}
print $OUT "$protein\n"; ## print last portion of reversed protein
}
close ($IN);
close ($OUT);
print "done\n";
这会按照你的要求去做
它从 FASTA 文件中构建一个散列 %fasta
,保持数组 @keys
以保持序列的顺序,然后打印出散列的每个元素
序列的每一行在添加到哈希之前使用reverse
反转,使用unshift
以相反的顺序添加序列的行
程序期望输入文件作为命令行参数,并将结果打印到STDOUT,可在命令行重定向
use strict;
use warnings 'all';
my (%fasta, @keys);
{
my $key;
while ( <> ) {
chomp;
if ( s/^>\K/REV/ ) {
$key = $_;
push @keys, $key;
}
elsif ( $key ) {
unshift @{ $fasta{$key} }, scalar reverse;
}
}
}
for my $key ( @keys ) {
print $key, "\n";
print "$_\n" for @{ $fasta{$key} };
}
输出
>REVseq1
GFEDCBA
>REVseq2
NMLKJIH
更新
如果您更喜欢重新包装序列以便短行在末尾,那么您只需要重写转储哈希的代码
该方案以原文件中最长行的长度为限,将反转后的序列重新包装成相同的长度。更清楚的是,指定一个显式长度而不是计算它会很简单
您需要在程序顶部添加 use List::Util 'max'
my $len = max map length, map @$_, values %fasta;
for my $key ( @keys ) {
print $key, "\n";
my $seq = join '', @{ $fasta{$key} };
print "$_\n" for $seq =~ /.{1,$len}/g;
}
给定原始数据,输出与上述解决方案相同。我用这个作为输入
>seq1
ABCDEFGHI
JKLMNOPQRST
UVWXYZ
>seq2
HIJKLMN
OPQRSTU
VWXY
有了这个结果。所有行都换行到 11 个字符 - 原始数据中最长 JKLMNOPQRST
行的长度
>REVseq1
ZYXWVUTSRQP
ONMLKJIHGFE
DCBA
>REVseq2
YXWVUTSRQPO
NMLKJIH
我不知道这是否仅适用于使用玩具数据集的 class 或大小可能达到千兆字节的实际研究 FASTA。如果是后者,那么不像您的程序和 Borodin 那样将整个数据集保存在内存中是有意义的,而是一次读取一个序列,将其打印出来并忘记它。以下代码执行此操作,还处理 may have asterisks as sequence-end markers 的 FASTA 文件,只要它们以 >
开头,而不是 ;
.
#!/usr/bin/perl
use strict;
use warnings;
my $COL_WIDTH = 80;
my $sequence = '';
my $seq_label;
sub print_reverse {
my $seq_label = shift;
my $sequence = reverse shift;
return unless $sequence;
print "$seq_label\n";
for(my $i=0; $i<length($sequence); $i += $COL_WIDTH) {
print substr($sequence, $i, $COL_WIDTH), "\n";
}
}
while(my $line = <>) {
chomp $line;
if($line =~ s/^>/>REV/) {
print_reverse($seq_label, $sequence);
$seq_label = $line;
$sequence = '';
next;
}
$line = substr($line, 0, -1) if substr($line, -1) eq '*';
$sequence .= $line;
}
print_reverse($seq_label, $sequence);