☰

如何在 FASTA 文件中找到“模式”并记录坐标和标题

tag-icon • tag-icon bash python shell-script awk perl
如何在 FASTA 文件中找到“模式”并记录坐标和标题

我正在寻找一种解决方案，用于在 FASTA 文件中搜索 17 个碱基对的字符串，即人类参考基因组。为了澄清起见，用简单的语言和资源，我尝试使用 grep 函数来计算“模式”的频率。

样本参考文件示例：

>chr1 CP068277.2 Homo sapiens isolate CHM13 chromosome 1
CaccctaaaccctaaccGTACATAAAATATGAAAcctaaccctaaccctaaccctaaccctaaccctaacccctaaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccct
aaccctaaccctaacccctaaccctaaccctaaccctaaccctaaccctaaccctGTACATAAAATATGAAAaaccctaaccctaaccctaaccctaaccctaacccGTACATAAAATATGAAAtaaccctaaccctaaccctaaccctaaccctaacccaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccct
aaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaa
ccctaaccctaaccctaaccctaGTACATAAAATATGAAAaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaGTACATAAAATATGAAAccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaa

我使用这种方法来计算模式是否存在，

grep -i "GTACATAAAATATGAAA" ref.fa |wc -l

此代码给出的结果是 49678。但这是完整搜索，在参考示例中，我有意添加了 5 个大写字母的延伸。

此外我使用，

grep -o -P -b -n -i "GTACATAAAATATGAAA" ref.fa

结果是：35659119：2888387418：gtacataaaatatgaaa

很可能，它是坐标和延伸本身，没有标题。

我需要的结果如下：

染色体 开始 结尾
chr1 x 位置 y 位置
chr1 x1 位置 y1 位置

所以这里，chr1 是一个示例标题，实际上，我的 fasta 文件将包含 23 条染色体和它自己的 FASTA 行。

我希望我能解释我的疑问，但如果不能，我会尽量详细说明。谢谢。

答案1

[编辑以添加不区分大小写的匹配]

尝试以下代码

#/usr/bin/perl

use strict;
use warnings;


# first commandline parameter is the pattern to match
my $want = shift @ARGV ;

my $chromosome ;

while (<>)
{
    # match lines that start with a ">" and remember the first word
    $chromosome = $1
        if /^>(\w+)/ ;

    while (/($want)/gi)
    {
        # remember the length of the matched string
        my $match = length $1;

        # "@-" stores the offsets to each of the submatches
        # "$-[0]" is the first one, namely the offset to the start that we need
        my $start = $-[0];

        # work out the offset for the end of the match
        my $end = $start + $match - 1 ;

        print "$chromosome\t$start\t$end\n";
    }

}

假设数据ref.fa包含

>chr1 CP068277.2 Homo sapiens isolate CHM13 chromosome 1
CaccctaaaccctaaccGTACATAAAATATGAAAcctaaccctaaccctaaccctaaccctaaccctaacccctaaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccct
aaccctaaccctaacccctaaccctaaccctaaccctaaccctaaccctaaccctGTACATAAAATATGAAAaaccctaaccctaaccctaaccctaaccctaacccGTACATAAAATATGAAAtaaccctaaccctaaccctaaccctaaccctaacccaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccct
aaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaa
ccctaaccctaaccctaaccctaGTACATAAAATATGAAAaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaGTACATAAAATATGAAAccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaaccctaa

运行代码得到

$ perl try.pl GTACATAAAATATGAAA ref.fa
chr1    17      33
chr1    55      71
chr1    107     123
chr1    23      39
chr1    137     153

如果您需要获取序列中的真实偏移量（即通过将所有序列行连接在一起），则可以按如下方式更改脚本

#/usr/bin/perl

use strict;
use warnings;

# enable paragraph mode
local $/ = '';

my $want = shift @ARGV ;

my $chromosome ;

while (<>)
{
    # get the lines in the paragraph
    my @lines = split "\n", $_;

    # first line is the title
    my $title = shift @lines;

    # only want paragraphs where the first lines begins with a ">"
    next
        if ! /^>(\w+)/ ;

    # remember the first word
    $chromosome = $1 ;

    # concatenate the "sequence" lines into a single string with no newlines
    my $full_sequence = join "", @lines;

    while ($full_sequence =~ /($want)/gi)
    {
        # remember the length of the matched string
        my $match = length $1;

        # @- srores the offsets to each of the submatches
        # $-[0] is the first one
        my $start = $-[0];

        # work out the offset for the end of the match
        my $end = $start + $match - 1 ;

        print "$chromosome\t$start\t$end\n";
    }
}

跑步

$ perl try.pl GTACATAAAATATGAAA ref.fa
chr1    17      33
chr1    232     248
chr1    284     300
chr1    554     570
chr1    668     684

相关内容