美文网首页
修改 VCF 文件中错误编码的基因型

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

作者: 风知秋 | 来源:发表于2024-01-30 17:13 被阅读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);

    相关文章

      网友评论

          本文标题:修改 VCF 文件中错误编码的基因型

          本文链接:https://www.haomeiwen.com/subject/hqfoodtx.html