在Perl中需要峰值信号检测方面的帮助

Question

大家好我有一些酵母菌群平板图像的强度值。我需要能够从强度值中找到峰值。下面是一个示例图像，显示了绘图时值的外观。

的一些价值观

 
5.7 
5.3 
8.2 
16.5 
34.2 
58.8 
**75.4** 
75 
65.9 
62.6 
58.6 
66.4 
71.4 
53.5 
40.5 
26.8 
14.2 
8.6 
5.9 
7.7 
14.9 
30.5 
49.9 
69.1 
**75.3** 
69.8 
58.8 
57.2 
56.3 
67.1 
69 
45.1 
27.6 
13.4 
8 
5 

这些值显示在75.4和75.3两个峰的例子，你可以看到，值增加后减小。这种变化并不总是一样的。

图强度值

http://lh4.ggpht.com/_aEDyS6ECO8s/THKTLgDPhaI/AAAAAAAAAio/HQW7Ut-HBhA/s400/peaks.png从 research

其中之一，我想这样做是为了保存每个组，即山中的哈希，然后寻找一个组中的最大价值的东西。其中一个如果我看到的问题是如何确定每个组的边界。

这里是我到目前为止的代码的链接： http://paste-it.net/public/y485822/

这里是一个链接到一个完整的数据集： http://paste-it.net/public/ub121b4/

我写我的代码在Perl。任何帮助将不胜感激。谢谢

Answer 1

你需要决定你想要的峰值的地方。这里的方法在数据的广泛区域内找到峰值和谷值。

use strict; 
use warnings; 

my @data = (
    5.7, 5.3, 8.2, 16.5, 34.2, 58.8, 75.4, 75, 65.9, 62.6, 
    58.6, 66.4, 71.4, 53.5, 40.5, 26.8, 14.2, 8.6, 5.9, 7.7, 
    14.9, 30.5, 49.9, 69.1, 75.3, 69.8, 58.8, 57.2, 56.3, 67.1, 
    69, 45.1, 27.6, 13.4, 8, 5, 
); 

# Determine mean. Or use Statistics::Descriptive. 
my $sum; 
$sum += $_ for @data; 
my $mean = $sum/@data; 

# Make a pass over the data to find contiguous runs of values 
# that are either less than or greater than the mean. Also 
# keep track of the mins and maxes within those groups. 
my $group = -1; 
my $gt_mean_prev = ''; 
my @mins_maxs; 
my $i = -1; 

for my $d (@data){ 
    $i ++; 
    my $gt_mean = $d > $mean ? 1 : 0; 

    unless ($gt_mean eq $gt_mean_prev){ 
     $gt_mean_prev = $gt_mean; 
     $group ++; 
     $mins_maxs[$group] = $d; 
    } 

    if ($gt_mean){ 
     $mins_maxs[$group] = $d if $d > $mins_maxs[$group]; 
    } 
    else { 
     $mins_maxs[$group] = $d if $d < $mins_maxs[$group]; 
    } 

    $d = { 
     i  => $i, 
     val  => $d, 
     group => $group, 
     gt_mean => $gt_mean, 
    }; 
} 

# A fun picture. 
for my $d (@data){ 
    printf 
     "%6.1f %2d %1s %1d %3s %s\n", 
     $d->{val}, 
     $d->{i}, 
     $d->{gt_mean} ? '+' : '-', 
     $d->{group}, 
     $d->{val} == $mins_maxs[$d->{group}] ? '==>' : '', 
     '.' x ($d->{val}/2), 
    ; 

}

输出：

5.7 0 - 0  .. 
    5.3 1 - 0 ==> .. 
    8.2 2 - 0  .... 
    16.5 3 - 0  ........ 
    34.2 4 - 0  ................. 
    58.8 5 + 1  ............................. 
    75.4 6 + 1 ==> ..................................... 
    75.0 7 + 1  ..................................... 
    65.9 8 + 1  ................................ 
    62.6 9 + 1  ............................... 
    58.6 10 + 1  ............................. 
    66.4 11 + 1  ................................. 
    71.4 12 + 1  ................................... 
    53.5 13 + 1  .......................... 
    40.5 14 - 2  .................... 
    26.8 15 - 2  ............. 
    14.2 16 - 2  ....... 
    8.6 17 - 2  .... 
    5.9 18 - 2 ==> .. 
    7.7 19 - 2  ... 
    14.9 20 - 2  ....... 
    30.5 21 - 2  ............... 
    49.9 22 + 3  ........................ 
    69.1 23 + 3  .................................. 
    75.3 24 + 3 ==> ..................................... 
    69.8 25 + 3  .................................. 
    58.8 26 + 3  ............................. 
    57.2 27 + 3  ............................ 
    56.3 28 + 3  ............................ 
    67.1 29 + 3  ................................. 
    69.0 30 + 3  .................................. 
    45.1 31 + 3  ...................... 
    27.6 32 - 4  ............. 
    13.4 33 - 4  ...... 
    8.0 34 - 4  .... 
    5.0 35 - 4 ==> .. 

Answer 2

my @data = ...; 

# filter out sequential duplicate values 
my @orig_index = 0; 
my @deduped = $data[0]; 
for my $index (1..$#data) { 
    if ($data[$index] != $data[$index-1]) { 
     push @deduped, $data[$index]; 
     push @orig_index, $index; 
    } 
} 

# add a sentinel (works for both ends) 
push @deduped, -9**9**9; 

my @local_maxima_indexes; 
for my $index (0..$#deduped-1) { 
    if ($deduped[$index] > $deduped[$index-1] && $deduped[$index] > $deduped[$index+1]) { 
     push @local_maxima_indexes, $orig_index[$index]; 
    } 
} 

注意，这考虑了第一值的局部最大值，也是值71.4和69.我不知道你怎么样正在区分你想要包括哪些。

Answer 3

你有控制数据集吗？如果是这样，我建议使用酵母强度和对照图像之间简单的对数比例来标准化数据。

然后，您可以使用perl port of ChiPOTle抢显著峰，这听起来方式比搜索局部/全局最大值更稳健等

的Chipotle“是用来分析芯片到芯片微阵列数据的峰值寻找算法“，但是我在其他许多应用程序中成功地使用了它（例如ChIP-seq，它确实比您的情况更接近它的原始目的）。

生成的日志（酵母/对照）负值将用于建立显着性估计的高斯背景模型。该算法然后使用错误发现率进行多个测试校正。

Here's the original paper。