标题匹配成功后将 FASTA 序列打印到输出文件的同一行答案

【问题标题】：Print FASTA sequence after successful header match onto the same line in output file标题匹配成功后将 FASTA 序列打印到输出文件的同一行
【发布时间】：2014-10-22 02:33:51
【问题描述】：

在上一个问题之后，我有一些代码几乎可以完成我打算做的事情，但并不完全。

我正在尝试将 FILE1（3750/126477 等）中的每个位置与 FILE2（517-1878,2156-3289 等）中的每个范围进行比较。如果它在一个范围内，我想将位置、范围、方向以及 FASTA 序列从下一行打印到输出文件的同一行。目前，如果多个位置位于同一范围内，那么它将所有位置分组到一个块中，然后仅在最后一行包含序列，当我希望每个匹配都包含相关基因序列时。

我的 FILE1 示例数据：

7065_8#10   3750    -   t
7065_8#10   126477  -   c
7065_8#10   1200    +   T
7065_8#10   3800    -   t

我的 FILE2 示例数据：

>SAEMRSA15_00010 dnaA_chromosomal_replication_initiator_protein_DnaA 517  1878 forward
ATGTCGGAAAAAGAAATTTGGGAAAAAGTGCTTGAAATTGCTCAAGAAAAATTATCAGCTGTAAGTTACTCAACTTTCCTAAA
>SAEMRSA15_00020 dnaN_DNA_polymerase_III,_beta_chain 2156  3289 forward
ATGATGGAATTCACTATTAAAAGAGATTATTTTATTACACAATTAAATGACACATTAAAAGCTATTTCACCAAGAACAACA
>SAEMRSA15_00030 conserved_hypothetical_protein 3670  3915 forward
GTGATTATTTTGGTTCAAGAAGTTGTAGTAGAAGGAGACATTAATTTAGGTCAATTTCTAAAAACAGAAGGGATTATTGAATCTGGTGGTCAAG

我的代码：

#!/usr/bin/perl 

use strict;
use warnings;
use autodie;

my $outputfile = "/Users/edwardtickle/Documents/CC22CDS.txt"; 

open FILE1, "/Users/edwardtickle/Documents/CC22indels.tab";

open FILE2, "/Users/edwardtickle/Documents/CC22_CDS_rmmge.aln";

open (OUTPUTFILE, ">$outputfile");
my @file1list=();

while (<FILE1>) {
if (/^\S+\s+(\d+)/) {
push @file1list, $1;
}
}

my $nextline = 0;
close FILE1;

while ( my $line = <FILE2> ) {
if ($nextline) {
    print OUTPUTFILE "$line\n";
    $nextline = '';
}
elsif ($line =~ /^>(\S+)\s+\S+\s+(\d+)\s+(\d+)\s+(\S+)/) {
    my $cds1 = $1;
    my $cds2 = $2;
    my $cds3 = $3;
    my $cds4 = $4;

    for my $cc22 (@file1list) {
        if ( $cc22 > $cds2 && $cc22 < $cds3 ) {
            $nextline++;
            print OUTPUTFILE "$cc22 $cds2 $cds3 $cds4\n";
        }
    }
    }
  }

close FILE2;

我的结果：

1200 517 1878 forward
ATGTCGGAAAAAGAAATTTGGGAAAAAGTGCTTGAAATTGCTCAAGAAAAATTATCAGCTGTAAGTTACTCAACTTTCCTAAA

3750 3670 3915 forward
3800 3670 3915 forward
 GTGATTATTTTGGTTCAAGAAGTTGTAGTAGAAGGAGACATTAATTTAGGTCAATTTCTAAAAACAGAAGGGATTATTGAATCTGGTGGTCAAG

我想要的结果：

1200 517 1878 forward
ATGTCGGAAAAAGAAATTTGGGAAAAAGTGCTTGAAATTGCTCAAGAAAAATTATCAGCTGTAAGTTACT

3750 3670 3915 forward    GTGATTATTTTGGTTCAAGAAGTTGTAGTAGAAGGAGACATTAATTTAGGTCAATTTCTAAAAACAGAAGGGATTATTGA
3800 3670 3915 forward
GTGATTATTTTGGTTCAAGAAGTTGTAGTAGAAGGAGACATTAATTTAGGTCAATTTCTAAAAACAGAAGGGATTATTGAATCTGGTGGTCAAG

我相信这是因为代码的第一部分在第二个 if 规则之前永远不会匹配，但我不知道如何在保持其功能的同时更改顺序。

或者，如果标题匹配包含字母 ATCG（它显然总是会），是否有一种方法可以在标题匹配后打印下一行。这会让我觉得这是一种更有效的方式，但我又不知道从哪里开始。

感谢您的帮助！

【问题讨论】：

您可以使用$nextline = <FILE2>; 拉入下一行，并在解析文件2 时省略if 循环的第一部分——删除$nextline++ 并将打印语句更改为print OUTPUTFILE "$cc22 $cds2 $cds3 $cds4 $nextline";。
再次感谢您的帮助。我试过了，但我无法让它工作，请问你能告诉我你建议更改的代码吗？我删除了第一个 if 循环，将 'elsif' 更改为 'if'，更改了 nextline = FILE2 并更改了打印语句。

标签： perl

【解决方案1】：

为了在不过多更改现有代码的情况下获得您想要的结果，您可以在处理标题行时获取序列行：

while ( my $line = <FILE2> ) {
    if ($line =~ /^>(\S+)\s+\S+\s+(\d+)\s+(\d+)\s+(\S+)/) {
        my $cds1 = $1;
        my $cds2 = $2;
        my $cds3 = $3;
        my $cds4 = $4;
        # fetch the next line from the file -- i.e. the sequence
        $nextline = <FILE2>;

        for my $cc22 (@file1list) {
            if ( $cc22 > $cds2 && $cc22 < $cds3 ) {
                print "$cc22 $cds2 $cds3 $cds4 $nextline";
            }
        }
    }
}

【讨论】：

效果很好，并且文件的格式与我的代码的下一阶段工作所需的完全相同，谢谢！我实际上也理解您所做的更改，这是一个巨大的帮助。

【解决方案2】：

您可以使用内部循环打印同一范围内的所有匹配项。

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

open my $IND, '<', 'file1' or die $!;
my @pos;
while (<$IND>) {
    push @pos, (split)[1];
}

@pos = sort { $a <=> $b } @pos;

open my $FST, '<', 'file2' or die $!;
while (<$FST>) {
    next unless /^>/;
    my ($from, $to, $direction) = (split)[2 .. 4];
    shift @pos while $pos[0] < $from;
    next if $pos[0] > $to;

    my $nextline = <$FST>;
    while ($pos[0] <= $to) {
        print "$pos[0] $from $to $direction\n";
        print $nextline;
        shift @pos;
    }
}

【讨论】：