使用 Perl 解析 FASTA 格式的序列

Question

我有以下文件，其中一些记录（每条记录由以“>”开头的行分隔）在一行中，而其他记录在多行中：

>NODE_1506886_length_92_cov_1.000000
CATGCGAACTTCCGCAAGGACATCGTCATCCATGGCTTCAGACATGCCATGCGTGACCGT
CTGAGAGCCGTTAGCTGCCCATCAGAGATGAT
>NODE_1506887_length_92_cov_1.000000
TATATTTTAATTCATTATTTGGTACATCAAACGATGAGGCTAAATATATAATAACTGGAA
ATAATCGAGACGTAAATCGATACCGATGGAAC
>k119_811274 flag=1 multi=3.0000 len=313
TAAAAAATTACTACCTTGCAGATAAGTCGCTCAGCCCATGGTTGTCAGGTTTATCTGCAATGGCAACAAATAACAGTGGCTACATGTTTATTGGTCTTATTGGATTCACCTATCTCAATGGCCTTTCCAGCATATGGCTAATGATTGGTTGGATATTAGGGGATTACCTAATTACAAAAAAATTATTTCCAA
>k119_405638 flag=1 multi=7.0000 len=562
CACAGTTAACATAAATTACCACGAACGACGCACGCTCAAATCCGAAGCAAAGTCAGCAATTGCTAAATCCTGTTCAAAAAAAATTCCAAACGGTAGTTCGATATTTATCAATATCGGGACATCGACTGAAGCCGTCGCTCAGGAATTAATGCAGCACAGTAACTTAATGGTTGTGACCAACAACATAAATGTTGCCAATATTTTATCGCCCAATGAGAATTGTGAGATTCTTTTAACTGGTGGTCAACTTAGACGTTCTGACGGAGGTCTCATCGGTAATTTGGCAGC

我只想匹配多行记录，而不是单行记录。我尝试了以下

perl -0777ne 's/(^[^>].+)\s[A-Z]/[]/s' file

但是没用。

Answer 1

您可以使用 > 作为记录分隔符，而不必将整个文件作为单个字符串读取

$ # to get records with more than one ACGT lines
$ perl -F'\n' -lane 'BEGIN{$/=">"; $\=""} print "$/$_" if $#F>1' ip.txt
>NODE_1506886_length_92_cov_1.000000
CATGCGAACTTCCGCAAGGACATCGTCATCCATGGCTTCAGACATGCCATGCGTGACCGT
CTGAGAGCCGTTAGCTGCCCATCAGAGATGAT
>NODE_1506887_length_92_cov_1.000000
TATATTTTAATTCATTATTTGGTACATCAAACGATGAGGCTAAATATATAATAACTGGAA
ATAATCGAGACGTAAATCGATACCGATGGAAC

$ # to get records with exactly one ACGT lines
$ perl -F'\n' -lane 'BEGIN{$/=">"; $\=""} print "$/$_" if $#F==1' ip.txt
>k119_811274 flag=1 multi=3.0000 len=313
TAAAAAATTACTACCTTGCAGATAAGTCGCTCAGCCCATGGTTGTCAGGTTTATCTGCAATGGCAACAAATAACAGTGGCTACATGTTTATTGGTCTTATTGGATTCACCTATCTCAATGGCCTTTCCAGCATATGGCTAATGATTGGTTGGATATTAGGGGATTACCTAATTACAAAAAAATTATTTCCAA
>k119_405638 flag=1 multi=7.0000 len=562
CACAGTTAACATAAATTACCACGAACGACGCACGCTCAAATCCGAAGCAAAGTCAGCAATTGCTAAATCCTGTTCAAAAAAAATTCCAAACGGTAGTTCGATATTTATCAATATCGGGACATCGACTGAAGCCGTCGCTCAGGAATTAATGCAGCACAGTAACTTAATGGTTGTGACCAACAACATAAATGTTGCCAATATTTTATCGCCCAATGAGAATTGTGAGATTCTTTTAACTGGTGGTCAACTTAGACGTTCTGACGGAGGTCTCATCGGTAATTTGGCAGC

-F'\n' 使用换行符作为字段分隔符，便于统计每条记录的行数
$/=">"; $\="" 将 > 设置为输入记录分隔符，将空字符串设置为输出记录分隔符
$#F>1或$#F==1所需条件

Answer 2

使用Bio::SeqIO module from BioPerl处理FASTA格式的序列。例如，这会将“.new”附加到 id，并将序列更改为小写：

cat in.fa | \
    perl -MBio::SeqIO -e '
my $in_seq = Bio::SeqIO->new( -fh => \*STDIN, -format => "fasta", );
my $out_seq = Bio::SeqIO->new( -fh => \*STDOUT, -format => "fasta", );
# Prevent sequence lines wrapping in bioperl by using arbitrary large width:
$out_seq->width(1e9);
while ( my $seq = $in_seq->next_seq ) {
    my $id = $seq->id;
    $seq->id( "$id.new" );
    $seq->seq( lc $seq->seq );
    $out_seq->write_seq($seq);
}
' > out.fa

head in.fa out.fa
==> in.fa <==
>seq1
ACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTG
>seq2
ACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTG
>seq3
ACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTGACTGCTG

==> out.fa <==
>seq1.new
actgctgactgctgactgctgactgctgactgctgactgctgactgctgactgctgactgctgactgctg
>seq2.new
actgctgactgctgactgctgactgctgactgctgactgctgactgctgactgctgactgctgactgctg
>seq3.new
actgctgactgctgactgctgactgctgactgctgactgctgactgctgactgctgactgctgactgctg

您可以安装 BioPerl，例如，使用 conda:

conda create --name bioperl perl-bioperl

使用 Perl 解析 FASTA 格式的序列

Parse sequences in FASTA format with Perl

regex

shell

perl

bioinformatics

fasta