如何保存匹配和不匹配的 FASTA 序列答案

【问题标题】：How to save matched and nonmatched FASTA sequences如何保存匹配和不匹配的 FASTA 序列
【发布时间】：2017-08-31 12:12:04
【问题描述】：

我有一个这样的 FASTA 序列文件

 >seq002
 ATGGTAAATGGTTTCTCAAATTGTGCACTGACAGACAAACCCCT
 >seq0009
 ATGGCGTCAAAGGTGATGCCGTCAGCGTCAACAACTAA
 >seq0001
 ATGGGAAATAGTGAGGACGGGAAATCTTTAG
 >seq0003
 ATGGGATCTTACTTGAACTTCAAGAATTGA
>seq00005
GCTAATTTTGAGGTTTACCCAGATAGCTG

我正在尝试提取以 ATG 开头并以 TAG/TGA/TAA 结尾的序列。我将此代码用于我的目的

#!/usr/bin/perl -w
# This script reads several sequences and print the sequence which don't strat with ATH and ends with TAG/TGA/TAA

use strict; 

my $infile = "id.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 = "full_length_seq.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 $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 =~m/^ATG[GTAC]+T(GA|AA|AG)$/g) { # This finds lines matching with pattern
        next;

      }

    print OUTFILE $line, "\n";
}

产生这样的结果

>seq002
>seq0009
>seq0001
>seq0003
>seq00005
GCTAATTTTGAGGTTTACCCAGATAGCTG

但我想像这样创建两个不同的文件

>seq002
 ATGGTAAATGGTTTCTCAAATTGTGCACTGACAGACAAACCCCT
 >seq0009
 ATGGCGTCAAAGGTGATGCCGTCAGCGTCAACAACTAA
 >seq0001
 ATGGGAAATAGTGAGGACGGGAAATCTTTAG
 >seq0003
 ATGGGATCTTACTTGAACTTCAAGAATTGA

和

>seq00005
GCTAATTTTGAGGTTTACCCAGATAGCTG

任何提示或帮助将不胜感激。谢谢。

【问题讨论】：

标签： regex perl

【解决方案1】：

所以诀窍是 - 你目前正在分割换行并逐行工作。

但您不必这样做 - 您可以改用 $/，并将其设置为合适的分隔符。

我会建议你这样做，你想要"\n>"，因为那样它就会把你的东西分成几块。

然后，你需要稍微改变你的模式匹配，因为现在每个“块”是两行。

类似这样的：

#!/usr/bin/env perl

use strict;
use warnings;

local $/ = "\n>";

open my $outfile, '>', 'full_length_seq.txt' or die $!;
open my $other_outfile, '>', 'everything_else.txt' or die $!;


while ( <DATA> ) { 
    chomp;
    s/^>//g; #remove leading >, because first line doesn't have a linefeed in front. 

    #just for some diagnostics - print what we're currently operating on. 
    print "New chunk:\n";
    print;
    print "\nEnd\n";


    if ( /\nATG[GTAC]+T(GA|AA|AG)$/ ) {
        print "**matches**\n";
        print {$other_outfile} ">",$_,"\n";
    }
    else {
        print {$outfile} ">", $_, "\n";
    }

}

__DATA__
>seq002
ATGGTAAATGGTTTCTCAAATTGTGCACTGACAGACAAACCCCT
>seq0009
ATGGCGTCAAAGGTGATGCCGTCAGCGTCAACAACTAA
>seq0001
ATGGGAAATAGTGAGGACGGGAAATCTTTAG
>seq0003
ATGGGATCTTACTTGAACTTCAAGAATTGA
>seq00005
GCTAATTTTGAGGTTTACCCAGATAGCTG

这给了我们一个文件：

>seq002
ATGGTAAATGGTTTCTCAAATTGTGCACTGACAGACAAACCCCT
>seq00005
GCTAATTTTGAGGTTTACCCAGATAGCTG

另一个是：

>seq0009
ATGGCGTCAAAGGTGATGCCGTCAGCGTCAACAACTAA
>seq0001
ATGGGAAATAGTGAGGACGGGAAATCTTTAG
>seq0003
ATGGGATCTTACTTGAACTTCAAGAATTGA

我使用上面的 __DATA__ 作为说明 - 您可能仍然应该读取输入文件，或者只使用 <> 来读取“STDIN 或命令行上命名的文件”（如 grep/sed 等。 )

另外，我建议使用带有词法文件句柄的 3 个参数 open 作为更好的样式。例如

open ( my $infile, '<', 'id.fasta' ) or die $!;

【讨论】：