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