标记之间插值填补

2024-01-31  本文已影响0人  风知秋

现在是这样的,遗传图谱中的标记是散在分布的,在某些计算中需要特定位点的基因型。

所以需要根据遗传图谱,填补两个 marker 之间所有位点的基因型。

原先的遗传图谱:

插值填补之后的图谱会补充1到308551之间的所有位点,按分位值进行补充。

#!/usr/bin/perl

use strict;

use warnings;

use List::Util qw(min max);

my $filename = shift;

my $output_filename = shift;

open(my $input_fh, '<', $filename) or die "can not open the file '$filename' $!";

open(my $output_fh, '>', $output_filename) or die "can not creat the file '$output_filename' $!";

my @data;

while (my $line = <$input_fh>) {

    chomp $line;

    my @fields = split(/\t/, $line);

    push @data, \@fields;

}

# 填充位点并计算分位数,并将结果写入新文件

for my $i (0..$#data-1) {

    my ($start_chr, $start_pos, $start_val) = @{$data[$i]};

    my ($end_chr, $end_pos, $end_val) = @{$data[$i+1]};

    # 计算填充的行数

    my $num_rows = $end_pos - $start_pos;

    # 计算分位数

    my $min_val = min($start_val, $end_val);

    my $max_val = max($start_val, $end_val);

    # 填充位点并写入新文件

    for my $j (0..$num_rows-1) {

        my $interpolated_val = $start_val + ($end_val - $start_val) * ($j / $num_rows);

        my $pos = $start_pos + $j;

        print $output_fh "$start_chr\t$pos\t$interpolated_val\n";

    }

}

my ($last_chr, $last_pos, $last_val) = @{$data[$#data]};

print $output_fh "$last_chr\t$last_pos\t$last_val\n";

close($input_fh);

close($output_fh);

上一篇 下一篇

猜你喜欢

热点阅读