使用Perl解析FASTA格式的序列

如何解决使用Perl解析FASTA格式的序列

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

>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]/[$1]/s' file

但它不起作用。

解决方法

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

$ # 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所需条件

使用Bio::SeqIO中的BioPerl模块来处理FASTA格式的序列。例如，这会在ID后面附加“ .new”，并将序列更改为小写：

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,);
# 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

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

conda create --name bioperl perl-bioperl