计算一行中的多个唯一字符串

Question

这是我的第一篇文章。我想写一个小脚本来计算一行中的多个唯一重复。文本是一个 DNA 序列enter link description here，因此文本将是四个字母的组合：A、T、G 和 C。 如果一个字符串出现两次，就会被计算两次，以此类推。

我要查找的唯一字符串是三个 AG、GA、CT 或 TC 的重复，分别是 (AG)3、(GA)3、(CT)3 和 (TC)3。我不希望程序计算四次或更多的重复次数。

要计数的字符串：

AGAGAG
GAGAGA
CTCTCT
TCTCTC

示例输入文件（由制表符分隔的两列）：

Sequence_1    AGAGAG                   
Sequence_2    AGAGAGT                  
Sequence_3    AGAGAGAG                 
Sequence_4    AGAGAT                   
Sequence_5    AGAGAGAGAGAGAGAGAGT      
Sequence_6    AGAGAGTAGAGAG 
Sequence_7    CTCTCTCTCTC  
Sequence_8    TAGAGAGAT                
Sequence_9    TAAGAGAGAAG              

期望的输出：

Sequence_1    AGAGAG                   1
Sequence_2    AGAGAGT                  1
Sequence_3    AGAGAGAG                 0
Sequence_4    AGAGAT                   0
Sequence_5    AGAGAGAGAGAGAGAGAG       0
Sequence_6    AGAGAGTAGAGAG            2
Sequence_7    CTCTCTCTCTCAAGAGAG       1 
Sequence_8    TAGAGAGAT                1
Sequence_9    TAAGAGAGAAG              1

我有一个用awk写的小单行，但我认为匹配字符串时并不具体：

awk '{if($1 ~ /AGAGAG/)x++; if($1 ~ /TCTCTC/)x++;if($1 ~ /GAGAGA/)x++;if($1 ~ /CTCTCT/)x++;print x;x=0}' inputfile.tab

非常感谢您的帮助。一切顺利，贝尔纳多

Answer 1

我认为您的描述以及示例输入和输出中存在一些不一致之处。所以这个脚本可能并不完美，但我希望它足够接近你可以弄清楚其余部分：

#!/usr/bin/perl -n

my ($seq, $dna) = split(/\s+/);
my @strings = qw/AG GA CT TC/;
my $count = 0;
foreach my $s (@strings) {
    my ($b, $e) = split(//, $s);
    @matches = $dna =~ m/(?<!$e)($s){3}(?!$b)/g;
    $count += scalar(@matches);
}
print join("\t", $seq, sprintf("%-20s", $dna), $count), "\n";

你可以使用它：

./script.pl < sample.txt

输入：

Sequence_1    AGAGAG
Sequence_2    AGAGAGT
Sequence_3    AGAGAGAG
Sequence_4    AGAGAT
Sequence_5    AGAGAGAGAGAGAGAGAGT
Sequence_6    AGAGAGTAGAGAG
Sequence_7    CTCTCTCTCTCAAGAGAG

它给出：

Sequence_1    AGAGAG                1
Sequence_2    AGAGAGT               1
Sequence_3    AGAGAGAG              0
Sequence_4    AGAGAT                0
Sequence_5    AGAGAGAGAGAGAGAGAGT   0
Sequence_6    AGAGAGTAGAGAG         2
Sequence_7    CTCTCTCTCTCAAGAGAG    1

它是如何工作的：

• 感谢 shebang 中的-n 标志，脚本针对来自stdin 的每一行执行
• @strings 是我们感兴趣的字符串列表
• 对于@strings 中的每个项目，我们计算匹配项
  • $s 取值 AG, GA, CT, TC
  • 表达式(?<!$s)($s){3}(?!$s) 匹配3 个连续的$s，后面没有$s，前面也没有$s
  • 表达式 (?<!$e)($s){3}(?!$b) 匹配 3 个连续的 $s，其后面没有 $s 的第一个字符，并且前面没有 $s 的第二个字符
  • $x =~ m///g 操作返回一个包含所有匹配项的数组
  • scalar(@matches) 是所有匹配数组的大小，我们将其添加到计数中