Mercurial > repos > geert-vandeweyer > varscan_wrapper

comparison varscan_somatic.pl @ 2:90c9b7a6e58b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author geert-vandeweyer
date Fri, 07 Mar 2014 06:18:29 -0500
parents
children
comparison
equal deleted inserted replaced
1:e5831cc9185f 2:90c9b7a6e58b
1 #!/usr/bin/perl
2
3
4 use strict;
5 use Cwd;
6
7 die qq(
8 Bad numbr of inputs
9
10 ) if(!@ARGV);
11
12 my $options ="";
13 my $normal="";
14 my $command="";
15 my $tumor="";
16 my $tumorbam = "";
17 my $output="";
18 my $snp="";
19 my $indel="";
20 my $working_dir = cwd();
21 my $log = '';
22
23 foreach my $input (@ARGV)
24 {
25 my @tmp = split "::", $input;
26 if($tmp[0] eq "COMMAND")
27 {
28 $command = $tmp[1];
29 }
30 elsif($tmp[0] eq "NORMAL")
31 {
32 $normal = $tmp[1];
33 }
34 elsif($tmp[0] eq "TUMOR")
35 {
36 $tumor = $tmp[1];
37 }
38 elsif($tmp[0] eq "TUMORBAM")
39 {
40 $tumorbam = $tmp[1];
41 }
42 elsif($tmp[0] eq "OPTION")
43 {
44 my @p = split(/\s+/,$tmp[1]);
45 if ($p[0] =~ m/validation|strand-filter|output-vcf/ && $p[1] == 0) {
46 next;
47 }
48 $options = "$options ${tmp[1]}";
49 }
50 elsif($tmp[0] eq "OUTPUT")
51 {
52 $output = $tmp[1];
53 }
54 elsif($tmp[0] eq "SNP")
55 {
56 $snp = $tmp[1];
57 }
58 elsif($tmp[0] eq "INDEL")
59 {
60 $indel = $tmp[1];
61 }
62
63 elsif($tmp[0] eq "LOG")
64 {
65 $log = $tmp[1];
66 }
67
68 else
69 {
70 die("Unknown Input: $input\n");
71 }
72 }
73
74 ## VCF OUTPUT
75 if ($output ne '') {
76 system ("$command $normal $tumor $options --output-snp $working_dir/out.snp --output-indel $working_dir/out.indel 2>$log");
77 my $indels = "$working_dir/out.indel.vcf";
78 my $snps = "$working_dir/out.snp.vcf";
79
80 system("grep -v '^\#' $indels | grep -v '^chrom position' >> $snps");
81 my @chr_ord = chromosome_order($tumorbam);
82
83 vs2vcf($snps, $output,\@chr_ord);
84 }
85 ## SNP/INDEL OUTPUT
86 else {
87 system ("$command $normal $tumor $options --output-snp $snp --output-indel $indel 2>$log");
88 }
89
90 sub vs2vcf
91 {
92
93 #
94 # G l o b a l v a r i a b l e s
95 #
96 my $version = '0.1';
97
98 #
99 # Read in file
100 #
101 my $input = shift;
102 my $output = shift;
103 my $chr_ord = shift;
104 open(IN, $input) or die "Can't open $input': $!\n";
105 open(OUT, ">$output") or die "Can't create $output': $!\n";
106 my %output;
107
108 while ( <IN> )
109 {
110 if ( /^#/ )
111 {
112 print OUT;
113 next;
114 }
115 chomp;
116 my $line = $_;
117
118 my @flds = split ( "\t", $line );
119 my $ref = $flds[3];
120 my $alt = $flds[4];
121 #
122 # Deletion of bases
123 #
124 if ( $alt =~ /^\-/ )
125 {
126 ($flds[3], $flds[4]) = ($ref.substr($alt,1), $ref);
127 }
128
129 #
130 # Insertion of bases
131 #
132 if ( $alt =~ /^\+/ )
133 {
134 $flds[4] = $ref.substr($alt,1);
135 }
136 ## insert dot for reference positions.
137 if ($flds[4] eq '') {
138 $flds[4] = '.';
139 }
140 print OUT join( "\t", @flds),"\n" unless defined $chr_ord;
141 $output{$flds[0]}{$flds[1]} = join( "\t", @flds)."\n" if defined $chr_ord;
142 }
143 close(IN);
144 # if chromosome order given return in sorted order
145 if(defined $chr_ord)
146 {
147 for my $chrom (@{ $chr_ord })
148 {
149 for my $pos (sort {$a<=>$b} keys %{ $output{$chrom} })
150 {
151 print OUT $output{$chrom}{$pos};
152 }
153 }
154 }
155 close(OUT);
156 }
157
158
159 sub chromosome_order
160 {
161 my $input = shift;
162 # calculate flagstats
163 my $COMM = "samtools view -H $input | grep '^\@SQ'";
164 my @SQ = `$COMM`;
165 chomp @SQ;
166 for(my $i = 0; $i <= $#SQ; $i++)
167 {
168 $SQ[$i] =~ s/^\@SQ\tSN:(.*?)\tLN:\d+$/$1/;
169 }
170 return(@SQ);
171 }
172
173