根据地图估算 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
我想根据已知标记 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