Mercurial > repos > fcaramia > varscan

comparison varscan_somatic.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
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 $output="";
17 my $working_dir = cwd();
18 my $snp = "$working_dir/output.snp.vcf";
19 my $indels = "$working_dir/output.indel.vcf";
20
21 foreach my $input (@ARGV)
22 {
23 my @tmp = split "::", $input;
24 if($tmp[0] eq "COMMAND")
25 {
26 $command = $tmp[1];
27 }
28 if($tmp[0] eq "NORMAL")
29 {
30 $normal = $tmp[1];
31 }
32 elsif($tmp[0] eq "TUMOR")
33 {
34 $tumor = $tmp[1];
35 }
36 elsif($tmp[0] eq "OPTION")
37 {
38 $options = "$options ${tmp[1]}";
39 }
40 elsif($tmp[0] eq "OUTPUT")
41 {
42 $output = $tmp[1];
43 }
44
45 else
46 {
47 die("Unknown Input: $input\n");
48 }
49 }
50
51 system ("$command $normal $tumor $options ");
52 system("grep -v '^\#' $indels | grep -v '^chrom position' >> $snp");
53
54 my @chr_ord = chromosome_order($tumor);
55
56 vs2vcf($snp, $output,\@chr_ord);
57
58
59 sub vs2vcf
60 {
61
62 #
63 # G l o b a l v a r i a b l e s
64 #
65 my $version = '0.1';
66
67 #
68 # Read in file
69 #
70 my $input = shift;
71 my $output = shift;
72 my $chr_ord = shift;
73 open(IN, $input) or die "Can't open $input': $!\n";
74 open(OUT, ">$output") or die "Can't create $output': $!\n";
75 my %output;
76
77 while ( <IN> )
78 {
79 if ( /^#/ )
80 {
81 print OUT;
82 next;
83 }
84 chomp;
85 my $line = $_;
86
87 my @flds = split ( "\t", $line );
88 my $ref = $flds[3];
89 my $alt = $flds[4];
90 #
91 # Deletion of bases
92 #
93 if ( $alt =~ /^\-/ )
94 {
95 ($flds[3], $flds[4]) = ($ref.substr($alt,1), $ref);
96 }
97
98 #
99 # Insertion of bases
100 #
101 if ( $alt =~ /^\+/ )
102 {
103 $flds[4] = $ref.substr($alt,1);
104 }
105 print OUT join( "\t", @flds),"\n" unless defined $chr_ord;
106 $output{$flds[0]}{$flds[1]} = join( "\t", @flds)."\n" if defined $chr_ord;
107 }
108 close(IN);
109 # if chromosome order given return in sorted order
110 if(defined $chr_ord)
111 {
112 for my $chrom (@{ $chr_ord })
113 {
114 for my $pos (sort {$a<=>$b} keys %{ $output{$chrom} })
115 {
116 print OUT $output{$chrom}{$pos};
117 }
118 }
119 }
120 close(OUT);
121 }
122
123
124 sub chromosome_order
125 {
126 my $input = shift;
127 # calculate flagstats
128 my $COMM = "samtools view -H $input | grep '^\@SQ'";
129 my @SQ = `$COMM`;
130 chomp @SQ;
131 for(my $i = 0; $i <= $#SQ; $i++)
132 {
133 $SQ[$i] =~ s/^\@SQ\tSN:(.*?)\tLN:\d+$/$1/;
134 }
135 return(@SQ);
136 }
137
138