修改 VCF 文件中错误编码的基因型

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

# 在计算的时候发现,VCF 文件中有些个体的基因型文件为  '.'  而不是  './.'   这会导致某些计算中报错,比如说:

Caused by: java.lang.IllegalArgumentException: ERROR: inconsistent number of alleles for sample NEP0204 at marker [chr01 79336 . T  A]

不清楚什么软件可以直接处理,就写了个脚本对此进行修改;

#!/usr/bin/perl

use strict;

use warnings;

open(my $fh, '<', 'chr01.recode.vcf') or die "无法打开文件: $!";

open(my $output_fh, '>', 'output.vcf') or die "无法创建输出文件: $!";

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

    chomp $line;

    if ($line =~ /^#/ || $line =~ /^\s*$/) {

        print $output_fh "$line\n";

        next;

    }

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

    my $format = $fields[8];

    my @format_fields = split(':', $format);

    my $gt_index = 0;

    for (my $i = 0; $i < scalar(@format_fields); $i++) {

        if ($format_fields[$i] eq 'GT') {

            $gt_index = $i;

            last;

        }

    }

    for (my $i = 9; $i < scalar(@fields); $i++) {

        my $gt = $fields[$i];

my @gt_fields = split(':', $gt);

if ($gt_fields[$gt_index] eq '.'){

$gt_fields[$gt_index] = './.';

}

my $new_gt = join(':', @gt_fields);

$fields[$i] = $new_gt;

    }

    my $new_line = join("\t", @fields);

    print $output_fh "$new_line\n";

}

close($fh);

close($output_fh);

上一篇下一篇

猜你喜欢

热点阅读