Mercurial > repos > jjohnson > crest

comparison countDiff.pl @ 0:acc8d8bfeb9a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jjohnson
date Wed, 08 Feb 2012 16:59:24 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:acc8d8bfeb9a
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