Perl 程序将完整的 .fasta 文件序列打印到文件中，但试图实现与基因相关的特定核苷酸计数

Question

到目前为止，我已经设法对 Perl 有了更多的了解，这让我松了一口气，我要感谢你们。我目前仍在研究另一个方面，我需要读取 .fasta 文件并找到所有 G 和 C 核苷酸，然后创建一个制表符分隔的文件。

这些是我过去几天的帖子，按时间顺序排列：

1. How do I average column values from a tab-separated data... （已解决）
2. Why do I see no computed results in my output file? （已解决）
3. Using a .fasta file to compute relative content of sequences
4. Reading .fasta sequences to extract nucleotide data, and then... （在此之前发布）

最后一个查询仍在处理中，但我已经取得了一些进展。

一些背景，.fasta 文件的内容如下：

>label
sequence
>label
sequence
>label
sequence

我不确定如何打开 .fasta 文件，所以我不确定哪些标签适用于哪个，但我知道基因应该标记为 gag、pol 或 env。我是否需要打开 .fasta 文件才能知道我在做什么，或者我可以通过上述格式“盲目”地做吗？

反正我现在的代码如下：

#!/usr/bin/perl -w
# This script reads several sequences and computes the relative content of G+C of each sequence.

use strict; 

my $infile = "Lab1_seq.fasta";                               # This is the file path
open INFILE, $infile or die "Can't open $infile: $!";        # This opens file, but if file isn't there it mentions this will not open
my $outfile = "Lab1_SeqOutput.txt";             # This is the file's output
open OUTFILE, ">$outfile" or die "Cannot open $outfile: $!"; # This opens the output file, otherwise it mentions this will not open

my $sequence = ();  # This sequence variable stores the sequences from the .fasta file
my $GC = 0;         # This variable checks for G + C content

my $line;                             # This reads the input file one-line-at-a-time

while ($line = <INFILE>) {
    chomp $line;                      # This removes "\n" at the end of each line (this is invisible)

    if($line =~ /^\s*$/) {         # This finds lines with whitespaces from the beginning to the ending of the sequence. Removes blank line.
        next;

    } elsif($line =~ qr(^\s*\#/)) {        # This finds lines with spaces before the hash character. Removes .fasta comment
        next; 
    } elsif($line =~ /^>/) {           # This finds lines with the '>' symbol at beginning of label. Removes .fasta label
        next;
    } else {
        $sequence = $line;
    }

    $sequence =~ s/\s//g;               # Whitespace characters are removed
    print OUTFILE $sequence;
}

代码现在将整个序列打印到文本文件中，没有空格。唯一的问题是，我不知道序列从哪里开始或结束，所以我不确定哪些序列适用于每个基因。虽然停止/开始密码子应该给我一个指示。考虑到这一点，我将如何修改/添加到代码中以计算序列中 G+C 的数量，然后将它们打印到一个制表符分隔的文件中，其中包含与它们各自的 G/C 内容相关的基因名称?

我期待听到有人可以提供一些指导，与上面发布的代码类似，关于如何找到 G/C，然后将各个计数制成表格。

Answer 1

以下链接可能会有所帮助。已经写了很多代码，Bio::SeqIO 和Bio::Seq 似乎经常被讨论。 BioPerl 有一个网站，但我不熟悉。那里有代码示例和其他信息。 FAQ 也很有帮助。

这是来自 Bio::SeqIO 文档的示例。

use Bio::SeqIO;

$in  = Bio::SeqIO->new(-file => "inputfilename" ,
                       -format => 'Fasta');
$out = Bio::SeqIO->new(-file => ">outputfilename" ,
                       -format => 'EMBL');

while ( my $seq = $in->next_seq() ) {
        $out->write_seq($seq);
}

Answer 2

我实际上自己使用FASTA 文件。所以，我感受到了你的痛苦。

回答您关于标签对每个序列的适用性的重复问题：如果文件格式正确，则序列信息之前的每个标签都应该用于后面的序列。因此，您应该按如下方式从头到尾解析文件：

>label1
sequence1
>label2
sequence2
>label3
sequence3
...

... 其中每个标签表示要遵循的新序列信息。您还需要忽略以分号 (;) 开头的行，因为这些行也表示旧版注释字段。

否则，您似乎在重排文件时正确地删除了空格。我建议使用换行符保持标签字段完整，因此您的输出文件看起来像上面提到的格式，删除了 cmets 和空格。

一旦你有了这个，这很简单，遍历重排文件，抓取你需要的序列片段，并在遇到新标签时重新启动计数器。