annotate GetHaplotypesFromPhasedVCF/GetHaplotypesFromPhasedVCF.pl @ 9:98c37a5d67f4 draft

Uploaded
author dereeper
date Wed, 07 Feb 2018 22:08:47 -0500
parents ebb0ac9b6fa9
children 88748d846a20
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
6
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
1 #!/usr/bin/perl
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
2
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
3 use strict;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
4
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
5 my $vcf = $ARGV[0];
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
6 my $out = $ARGV[1];
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
7
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
8 open(O1,">$out.haplo.fas");
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
9 open(O2,">$out.distinct_haplotypes.txt");
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
10
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
11 my %indiv;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
12 my %genes;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
13 my $nb_cols = 0;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
14 my %phasing;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
15 open(my $V,$vcf);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
16 while(<$V>){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
17 my $line = $_;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
18 $line =~s/\n//g;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
19 $line =~s/\r//g;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
20 my @infos = split(/\t/,$line);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
21 if (/#CHROM/){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
22 for (my $i = 9; $i <= $#infos; $i++){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
23 $indiv{$i} = $infos[$i];
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
24 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
25 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
26 else{
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
27 my $gene = $infos[0];
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
28 my $ref = $infos[3];
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
29 my $alt = $infos[4];
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
30 $nb_cols = $#infos;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
31 for (my $i = 9; $i <= $#infos; $i++){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
32 my ($alleles,$depth) = split(":",$infos[$i]);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
33 $alleles =~s/0/$ref/g;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
34 $alleles =~s/1/$alt/g;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
35 my $ind = $indiv{$i};
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
36 #if ($alleles =~/\// && $genes{$gene}){next;}
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
37 $phasing{$gene}{$i}.= $alleles ." ";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
38 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
39 $genes{$gene} = 1;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
40 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
41 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
42 close($V);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
43
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
44 my %haplos;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
45 my $gene_ok = 0;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
46 my $gene_shared = 0;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
47 my $gene_all_shared = 0;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
48 my $nb_gene = 0;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
49 my %haplotypes2;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
50 foreach my $gene(keys(%phasing))
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
51 {
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
52 my $list = "";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
53 my $ok = 0;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
54 my %haplotypes;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
55 for (my $i=9;$i <= $nb_cols;$i++){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
56
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
57 my $alleles = $phasing{$gene}{$i};
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
58 if ($alleles eq "" or $alleles =~ /\./ or $alleles =~/ \w+\//){next;}
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
59 my @snps = split(" ",$alleles);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
60 if (scalar @snps < 2){next;}
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
61 $alleles =~s/\//\|/g;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
62 my @al = split(" ",$alleles);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
63 my $haplo1 = "";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
64 my $haplo2 = "";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
65 foreach my $a(@al){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
66 my ($a1,$a2) = split(/\|/,$a);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
67 $haplo1.= $a1;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
68 $haplo2.= $a2;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
69 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
70 $haplotypes{$haplo1}++;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
71 $haplotypes{$haplo2}++;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
72 my $ind = $indiv{$i};
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
73 $haplotypes2{$gene}{$haplo1}.=$ind."_1,";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
74 $haplotypes2{$gene}{$haplo2}.=$ind."_2,";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
75 $haplos{$gene}{$haplo1}++;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
76 $haplos{$gene}{$haplo2}++;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
77 $list .= ">".$gene."_".$ind."_1\n$haplo1\n";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
78 $list .= ">".$gene."_".$ind."_2\n$haplo2\n";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
79 $ok++;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
80 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
81 #print "$gene $nb_cols $ok\n";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
82 #if ($ok > 13 && $found_M1_M2 == 2){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
83 if ($ok == ($nb_cols - 8)){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
84 $nb_gene++;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
85 print O1 $list;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
86 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
87 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
88 foreach my $gene(sort {$a<=>$b} keys(%haplos)){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
89 print O2 "===$gene===\n";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
90 my $ref_hash = $haplos{$gene};
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
91 my %hash = %$ref_hash;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
92 my $num_haplo = 0;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
93 foreach my $haplo(keys(%hash)){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
94 $num_haplo++;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
95 my $haplo_name = "haplo".$num_haplo;
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
96 my $nb = $haplos{$gene}{$haplo};
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
97 my $ind = $haplotypes2{$gene}{$haplo};
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
98 print O2 $haplo_name.":$nb:".$haplotypes2{$gene}{$haplo}."\n".$haplo."\n";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
99 if ($nb > 1){
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
100 #print "$nb \n";
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
101 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
102 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
103 }
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
104
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
105 close(O1);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
106 close(O1);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
107 #print scalar keys(%haplos);
ebb0ac9b6fa9 planemo upload
gandres
parents:
diff changeset
108