【问题标题】:Selectively joining elements of an array into fewer elements of a new array有选择地将数组的元素连接到新数组的较少元素中
【发布时间】:2012-05-19 13:09:33
【问题描述】:

我在处理 .fasta 格式的 DNA 序列数据数组时遇到了一些问题。我特别想做的是将一个包含几千个序列的文件和文件中每个序列的相邻序列数据放到文件中的单行上。 [Fasta 格式是这样的:序列 ID 以 > 开头,之后该行的所有内容都是描述。在下一行中,存在与此 ID 对应的序列。这可以无限期地持续到以 > 开头的下一行,这是文件中下一个序列的 id] 所以,在我的特定文件中,我的大部分序列都在多行上,所以我想做的基本上是删除换行符,但只删除序列数据之间的新行,而不是序列数据和序列 ID 行(以 > 开头)之间的新行。

我这样做是因为我希望能够获得每个序列的序列长度(通过长度,我相信是最简单的方法),然后获得整个文件中所有序列的平均序列长度。

到目前为止,这是我的脚本,它似乎不想工作:

#!/usr/bin/perl -w


##Subroutine
sub get_file_data1 { 
    my($filename) = $_[0];
    my @filedata = ();
    unless( open(GET_FILE_DATA, $filename)) {
    print STDERR "Cannot open file \"$filename\"\n\n";
    exit;
    }
    @filedata = <GET_FILE_DATA>;
    close GET_FILE_DATA;
    return @filedata;
}



##Opening files
my $fsafile = $ARGV[0];
my @filedata = &get_file_data1($fsafile);


##Procedure
my @count;
my @ids;
my $seq;

foreach $seq (@filedata){
        if ($seq =~ /^>/) {push @ids, $seq;
                                 push @count, "\n";
    }
        else {push @count, $seq;
    }
}


foreach my $line (@count) {
    if ($line =~ /^[AGTCagtc]/){
         $line =~ s/^([AGTCagtc]*)\n/$1/;
    }
}

##Make a text file to have a look
open FILE3, "> unbrokenseq.txt" or die "Cannot open output.txt: $!";

foreach (@count)
{
    print FILE3 "$_\n"; # Print each entry in our array to the file
}
close FILE3;


__END__
##Creating array of lengths
my $number;
my @numberarray;
foreach $number (@count) {
                push @numberarray, length($number);
                }
print @numberarray;


__END__
use List::Util qw(sum);

sub mean {
    return sum(@numberarray)/@numberarray;
}

“过程”部分的第二行 foreach 有问题,我似乎无法弄清楚它是什么。请注意,我什至还没有尝试过 END 行之后的代码,因为我似乎无法在过程步骤中获取代码来执行我想要的操作。知道如何获得一个包含完整序列元素的漂亮数组(我选择只从新数组中删除序列 ID 行..)吗?什么时候我可以得到一个长度数组,然后我可以平均?

最后我很遗憾地承认我无法让 Bio::Perl 在我的计算机上运行,​​我已经尝试了几个小时,但这些错误超出了我的能力范围。我将与希望能帮助我解决 Bio::perl 问题的人交谈。但现在我只好在没有它的情况下继续前进。

谢谢!抱歉这篇文章太长了,感谢您的帮助。

安德鲁

【问题讨论】:

    标签: perl sequence bioinformatics fasta bioperl


    【解决方案1】:

    第二个循环的问题在于,您实际上并未更改 @count 中的任何内容,因为 $line 包含 @count 中值的副本。

    但是,如果您只想在第二个循环中删除末尾的换行符,请使用chomp 函数。有了这个,你就不需要你的第二个循环了。 (而且它也会比使用正则表达式更快。)

    # remove newlines for all array elements before doing anything else with it
    chomp @filedata;
    
    # .. or you can do it in your first loop
    foreach $seq (@filedata){
        chomp $seq;
        if ($seq =~ /^>/) {
        ...
    }
    

    附加提示:如果文件很大,使用get_file_data1 将整个文件读入数组可能会很慢。在这种情况下,最好在执行过程中遍历文件:

    open my $FILE_DATA, $filename or die "Cannot open file \"$filename\"\n";
    while (my $line = <$FILE_DATA>) {
        chomp $line;
        # process the record as in your Procedure section
        ...
    }
    close $FILE_DATA;
    

    【讨论】:

    • 非常感谢您的回复,我认为 chomp'ing 这绝对是要走的路。尽管在考虑了更多之后,我认为我实际上在这里遇到的一个大问题是无法选择性地将数组的多个序列元素压缩为一个元素。我正在尝试将与序列 ID 元素(以 > 开头)对应的所有序列数据元素压缩为一个元素。我天真地以为我可以通过删除序列元素之后的新行来做到这一点,但这不起作用..所以我认为我必须完全重建它。
    【解决方案2】:

    您的正则表达式专门捕获 $1 但您正在将 $_ 打印到文件中。结果很可能不是您想要的。

    【讨论】:

    • 不,那部分代码是正确的。 $_@count 的每个元素的别名。
    【解决方案3】:

    小心使用 s/// 中角色组的“*”或“贪婪”修饰符。您通常需要“+”。 '*' 也将匹配不包含任何字符的行。

    带有“g”修饰符的搜索表达式也可以计算字符数。像这样:

    $perl -e '$a="aggaacaat"; $b = $a =~ s/[a]//g; print $b; '
    5
    

    很酷啊!或者,在您的代码中,您可以针对 $1 调用 length()。

    看到您的正则表达式中转义的 '/n' 让我大吃一惊。虽然它工作正常,但常见的“行尾”搜索词是“$”。这更便于携带,并且不会弄乱您的字符数。

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2023-04-03
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2017-03-19
      • 2011-11-03
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多