Mercurial > repos > fcaramia > varscan

comparison varscan_mpileup.pl @ 2:51969e284317 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author fcaramia
date Thu, 20 Jun 2013 00:00:22 -0400
parents
children
comparison
equal deleted inserted replaced
1:893954763c0e 2:51969e284317
1 #!/usr/bin/perl
2
3 use strict;
4 use Cwd;
5
6 die qq(
7 Bad numbr of inputs
8
9 ) if(!@ARGV);
10
11 my $options ="";
12 my $file="";
13 my $command="";
14 my $output="";
15 my $working_dir = cwd();
16 my $temp_vcf = "$working_dir/temp";
17 my $log="";
18
19 foreach my $input (@ARGV)
20 {
21 my @tmp = split "::", $input;
22 if($tmp[0] eq "COMMAND")
23 {
24 $command = $tmp[1];
25 }
26 elsif($tmp[0] eq "INPUT")
27 {
28 $file = $tmp[1];
29 }
30 elsif($tmp[0] eq "OPTION")
31 {
32 $options = "$options ${tmp[1]}";
33 }
34 elsif($tmp[0] eq "OUTPUT")
35 {
36 $output = $tmp[1];
37 }
38 elsif($tmp[0] eq "LOG")
39 {
40 $log = $tmp[1];
41 }
42 else
43 {
44 die("Unknown Input: $input\n");
45 }
46 }
47
48 system ("$command $file $options 1>$temp_vcf 2>$log");
49
50 vs2vcf($temp_vcf, $output);
51
52
53 sub vs2vcf
54 {
55
56 #
57 # G l o b a l v a r i a b l e s
58 #
59 my $version = '0.1';
60
61 #
62 # Read in file
63 #
64 my $input = shift;
65 my $output = shift;
66 my $chr_ord = shift;
67 open(IN, $input) or die "Can't open $input': $!\n";
68 open(OUT, ">$output") or die "Can't create $output': $!\n";
69 my %output;
70
71 while ( <IN> )
72 {
73 if ( /^#/ )
74 {
75 print OUT;
76 next;
77 }
78 chomp;
79 my $line = $_;
80
81 my @flds = split ( "\t", $line );
82 my $ref = $flds[3];
83 my $alt = $flds[4];
84 #
85 # Deletion of bases
86 #
87 if ( $alt =~ /^\-/ )
88 {
89 ($flds[3], $flds[4]) = ($ref.substr($alt,1), $ref);
90 }
91
92 #
93 # Insertion of bases
94 #
95 if ( $alt =~ /^\+/ )
96 {
97 $flds[4] = $ref.substr($alt,1);
98 }
99 print OUT join( "\t", @flds),"\n" unless defined $chr_ord;
100 $output{$flds[0]}{$flds[1]} = join( "\t", @flds)."\n" if defined $chr_ord;
101 }
102 close(IN);
103 # if chromosome order given return in sorted order
104 if(defined $chr_ord)
105 {
106 for my $chrom (@{ $chr_ord })
107 {
108 for my $pos (sort {$a<=>$b} keys %{ $output{$chrom} })
109 {
110 print OUT $output{$chrom}{$pos};
111 }
112 }
113 }
114 close(OUT);
115 }
116