Mercurial > repos > mcharles > rapsodyn

comparison rapsodyn/PileupVariant.pl @ 29:7b8646f46010 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mcharles
date Wed, 08 Oct 2014 09:06:53 -0400
parents e8e6b962c1f2
children
comparison
equal deleted inserted replaced
28:fefe6d89dbf1 29:7b8646f46010
1 #!/usr/bin/perl 1 #!/usr/bin/perl
2 #V1.0.1 added log, option parameters
2 use strict; 3 use strict;
4 use warnings;
5 use Getopt::Long;
3 6
4 my $inputfile = $ARGV[0]; 7 my $input_pileup_file;
5 open(IF, $inputfile) or die("Can't open $inputfile\n"); 8 my $output_pileup_file;
9 my $log_file;
10
11 my $nb_base_covered=0;
12 my $nb_variant=0;
13 GetOptions (
14 "input_pileup_file=s" => \$input_pileup_file,
15 "log_file=s" => \$log_file
16 ) or die("Error in command line arguments\n");
17
18 open(IN, $input_pileup_file) or die ("Can't open $input_pileup_file\n");
6 19
7 #Extraction des variants 20 #Extraction des variants
8 my $nb_line=0; 21 my $nb_line=0;
9 while (my $line=<IF>){ 22 while (my $line=<IN>){
10 my $test = $line; 23 #print $line;
11 $test =~ s/\$//g; #the read start at this position 24 $nb_base_covered++;
12 $test =~ s/\^.//g; #the read end at this position followed by quality char 25 $line =~ s/\$//g; #the read start at this position
13 my @field = split(/\s+/,$test); 26 $line =~ s/\^.//g; #the read end at this position followed by quality char
27 #print $line;
14 28
29 my @field = split(/\s+/,$line);
15 if ($field[4]=~/[ATGCN]/i){ 30 if ($field[4]=~/[ATGCN]/i){
16 print $line; 31 print $line;
32 $nb_variant++;
17 } 33 }
18 } 34 }
19 close(IF); 35 close(IN);
36
37 open (LF,">$log_file") or die("Can't open $log_file\n");
38 print LF "\n####\t Variant extraction \n";
39 print LF "Position covered :\t$nb_base_covered\n";
40 print LF "Variant detected :\t$nb_variant\n";
41 close (LF);