|
0
|
1 #!/nfs_exports/apps/64-bit/gnu-apps/perl5.8.9/bin/perl -w
|
|
|
2 use strict;
|
|
|
3 use Carp;
|
|
|
4 use Getopt::Long;
|
|
|
5 use English;
|
|
|
6 use Pod::Usage;
|
|
|
7 use Data::Dumper;
|
|
|
8 use Bio::DB::Sam;
|
|
|
9 use Bio::DB::Sam::Constants;
|
|
|
10 use Bio::SearchIO;
|
|
|
11 use Bio::SeqIO;
|
|
|
12 use File::Temp qw/ tempfile tempdir /;
|
|
|
13 use File::Spec;
|
|
|
14 use File::Path;
|
|
|
15 use File::Copy;
|
|
|
16 use File::Basename;
|
|
|
17 use Cwd;
|
|
|
18
|
|
|
19 # bam related variables
|
|
|
20 # input/output
|
|
|
21 my($file_d, $file_g);
|
|
|
22
|
|
|
23 my ( $help, $man, $version, $usage );
|
|
|
24 my $optionOK = GetOptions(
|
|
|
25 'd|diagnostic=s' => \$file_d,
|
|
|
26 'g|germline=s' => \$file_g,
|
|
|
27 'h|help|?' => \$help,
|
|
|
28 'man' => \$man,
|
|
|
29 'usage' => \$usage,
|
|
|
30 'v|version' => \$version,
|
|
|
31 );
|
|
|
32
|
|
|
33 pod2usage(-verbose=>2) if($man or $usage);
|
|
|
34 pod2usage(1) if($help or $version );
|
|
|
35
|
|
|
36 my $start_dir = getcwd;
|
|
|
37
|
|
|
38 # figure out input file
|
|
|
39 if(!$file_d or !$file_g ) {
|
|
|
40 croak "you need specify diagnostic and germline input file";
|
|
|
41 }
|
|
|
42 my ( @cover_D, @cover_diff);
|
|
|
43 for( my $i = 1; $i < 1000; $i++) {
|
|
|
44 $cover_D[$i] = 0;
|
|
|
45 $cover_diff[$i] = 0;
|
|
|
46 }
|
|
|
47 my($prefix, $path, $suffix) = fileparse($file_d);
|
|
|
48 my $diff_file = $prefix . ".somatic.cover";
|
|
|
49 $diff_file = File::Spec->catfile($path, $diff_file);
|
|
|
50 open my $SOMATIC, ">$diff_file" or croak "can't open $diff_file:$OS_ERROR";
|
|
|
51 open my $IN_D, "<$file_d" or croak "can't open $file_d:$OS_ERROR";
|
|
|
52 open my $IN_G, "<$file_g" or croak "can't open $file_g:$OS_ERROR";
|
|
|
53 my %g_sclip;
|
|
|
54 while( my $line = <$IN_G> ) {
|
|
|
55 chomp $line;
|
|
|
56 my ($chr, $pos, $ort, $cover) = split /\t/, $line;
|
|
|
57 $g_sclip{$chr} = {} if(!exists($g_sclip{$chr}));
|
|
|
58 $g_sclip{$chr}->{$pos} = $ort;
|
|
|
59 }
|
|
|
60 close($IN_G);
|
|
|
61 while( my $line = <$IN_D> ) {
|
|
|
62 chomp $line;
|
|
|
63 my ($chr, $pos, $ort, $cover) = split /\t/, $line;
|
|
|
64 $cover_D[$cover]++;
|
|
|
65 next if(exists ($g_sclip{$chr}->{$pos}) && $g_sclip{$chr}->{$pos} eq $ort);
|
|
|
66 $cover_diff[$cover]++;
|
|
|
67 print $SOMATIC $line, "\n";
|
|
|
68 }
|
|
|
69 close($IN_D);
|
|
|
70 close $SOMATIC;
|
|
|
71 for( my $i = 1; $i < 1000; $i++) {
|
|
|
72 print join("\t", $i, $cover_D[$i], $cover_diff[$i]), "\n";
|
|
|
73 }
|
|
|
74
|
