如何通过知道坐标来提取子字符串答案

【问题标题】：how to extract substrings by knowing the coordinates如何通过知道坐标来提取子字符串
【发布时间】：2013-08-29 15:37:43
【问题描述】：

非常抱歉在几个问题上打扰您，但我需要解决它...

我想从一个包含字符串的文件中提取几个子字符串，方法是使用另一个文件，其中包含我要提取的每个子字符串的开头和结尾。第一个文件是这样的：

>scaffold30     24194
CTTAGCAGCAGCAGCAGCAGTGACTGAAGGAACTGAGAAAAAGAGCGAGCTGAAAGGAAGCATAGCCATTTGGGAGTGCCAGAGAGTTGGGAGG GAGGGAGGGCAGAGATGGAAGAAGAAAGGCAGAAATACAGGGAGATTGAGGATCACCAGGGAG.........
.................

（字符串必须是文件中除第一行以外的所有内容），坐标文件如下：

44801988    44802104
44846151    44846312
45620133    45620274
45640443    45640543
45688249    45688358
45729531    45729658
45843362    45843490
46066894    46066996
46176337    46176464
.....................

我的脚本是这样的：

my $chrom = $ARGV[0];
my $coords_file = $ARGV[1];

#finds  subsequences: fasta files



open INFILE1, $chrom or die "Could not open $chrom: $!";
my $count = 0;

while(<INFILE1>) {
    if ($_ !~ m/^>/) {

    local $/ = undef;
    my $var = <INFILE1>;

    open INFILE, $coords_file or die "Could not open $coords_file: $!";
           my @cline = <INFILE>;
    foreach my $cline (@cline) {
    print "$cline\n";
            my@data = split('\t', $cline);
            my $start = $data[0];
            my $end = $data[1];
            my $offset = $end - $start;
           $count++;
           my $sub = substr ($var, $start, $offset);
           print ">conserved $count\n";
           print "$sub\n";

    }
    close INFILE;
    }
}

当我运行它时，它看起来只进行了一次迭代，并打印出第一个文件的开头。似乎 foreach 循环不起作用。 substr 似乎也不起作用。当我退出打印 cline 以检查循环时，它会打印文件的所有行以及坐标。

如果我变得烦人我很抱歉，但我必须完成它，我有点绝望......

再次感谢您。

【问题讨论】：

你考虑过使用哈希吗？
我假设$chrom 文件中只有 1 个序列 - 对吧？ id 为>scaffold30 24194 .
糟糕，我想在再次阅读您的帖子后 - （字符串必须是文件中除第一行之外的所有内容）

标签： perl

【解决方案1】：

这一行

local $/ = undef;

为整个封闭块更改$/，其中包括您在第二个文件中读取的部分。 $/ 是输入记录分隔符，它本质上定义了“行”是什么（默认为换行符，详见perldoc perlvar）。当您使用<> 从文件句柄中读取时，$/ 用于确定停止读取的位置。例如，以下程序依赖于默认的换行行为，因此只读取到第一个换行符：

my $foo = <DATA>;
say $foo;
# Output:
# 1

__DATA__
1
2
3

而这个程序一直读取到 EOF：

local $/;
my $foo = <DATA>;
say $foo;
# Output:
# 1
# 2
# 3

__DATA__
1
2
3

这意味着你的@cline 数组只有一个元素，它是一个包含整个坐标文件文本的字符串。您可以使用Data::Dumper 看到这一点：

use Data::Dumper;

print Dumper(\@cline);

在您的情况下会输出如下内容：

$VAR1 = [
          '44801988    44802104
44846151    44846312
45620133    45620274
45640443    45640543
45688249    45688358
45729531    45729658
45843362    45843490
46066894    46066996
46176337    46176464
'
        ];

请注意您的数组（在本例中为 arrayref），由[ 和] 描述，它只包含一个元素，它是一个包含换行符的字符串（由单引号描述）。

让我们浏览一下代码的相关部分：

while(<INFILE1>) {
    if ($_ !~ m/^>/) {
        # Enable localized slurp mode. Stays in effect until we leave the 'if'
        local $/ = undef;

        # Read the rest of INFILE1 into $var (from current line to EOF)
        my $var = <INFILE1>;

        open INFILE, $coords_file or die "Could not open $coords_file: $!";

        # In list context, return each block until the $/ character as a
        # separate list element. Since $/ is still undef, this will read
        # everything until EOF into our first list element, resulting in
        # a one-element array
        my @cline = <INFILE>;

        # Since @cline only has one element, the loop only has one iteration
        foreach my $cline (@cline) {

附带说明一下，您的代码可以稍微清理一下。您为文件句柄选择的名称有待改进，无论如何您都应该使用词法文件句柄（以及open 的三参数形式）：

open my $chromosome_fh,  "<", $ARGV[0] or die $!;
open my $coordinates_fh, "<", $ARGV[1] or die $!;

此外，在这种情况下，您不需要嵌套循环，它只会让您的代码更加复杂。首先将染色体文件的相关部分读入一个变量（命名为比var 更有意义的名称）：

# Get rid of the `local $/` statement, we don't need it
my $chromosome;
while (<$chromosome_fh>) {
    next if /^>/;
    $chromosome .= $_;
}

然后读入你的坐标文件：

my @cline = <$coordinates_fh>;

或者，如果您只需要使用坐标文件的内容一次，请使用 while 循环处理每一行：

while (<$coordinates_fh>) {
    # Do something for each line here
}

【讨论】：

问题是我使用这一行是因为我想加载 $var 上的所有文件，以便使用坐标来提取子字符串
@Vasilis 很好，但是如果你想在换行符上拆分，你需要在读取第二个文件之前将$/ 改回\n。
@ThisSuitlsBlackNot 谢谢你。现在我明白了！但是我怎样才能把 $/ 改回正常呢？
@Vasilis 您可以使用local $/ = "\n" 将其更改回同一块中的换行符，但我真的建议您清理您的代码，在这种情况下您甚至不需要弄乱首先是$/。请参阅我答案的最后一部分。

【解决方案2】：

正如“ThisSuitIsBlackNot”所建议的，您的代码可以稍微清理一下。这是一个可能的解决方案，可能是您想要的。

#!/usr/bin/perl
use strict;
use warnings;

my $chrom = $ARGV[0];
my $coords_file = $ARGV[1];

#finds  subsequences: fasta files

open INFILE1, $chrom or die "Could not open $chrom: $!";
my $fasta;

<INFILE1>; # get rid of the first line - '>scaffold30     24194'

while(<INFILE1>) {
    chomp;
    $fasta .= $_;
}
close INFILE1 or die "Could not close '$chrom'. $!";

open INFILE, $coords_file or die "Could not open $coords_file: $!";
my $count = 0;

while(<INFILE>) {
    my ($start, $end) = split;

    # Or, should this be: my $offset = $end - ($start - 1);
    # That would include the start fasta
    my $offset = $end - $start;

    $count++;
    my $sub = substr ($fasta, $start, $offset);
    print ">conserved $count\n";
    print "$sub\n";
}
close INFILE or die "Could not close '$coords_file'. $!";

【讨论】：