-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathVcf.phase2fasta.pl
75 lines (66 loc) · 1.5 KB
/
Vcf.phase2fasta.pl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#! perl
use warnings;
use strict;
my $vcf = shift;
my $ref = shift;
my %r;
%r = &read_ref($ref) if ($ref);
my %fa;
my %fa2;
open IN,'<',$vcf or die "$!";
my @head;
push @head, 0..8;
D:while(<IN>){
next if /^##/;
chomp;
if(/^#C/){
my @line = (split/\t/,$_);
for(my $i = 9;$i< @line;$i += 1){
push @head, $line[$i];
}
}else{
my @line = split/\t/,$_;
my $p;
if($ref){
next D if (!exists $r{$line[0]}{$line[1]});
}
my @base = ($line[3] , $line[4]); # 0 for ref, 1.. for ale
for(my $i = 9;$i< @line;$i += 1){
if ($line[$i] =~ /^\.\/\./){
$p = "--";
$fa2{$head[$i]}{$line[0]}{$line[1]} = $p;
#push @{$fa{$head[$i]}}, $p;
}else{
my @info = split/:/,$line[$i];
$info[0] =~ /(\d+)\/(\d+)/;
(my $GT1,my $GT2) = ($1,$2);
$p = "$base[$GT1]$base[$GT2]";
$fa2{$head[$i]}{$line[0]}{$line[1]} = $p;
}
}
}
}
for my $sam (sort {$a cmp $b} keys %fa2){
print ">$sam\n";
for my $chr (sort {$a cmp $b} keys %{$fa2{$sam}}){
for my $loc (sort {$a <=> $b} keys %{$fa2{$sam}{$chr}}){
print "$fa2{$sam}{$chr}{$loc}";
}
}
print "\n";
}
sub read_ref{
my $f = shift @_;
open R,'<', $f or die "$!";
my %h;
while(<R>){
chomp;
my @l=split /\s+/,$_;
#(my $chr,my $loc)=(split /\t/,$_ )[0,1];
my $chr=$l[0];
my $loc=$l[1];
$h{$chr}{$loc} = 1;
}
close R;
return %h;
}