1

那是我的第一篇文章。我想写一个小脚本来计算一行中的多个唯一重复。文本是 DNA 序列，在此处输入链接描述，因此文本将是四个字母的组合：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

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

4

1 回答 1

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

这个怎么运作：

  • 由于-nshebang 中的标志，脚本针对来自的每一行执行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)是所有匹配的数组的大小，我们将其添加到计数中
于 2013-09-28T11:10:21.697 回答