|
18
|
1
|
|
|
2 #!/usr/bin/perl
|
|
|
3
|
|
|
4 use strict;
|
|
|
5 use Switch;
|
|
|
6 use Getopt::Long;
|
|
|
7 use Bio::SeqIO;
|
|
|
8
|
|
|
9 my $usage = qq~Usage:$0 <args> [<opts>]
|
|
|
10
|
|
|
11 where <args> are:
|
|
|
12
|
|
|
13 -i, --input <VCF input>
|
|
|
14 -o, --out <output basename>
|
|
|
15 ~;
|
|
|
16 $usage .= "\n";
|
|
|
17
|
|
|
18 my ($input,$out);
|
|
|
19
|
|
|
20 GetOptions(
|
|
|
21 "input=s" => \$input,
|
|
|
22 "out=s" => \$out
|
|
|
23 );
|
|
|
24
|
|
|
25
|
|
|
26 die $usage
|
|
|
27 if ( !$input);
|
|
|
28
|
|
|
29
|
|
|
30
|
|
|
31 my $nb_gene = `grep -c mRNA $input`;
|
|
|
32 $nb_gene =~s/\n//g;
|
|
|
33 my $nb_intergenic = `grep -c INTERGENIC $input`;
|
|
|
34 $nb_intergenic =~s/\n//g;
|
|
|
35
|
|
|
36 my $nb_intron = `grep -c INTRON $input`;
|
|
|
37 $nb_intron =~s/\n//g;
|
|
|
38 my $nb_UTR = `grep -c UTR $input`;
|
|
|
39 $nb_UTR =~s/\n//g;
|
|
|
40 my $nb_exon = $nb_gene - $nb_intron - $nb_UTR;
|
|
|
41
|
|
|
42 my $nb_ns = `grep -c NON_SYNONYMOUS_CODING $input`;
|
|
|
43 $nb_ns =~s/\n//g;
|
|
|
44 my $nb_s = $nb_exon - $nb_ns;
|
|
|
45
|
|
|
46
|
|
|
47
|
|
|
48
|
|
|
49 #system("$VCFTOOLS_EXE --vcf $input --remove-filtered-all --out $out --hardy >>vcftools.log 2>&1");
|
|
|
50 system("vcftools --vcf $input --remove-filtered-all --out $out --het >>vcftools.log 2>&1");
|
|
|
51 system("vcftools --vcf $input --remove-filtered-all --out $out --TsTv-summary >>vcftools.log 2>&1");
|
|
|
52 system("vcftools --vcf $input --remove-filtered-all --out $out --missing-indv >>vcftools.log 2>&1");
|
|
|
53
|
|
|
54 open(my $G,">$out.annotation");
|
|
|
55 print $G "Genic $nb_gene\n";
|
|
|
56 print $G "Intergenic $nb_intergenic\n";
|
|
|
57 print $G "========\n";
|
|
|
58 print $G "Intron $nb_intron\n";
|
|
|
59 print $G "Exon $nb_exon\n";
|
|
|
60 print $G "UTR $nb_UTR\n";
|
|
|
61 print $G "========\n";
|
|
|
62 print $G "Non-syn $nb_ns\n";
|
|
|
63 print $G "Synonym $nb_s\n";
|
|
|
64 close($G);
|
|
|
65
|
|
|
66
|
|
|
67
|
|
|
68
|
|
|
69
|
|
|
70
|
|
|
71
|
