|
2
|
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
|
