根据地图估算 class 个值

Impute class values based on map

我想根据已知标记 classes 的接近程度来估算标记 classes(class A 或 class B)。因此,例如,如果我知道 M1 和 M4 是 class A,那么位于 M1 和 M4 之间的地图中的所有标记也可以 class 化为 A.

如果我知道标记 M4 是 class A,它的位置是 chr1 13,标记 M7 是 B,位置 16,那么我们可以 class 化所有位置小于等于的标记(13+16)/2=14.5 作为 A,14.5 到 16 之间的所有内容作为 B 在同一条染色体上。所以M5是A,M6可以class化为B。

我有一张标记排序位置的地图

M0  chr1    9
M1  chr1    10
M2  chr1    11
M3  chr1    12
M4  chr1    13
M5  chr1    14
M6  chr1    15
M7  chr1    16
M8  chr2    1
M9  chr2    2
M10 chr2    3
M11 chr2    4

所以给出一个简单的 backbone 的

M1  A
M4  A
M7  B
M8  B
M10 A

如果可能的话,我想估算地图上的其余标记。

所以我想要的输出是

M1  A
M2  A
M3  A
M4  A
M5  A
M6  B
M7  B
M8  B
M9  B
M10 A

我是一名生物学家,正在尝试学习一点点 awk,但我认为这可能只是一个计算问题,我不确定从哪里开始。请帮忙。我可以访问 运行 awk 和 perl 的 unix 集群。请注意,只能在映射到同一染色体的标记之间进行正确插补。

你从来没有回答过我的任何问题,所以这里有一个 Perl 解决方案,可以进行很多猜测

use strict;
use warnings 'all';
use autodie;

my (@markers, %markers);
{
    open my $fh, '<', 'markers.txt';

    while ( <$fh> ) {
        my @marker = split;
        push @markers, \@marker;
        $markers{$marker[0]} = $#markers;
    }
}

my ($i0, $i1);

open my $fh, '<', 'classes.txt';

while ( <$fh> ) {

    my ($marker, $class) = split;

    $i1 = $markers{$marker};
    my $m1 = $markers[$i1];
    push @$m1, $class;

    next unless defined $i0;

    my $m0 = $markers[$i0];

    next if $m0->[1] ne $m1->[1];          # Different chromosomes

    my $mid = ( $m0->[2] + $m1->[2] ) / 2; # Mid point between markers

    for my $m ( @markers[ $i0+1 .. $i1-1 ] ) {
        push @$m, $m->[2] <= $mid ? $m0->[3] : $m1->[3];
    }
}
continue {
    $i0 = $i1;
}

printf "%-4s%-8s%-4d%-s\n", @{$_}[0..2], $_->[3] // '' for @markers;

输出

M0  chr1    9   
M1  chr1    10  A
M2  chr1    11  A
M3  chr1    12  A
M4  chr1    13  A
M5  chr1    14  A
M6  chr1    15  B
M7  chr1    16  B
M8  chr2    1   B
M9  chr2    2   B
M10 chr2    3   A
M11 chr2    4