2
|
1
|
|
2 #!/usr/bin/perl
|
|
3
|
|
4 use strict;
|
|
5 use Getopt::Long;
|
3
|
6 use Bio::SeqIO;
|
2
|
7
|
|
8 my $usage = qq~Usage:$0 <args> [<opts>]
|
|
9
|
|
10 where <args> are:
|
|
11
|
|
12 -i, --input <VCF input>
|
|
13 -o, --out <output basename>
|
|
14 -s, --step <step size>
|
|
15 -w, --window <window size>
|
|
16
|
|
17 <opts> are:
|
|
18 -g, --group <group file>
|
|
19 ~;
|
|
20 $usage .= "\n";
|
|
21
|
|
22 my ($input,$out);
|
|
23 my $window = 200000;
|
|
24 my $step = 200000;
|
|
25 my $groupfile;
|
|
26 GetOptions(
|
|
27 "input=s" => \$input,
|
|
28 "out=s" => \$out,
|
|
29 "step=s" => \$step,
|
|
30 "window=s" => \$window,
|
|
31 "group=s" => \$groupfile
|
|
32 );
|
|
33
|
|
34
|
|
35
|
|
36
|
|
37 die $usage
|
|
38 if ( !$input);
|
|
39
|
|
40 my %hash;
|
|
41 my $cmd_part = "";
|
|
42 if ($groupfile && -e $groupfile)
|
|
43 {
|
|
44 open(my $G,$groupfile) or die "Cannot open $groupfile: $!";
|
|
45 while(<$G>)
|
|
46 {
|
|
47 my $line = $_;
|
|
48 chomp($line);
|
|
49 $line=~s/\r//g;
|
|
50 $line=~s/\n//g;
|
|
51 my @infos = split(/;/,$line);
|
|
52 if ($infos[0] && $infos[1])
|
|
53 {
|
|
54 $hash{$infos[1]} .= " --indv " . $infos[0];
|
|
55 $cmd_part .= " --indv " . $infos[0];
|
|
56 }
|
|
57 }
|
|
58 close($G);
|
|
59 }
|
|
60
|
|
61
|
|
62 if ($step =~/^(\d+)\s*$/){
|
|
63 $step = $1;
|
|
64 }
|
|
65 else{
|
|
66 die "Error: step size must be an integer\n";
|
|
67 }
|
|
68 if ($window =~/^(\d+)\s*$/){
|
|
69 $window = $1;
|
|
70 }
|
|
71 else{
|
|
72 die "Error: window size must be an integer\n";
|
|
73 }
|
|
74
|
|
75
|
3
|
76 my $input_format = "--vcf $input.sorted";
|
|
77 if ($input =~/\.vcf/){
|
|
78 system("vcf-sort $input >$input.sorted");
|
|
79 $input_format = "--vcf $input.sorted";
|
|
80 }
|
|
81 elsif ($input =~/\.bcf/){
|
|
82 $input_format = "--bcf $input";
|
|
83 }
|
|
84 system("vcftools $input_format --out $out --window-pi $window --window-pi-step $step >>$out.vcftools.log 2>&1");
|
|
85 system("vcftools $input_format --out $out --TajimaD $window >>$out.vcftools.log 2>&1");
|
|
86 system("vcftools $input_format --out $out --TsTv $window >>$out.vcftools.log 2>&1");
|
|
87 system("vcftools $input_format --out $out --SNPdensity $window >>$out.vcftools.log 2>&1");
|
2
|
88
|
|
89
|
|
90
|
3
|
91 ##################################################
|
|
92 # Tajima distribution
|
|
93 ##################################################
|
|
94 my %counts3;
|
|
95 my @values;
|
|
96 open(my $TAJIMA,"$out.Tajima.D");
|
|
97 <$TAJIMA>;
|
|
98 my $max = 0;
|
|
99 while(<$TAJIMA>){
|
|
100 my $line = $_;
|
|
101 $line =~s/\n//g;
|
|
102 $line =~s/\r//g;
|
|
103 my @inf = split(/\t/,$line);
|
|
104 my $val = $inf[3];
|
|
105 my $arrondi = sprintf("%.1f",$val);
|
|
106 $counts3{$arrondi}++;
|
|
107 }
|
|
108 close($TAJIMA);
|
|
109
|
|
110
|
|
111
|
|
112 open(my $M3,">$out.tajima_distrib");
|
|
113 print $M3 "Values nb of values\n";
|
|
114 for (my $i = -2.5; $i <= 2.5; $i=$i+0.1)
|
|
115 {
|
|
116 my $arrondi = sprintf("%.1f",$i);
|
|
117 my $nb = 0;
|
|
118 if ($counts3{$arrondi})
|
|
119 {
|
|
120 $nb = $counts3{$arrondi};
|
|
121 }
|
|
122 print $M3 "Tajima $arrondi $nb\n";
|
|
123 }
|
|
124 close($M3);
|
|
125
|
2
|
126
|
|
127 if (keys(%hash) > 0)
|
|
128 {
|
|
129 my $files_pi = "";
|
|
130 my $files_dtajima = "";
|
3
|
131 my $cmd_fst = "--fst-window-size $window";
|
|
132 my $cmd_fst_by_marker = "";
|
2
|
133 foreach my $pop(sort(keys(%hash)))
|
|
134 {
|
3
|
135 my $cmd_part = $hash{$pop};
|
|
136
|
|
137
|
|
138 open(my $POP,">$input.$pop.txt");
|
|
139 my @ind_of_pop = split(" --indv ",$cmd_part);
|
|
140 print $POP join("\n",@ind_of_pop);
|
|
141 close($POP);
|
|
142 $cmd_fst .= " --weir-fst-pop $input.$pop.txt";
|
|
143 $cmd_fst_by_marker .= " --weir-fst-pop $input.$pop.txt";
|
|
144
|
2
|
145 my $cmd_part = $hash{$pop};
|
3
|
146 system("vcftools $input_format --remove-filtered-all --out $out.$pop --window-pi $window --window-pi-step $step $cmd_part --maf 0.001 >>$out.vcftools.log 2>&1");
|
2
|
147 my $sed_cmd = "sed -i \"s\/PI\/$pop\/g\" $out.$pop.windowed.pi";
|
|
148 system($sed_cmd);
|
|
149 $files_pi .= "$out.$pop.windowed.pi ";
|
|
150
|
3
|
151 system("vcftools $input_format --remove-filtered-all --out $out.$pop --SNPdensity $window $cmd_part --maf 0.001 >>$out.vcftools.log 2>&1");
|
2
|
152
|
3
|
153 system("vcftools $input_format --remove-filtered-all --out $out.$pop --TajimaD $window $cmd_part --maf 0.001 >>$out.vcftools.log 2>&1");
|
2
|
154 my $sed_cmd = "sed -i \"s\/TajimaD\/$pop\/g\" $out.$pop.Tajima.D";
|
|
155 system($sed_cmd);
|
|
156 $sed_cmd = "sed -i \"s/nan/0/g\" $out.Tajima.D";
|
|
157 system($sed_cmd);
|
|
158 $files_dtajima .= "$out.$pop.Tajima.D ";
|
|
159
|
3
|
160 system("vcftools $input_format --remove-filtered-all --out $out.$pop --TsTv $window $cmd_part --maf 0.001 >>$out.vcftools.log 2>&1");
|
2
|
161 }
|
|
162 system("paste $files_pi >>$out.combined.pi");
|
|
163 my $awk_cmd = "awk {'print \$1\"\t\"\$2\"\t\"\$5\"\t\"\$10'} $out.combined.pi >$out.combined.pi.txt";
|
|
164 system($awk_cmd);
|
|
165
|
|
166 system("paste $files_dtajima >>$out.combined.dtajima");
|
3
|
167 $awk_cmd = "awk {'print \$1\"\t\"\$2\"\t\"\$4\"\t\"\$8'} $out.combined.dtajima >$out.combined.dtajima.txt";
|
|
168 system($awk_cmd);
|
|
169
|
|
170 system("vcftools $input_format --remove-filtered-all --out $out.fst $cmd_fst --maf 0.001 >>$out.vcftools.log 2>&1");
|
|
171 $awk_cmd = "awk {'print \$1\"\t\"\$2\"\t\"\$5'} $out.fst.windowed.weir.fst >$out.fst.txt";
|
2
|
172 system($awk_cmd);
|
3
|
173
|
|
174 system("vcftools $input_format --remove-filtered-all --out $out.fst $cmd_fst_by_marker --maf 0.001 >>$out.vcftools.log 2>&1");
|
|
175
|
|
176 my %genes;
|
|
177 open(my $V,$input);
|
|
178 my $go = 0;
|
|
179 while(<$V>){
|
|
180 if ($go){
|
|
181 my @infos = split(/\t/,$_);
|
|
182 my $chr = $infos[0];
|
|
183 my $pos = $infos[1];
|
|
184 my $annot = $infos[7];
|
|
185 my $gene = "#";
|
|
186 if ($annot =~/EFF=\w+\((.*)\)/){
|
|
187 my @parts = split(/\|/,$1);
|
|
188 $gene = $parts[4];
|
|
189 }
|
|
190 $genes{"$chr:$pos"} = $gene;
|
|
191 }
|
|
192 if (/^#CHROM/){$go = 1;}
|
|
193 }
|
|
194
|
|
195 my $cumul = 0;
|
|
196 my $previous_pos;
|
|
197 my %chrs;
|
|
198 open(my $F,"$out.fst.weir.fst");
|
|
199 open(my $O,">$out.fst.by_marker.txt");
|
|
200 open(my $O2,">$out.fst.by_marker.genes.txt");
|
|
201 print $O "Chromosome Marker Position Fst\n";
|
|
202 print $O2 "Chromosome Marker Position Fst Gene\n";
|
|
203 while(<$F>){
|
|
204 my $line = $_;
|
|
205 chomp($line);
|
|
206 my ($chr,$pos,$val) = split(/\t/,$line);
|
|
207 if (!$chrs{$chr}){
|
|
208 $cumul = $previous_pos;
|
|
209 $chrs{$chr} = 1;
|
|
210 }
|
|
211 if ($val > 0.6 && $val !~/-/ && $val ne 'nan'){
|
|
212 my $x = $pos+$cumul;
|
|
213 print $O2 "$chr $chr:$pos $pos $val " . $genes{"$chr:$pos"} . "\n";
|
|
214 print $O "$chr $chr:$pos $x $val\n";
|
|
215 $previous_pos = $x;
|
|
216 }
|
|
217
|
|
218 }
|
|
219 close($F);
|
|
220 close($O);
|
|
221 close($O2);
|
|
222
|
|
223 #$awk_cmd = "awk {'if (\$3>0.6 && \$3<=1 && \$3!~\"nan\"){print \$1\"\t\"\$2\"\t\"\$2\"\t\"\$3}'} $out.fst.weir.fst >$out.fst.by_marker.txt";
|
|
224 #system($awk_cmd);
|
2
|
225 }
|
|
226
|
|
227
|
|
228
|
|
229
|
|
230
|
|
231
|