Mercurial > repos > mcharles > rapsosnp
diff rapsodyn/CreateMatrix.pl @ 14:93e6f2af1ce2 draft
Uploaded
author | mcharles |
---|---|
date | Mon, 26 Jan 2015 18:10:52 -0500 |
parents | 0a6c1cfe4dc8 |
children |
line wrap: on
line diff
--- a/rapsodyn/CreateMatrix.pl Mon Jan 19 10:38:29 2015 -0500 +++ b/rapsodyn/CreateMatrix.pl Mon Jan 26 18:10:52 2015 -0500 @@ -5,16 +5,166 @@ use Getopt::Long; my $input_pileup_file; -my $input_variant_file; my $input_variant_unique_file; +my $input_name; GetOptions ( +"input_name=s" => \$input_name, "input_pileup_file=s" => \$input_pileup_file, -"input_variant_file=s" => \$input_variant_file, "input_variant_unique_file=s" => \$input_variant_unique_file ) or die("Error in command line arguments\n"); +print "Chrom\tPos\tRef\t$input_name\n"; -print "$input_pileup_file\n$input_variant_file\n$input_variant_unique_file\n"; +my %variant_unique; +open (VU,$input_variant_unique_file) or die ("Can't open variant unique file\n"); +while (my $line=<VU>){ + if ($line!~/^\s*$/){ + my @fields = split (/\t+/,$line); + my $chr; + my $pos; + if ($fields[0]=~/^\s*([\w\-]+)\s*$/){ + $chr = $1; + } + else { + print STDERR "Unable to detect chromosome in pileup file\n$line",$fields[0],"\n",$fields[1],"\n"; + exit(0); + } + if ($fields[1]=~/^\s*(\d+)\s*$/){ + $pos = $1; + } + else { + print STDERR "Unable to detect position in pileup file\n$line",$fields[0],"\n",$fields[1],"\n"; + exit(0); + } + if ($#fields>=4){ + $variant_unique{"$chr#$pos"}=$fields[4]; + } + else { + print STDERR "Unrecognized pileup format, need at least 5 columns\n$line\n"; + } + } +} +close VU; + +my %covered; +my $compt=0; +open (PU,$input_pileup_file) or die ("Can't open pileup file\n"); +while (my $line=<PU>){ + #$compt++; + #if ($compt>200000){last;} + if ($line!~/^\s*$/){ + my @fields = split (/\t+/,$line); + my $chr; + my $pos; + my $ref; + my $depth; + my $pile; + + if ($#fields>=4){ + if ($fields[0]=~/^\s*([\w\-]+)\s*$/){ + $chr = $1; + } + else { + print STDERR "Unable to detect chromosome in pileup file\n$line\n"; + exit(0); + } + if ($fields[1]=~/^\s*(\d+)\s*$/){ + $pos = $1; + } + else { + print STDERR "Unable to detect position in pileup file\n$line\n"; + exit(0); + } + + if ($fields[2]=~/^\s*([ATGCNX])\s*$/i){ + $ref = $1; + } + else { + print STDERR "Unable to detect reference base in pileup file\n$line\n"; + exit(0); + } + + if ($fields[3]=~/^\s*(\d+)\s*$/i){ + $depth = $1; + } + else { + print STDERR "Unable to detect depth in pileup file\n$line\n"; + exit(0); + } + + $pile = $fields[4]; + + + if ($variant_unique{"$chr#$pos"}) { + $pile = $variant_unique{"$chr#$pos"}; + $pile =~ s/\$//g; #the read start at this position + $pile =~ s/\^.//g; #the read end at this position followed by quality char + my $ori_pile = $pile; + + my @detail = split(//,$pile); + my $current_in=""; + my $current_del=""; + my $current_length=0; + my $noindel_pile=""; + my @IN; + my @DEL; + + $noindel_pile = $pile; + while ($noindel_pile =~ /\+(\d+)/){ + my $length = $1; + $noindel_pile =~ s/(\+\d+[ATGCNX]{$length})//i ; + push(@IN,$1); + #print "test : $pile / $noindel_pile\n"; + } + while ($noindel_pile =~ /\-(\d+)/){ + my $length = $1; + $noindel_pile =~ s/(\-\d+[ATGCNX]{$length})//i ; + push(@DEL,$1); + #print "test : $pile / $noindel_pile\n"; + } + + my %variant_type; + + my $indel_pile = $ori_pile; + $variant_type{"IN"} = ($indel_pile =~ s/\+//gi); + $variant_type{"DEL"} = ($indel_pile =~ s/\-//gi); + + my $base_pile = $noindel_pile; + $variant_type{"A"} = ($base_pile =~ s/a//gi); + $variant_type{"T"} = ($base_pile =~ s/t//gi); + $variant_type{"C"} = ($base_pile =~ s/c//gi); + $variant_type{"G"} = ($base_pile =~ s/g//gi); + $variant_type{"N"} = ($base_pile =~ s/n//gi); + + my $top_number=0; + my $top_variant=""; + foreach my $key (sort{$variant_type{$b}<=>$variant_type{$a}} keys %variant_type){ + if ($variant_type{$key} > 0){ + if ($variant_type{$key} > $top_number){ + $top_number = $variant_type{$key}; + $top_variant = $key; + } + elsif ($variant_type{$key} == $top_number){ + $top_variant .= "/$key"; + } + } + + } + print "$chr\t$pos\t$ref\t$top_variant\n"; + + } + else { + print "$chr\t$pos\t$ref\t$ref\n"; + } + } + else { + print STDERR "Unrecognized pileup format, need at least 5 columns\n$line\n"; + } + } +} +close PU; + +