使用 Perl 解析 FASTA 格式的序列 [关闭]答案

【问题标题】：Parse sequences in FASTA format with Perl [closed]使用 Perl 解析 FASTA 格式的序列 [关闭]
【发布时间】：2021-01-02 01:55:57
【问题描述】：

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

>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

但它不起作用。

【问题讨论】：

预期输出是什么？
试试perl -0777 -pe 's/^(?!>).+(?:\R[A-Z].*)+/[$&]/gm' file
@Wiktor Stribiżew 输出可以是任何东西，我只想检查是否匹配多行记录而不是单行记录。感谢您的代码，但它似乎不起作用。
见ideone.com/BqoHnb 和regex101.com/r/7CdBAk/1，什么不起作用？
有些模块可以读取这些 dna-string 文件。找到这样一个模块并使用它可能是值得的。

标签： regex shell perl bioinformatics fasta

【解决方案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 所需条件

【讨论】：

【解决方案2】：

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

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

【讨论】：