# HG changeset patch
# User mcharles
# Date 1422313852 18000
# Node ID 93e6f2af1ce25545d5ffcf5e6d9773ed92d10f60
# Parent  827da1a9a32626679f945dbb54ef3418242b408c
Uploaded

diff -r 827da1a9a326 -r 93e6f2af1ce2 rapsodyn/CreateMatrix.pl
--- 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;
+
+
diff -r 827da1a9a326 -r 93e6f2af1ce2 rapsodyn/CreateMatrix.xml
--- a/rapsodyn/CreateMatrix.xml	Mon Jan 19 10:38:29 2015 -0500
+++ b/rapsodyn/CreateMatrix.xml	Mon Jan 26 18:10:52 2015 -0500
@@ -1,11 +1,11 @@
 <tool id="CreateMatrix" name="CreateMatrix" version="1.00">
 <description>Create Genotyping Matrix</description>
 <command interpreter="perl">
-	CreateMatrix.pl -input_pileup_file $input_pileup_file -input_variant_file $input_variant_file -input_variant_unique_file $input_variant_unique_file > $output_file 
+	CreateMatrix.pl -input_name $input_name -input_pileup_file $input_pileup_file -input_variant_unique_file $input_variant_unique_file > $output_file 
 </command>
 <inputs>
+	<param name="input_name" size="30" type="text" label="Genotype name"/>
 	<param name="input_pileup_file"  type="data" format="pileup" label="Select suitable PILEUP files from your history"/>
-	<param name="input_variant_file"  type="data" format="pileup" label="Select suitable VARIANT files from your history"/>
 	<param name="input_variant_unique_file"  type="data" format="pileup" label="Select suitable VARIANT UNIQUE files from your history"/>
 	
 </inputs>