fcmp

#!/usr/bin/perl

use warnings;
use strict;
use Getopt::Long;
use fralib;
use POSIX qw(ceil floor);
use File::Path;
use File::Basename;

use Pod::Usage;

=head1 NAME

fcmp

=head1 SYNOPSIS

 fcmp [options] genotype-file-a genotype-file-b
	
  -v  verbose
  -d  debug
  -c  compare list: contains a list of SNPs or samples
      (Only 1 genotype file is expected if this list is given)
      a)sample-id-a
      b)sample-id-b
      OR
      a)snp-id-a
      b)snp-id-b
      e.g. snp-id-a        snp-id-b
           SNP_A-4258036   AFFX-SNP_12264458
           SNP_A-4301074   AFFX-SNP_9467511
           SNP_A-1932641   AFFX-SNP_12305905
  -o  output file

 example: fcmp -o pscalare-paltum.txt pscalare.gt paltum.gt
          fcmp -c pscalare-compare.txt pscalare.gt paltum.gt
 
 Compares the difference in genotype between common samples(or defined sample pairs) and markers(or defined marker pairs) in 1 or 2 genotype files.
 
 It will be a good idea to perhaps separate the functionalities of this progamme.  Mostly because the code is in a mess now. (Working though)
 
=head1 DESCRIPTION

=cut

#option variables
my $help;
my $outFile;
my $genotypeFileA;
my $genotypeFileB;
my $isGt;
my $compareListFile;
my %COMPARE;
my @COMPARE_A;
my @COMPARE_B;
my $colNo;
my %label2col;
my $fposA;
my $colNoA;
my @fieldsA;
my $fposB;
my $colNoB;
my @fieldsB;
my %A_COL_ELEMENTS;
my %A_ROW_ELEMENTS;
my %B_COL_ELEMENTS;
my %B_ROW_ELEMENTS;

#initialize options
Getopt::Long::Configure ('bundling');

if(!GetOptions ('h'=>\$help,'o=s'=>\$outFile,'c=s'=>\$compareListFile) 
   || (defined($compareListFile) && scalar(@ARGV)!=1)
   || (!defined($compareListFile) && scalar(@ARGV)!=2))
{
    if ($help)
    {
        pod2usage(-verbose => 2);
    }
    else
    {
        pod2usage(1);
    }
}

if (!defined($compareListFile))
{
	$genotypeFileA = $ARGV[0];
	$genotypeFileB = $ARGV[1];
}
else
{
	$genotypeFileA = $ARGV[0];
	$genotypeFileB = $ARGV[0];
}

#checks if input is genotype file 
if(isGt($genotypeFileA) && isGt($genotypeFileB))
{
	$isGt = 1;
}
elsif (isTg($genotypeFileA) && isTg($genotypeFileB))
{
	$isGt = 0;
}
else
{
	die "Input files should be genotype files in the same orientation.";
}

open(A, $genotypeFileA) || die "Cannot open $genotypeFileA\n";
open(B, $genotypeFileB) || die "Cannot open $genotypeFileB\n";

if(!defined($outFile))
{
    my ($name, $path, $ext) = fileparse($genotypeFileA, '\..*');
	my $fileA = $name;
	($name, $path, $ext) = fileparse($genotypeFileB, '\..*');
	my $fileB = $name;
	
	if (defined($compareListFile))
	{
		$outFile = "$fileA-call-comparison.txt";	
	}
	else
	{
		$outFile = "$fileA-$fileB-call-comparison.txt";	
	}
} 
	
#open genotype files A and B and parses the first header row
$_ = <A>;
$fposA = length($_);
s/\r?\n?$//;
$colNoA = s/\t/\t/g + 1;
@fieldsA = split('\t', $_, $colNoA);

$_ = <B>;
$fposB = length($_);
s/\r?\n?$//;
$colNoB = s/\t/\t/g + 1;
@fieldsB = split('\t', $_, $colNoB);

#populates
#%A_COL_ELEMENTS;
#%A_ROW_ELEMENTS;
#%B_COL_ELEMENTS;
#%B_ROW_ELEMENTS;
if (!defined($compareListFile))
{	
	map {$A_COL_ELEMENTS{$fieldsA[$_]} = $_} (1..$#fieldsA);
	
	while (<A>)
	{
		s/\r?\n?$//;
		@fieldsA = split('\t', $_, 2);
		$A_ROW_ELEMENTS{$fieldsA[0]} = $fposA;
		$fposA += length($_) + 1;
	}
	
	#check for intersection for SNPs and samples from genotype file 2
	#populate common elements and also to keep reference of the positions of the common column elements in %ROW_ELEMENT
	map {$B_COL_ELEMENTS{$fieldsB[$_]}=$_ if exists($A_COL_ELEMENTS{$fieldsB[$_]})} (1..$#fieldsB);
	
	while (<B>)
	{
	    s/\r?\n?$//;
		@fieldsB = split('\t', $_, 2);
		if(exists($A_ROW_ELEMENTS{$fieldsB[0]}))
		{
			$B_ROW_ELEMENTS{$fieldsB[0]} = $fposB;
			push(@COMPARE_A, $fieldsB[0]);
		    push(@COMPARE_B, $fieldsB[0]);
		}
		$fposB += length($_) + 1;
	}
	
	if(scalar(keys(%B_COL_ELEMENTS))==0 || scalar(keys(%B_ROW_ELEMENTS))==0)
	{
		die "There are no common elements to compare in $genotypeFileA annd $genotypeFileB";
	}
}
else
{
	open(COMPARE_LIST, $compareListFile) || die "Cannot open compare list file";
	
	while(<COMPARE_LIST>)
	{
	    s/\r?\n?$//;
	    
	    if($.==1)
	    {
	        $colNo = s/\t/\t/g + 1;
	        
	        my @fields = split('\t', $_, $colNo);
	        	        
	        SEARCH_LABEL: for my $label ('sample-id-a', 'sample-id-b', 'snp-id-a', 'snp-id-b')
	        {
	            for my $col (0 .. $#fields)
	            {
	                if ($fields[$col] eq $label)
	                {           
	                    $label2col{$label}=$col;
	                    next SEARCH_LABEL;
	                }
	            }
	        }
	        
	        my @keys = keys(%label2col);
	        
	        unless ((scalar(@keys)==2 && $isGt && exists($label2col{'sample-id-a'}) && exists($label2col{'sample-id-b'}))
	          ||(scalar(@keys)==2 && !$isGt && exists($label2col{'snp-id-a'}) && exists($label2col{'snp-id-b'})))
	        {
	        	die "Compare list is expected to have either (sample-id-a, sample-id-b) for a gt-file or (snp-id-a, snp-id-b) for a tg-file";
	        }
	    }
	    else
	    {
	    	my @fields = split('\t', $_, $colNo);
	    	
            if ($isGt)
            {
                if (defined($fields[$label2col{'sample-id-a'}]) && defined($fields[$label2col{'sample-id-b'}]))
                {
    		        $COMPARE{$fields[$label2col{'sample-id-a'}]} = $fields[$label2col{'sample-id-b'}];
    		        push(@COMPARE_A, $fields[$label2col{'sample-id-a'}]);
    		        push(@COMPARE_B, $fields[$label2col{'sample-id-b'}]);
			    }
			    else
			    {
			        die "Compare list not correctly defined";
			    }
			}
			else
			{
			    if (defined($fields[$label2col{'sample-id-a'}]) && defined($fields[$label2col{'sample-id-b'}]))
                {
    		        $COMPARE{$fields[$label2col{'snp-id-a'}]} = $fields[$label2col{'snp-id-b'}];
    		        push(@COMPARE_A, $fields[$label2col{'snp-id-a'}]);
    		        push(@COMPARE_B, $fields[$label2col{'snp-id-b'}]);
    		    }
			    else
			    {
			        die "Compare list not correctly defined";
			    }
			}
	    }
	}
	close(COMPARE_LIST);
	
	#search for elements
	map {$A_COL_ELEMENTS{$fieldsA[$_]} = $_} (1..$#fieldsA);
	map {$B_COL_ELEMENTS{$fieldsB[$_]} = $_} (1..$#fieldsB);
		
	while (<A>)
	{
		s/\r?\n?$//;
		@fieldsA = split('\t', $_, 2);
		$A_ROW_ELEMENTS{$fieldsA[0]} = $fposA;
		$fposA += length($_) + 1;
	}
	
	for (my $i=$#COMPARE_A; $i>=0; --$i)
	{
		my $rowElementA = $COMPARE_A[$i];
		my $rowElementB = $COMPARE_B[$i];
			
		if(exists($A_ROW_ELEMENTS{$rowElementA}))
		{		
			if(exists($A_ROW_ELEMENTS{$rowElementB}))
			{
				$B_ROW_ELEMENTS{$rowElementB} = $A_ROW_ELEMENTS{$rowElementB};
			}
			else
			{
				splice(@COMPARE_A, $i, 1);
				splice(@COMPARE_B, $i, 1);
				warn "$rowElementB not in $genotypeFileB, ($rowElementA, $rowElementB) pair discarded";
			}
		}
		else
		{
			splice(@COMPARE_A, $i, 1);
			splice(@COMPARE_B, $i, 1);
			warn "$rowElementA not in $genotypeFileA, ($rowElementA, $rowElementB) pair discarded";
		}
	}
}

#############
#COMMON CODE#
#############
#scan through A and B for common elements.
open(OUT, ">$outFile") || die "Cannot open $outFile\n";
$_ = <A>;
@fieldsA = split('\t', $_, 2);
if (defined($compareListFile))
{
	print OUT "$fieldsA[0]\tduplicate-$fieldsA[0]\tsimilarity\tN\t0A\t1A\t2A\t-1A\t0B\t1B\t2B\t-1B";
}
else
{
	print OUT "$fieldsA[0]\tsimilarity\tN\t0A\t1A\t2A\t-1A\t0B\t1B\t2B\t-1B";	
}
print OUT "\t0->0\t0->1\t0->2\t0->-1";
print OUT "\t1->0\t1->1\t1->2\t1->-1";
print OUT "\t2->0\t2->1\t2->2\t2->-1";
print OUT "\t-1->0\t-1->1\t-1->2\t-1->-1\n";

for my $i (0 .. $#COMPARE_A)
{
	my $rowElementA = $COMPARE_A[$i];
	my $rowElementB = $COMPARE_B[$i];
	
	seek(A, $A_ROW_ELEMENTS{$rowElementA}, 0);
	$_ = <A>;
	s/\r?\n?$//;
	@fieldsA = split('\t', $_, $colNoA);
	

	seek(B, $B_ROW_ELEMENTS{$rowElementB}, 0);
	$_ = <B>;
	s/\r?\n?$//;
	@fieldsB = split('\t', $_, $colNoB);
	
	my %A_GENO;
	map {$A_GENO{$_}=0} (0,1,2,-1);
	my %B_GENO;
	map {$B_GENO{$_}=0} (0,1,2,-1);
	my $similarity = 0;
	my $valid = 0;
	my %TRANSITION;
	map {$TRANSITION{$_}=0} ('00', '01', '02', '0-1', '10', '11', '12', '1-1' , '20', '21', '22', '2-1', '-10', '-11', '-12', '-1-1');

	for my $colElement (keys(%B_COL_ELEMENTS))
	{
		my $aGeno = $fieldsA[$A_COL_ELEMENTS{$colElement}];
		my $bGeno = $fieldsB[$B_COL_ELEMENTS{$colElement}];
		
		if ($aGeno!=-1 && $bGeno !=-1)
		{
			$valid++;
			
			if ($aGeno==$bGeno)
			{
				$similarity++;
			}	
		}
		
		$TRANSITION{"$aGeno$bGeno"}++;
		
		$A_GENO{$aGeno}++;
		$B_GENO{$bGeno}++;
	}

	if ($valid!=0)
	{
		printf " $similarity/$valid = %f - total: %d\n", $similarity/$valid, scalar(keys(%A_COL_ELEMENTS));

		if (defined($compareListFile))
		{	
			printf OUT "$rowElementA\t$rowElementB\t%f\t%d\t$A_GENO{0}\t$A_GENO{1}\t$A_GENO{2}\t$A_GENO{-1}\t$B_GENO{0}\t$B_GENO{1}\t$B_GENO{2}\t$B_GENO{-1}", $similarity/$valid, $valid;
		}
		else
		{
			printf OUT "$rowElementA\t%f\t%d\t$A_GENO{0}\t$A_GENO{1}\t$A_GENO{2}\t$A_GENO{-1}\t$B_GENO{0}\t$B_GENO{1}\t$B_GENO{2}\t$B_GENO{-1}", $similarity/$valid, $valid;			
		}
		print OUT "\t$TRANSITION{'00'}\t$TRANSITION{'01'}\t$TRANSITION{'02'}\t$TRANSITION{'0-1'}";
		print OUT "\t$TRANSITION{'10'}\t$TRANSITION{'11'}\t$TRANSITION{'12'}\t$TRANSITION{'1-1'}";
		print OUT "\t$TRANSITION{'20'}\t$TRANSITION{'21'}\t$TRANSITION{'22'}\t$TRANSITION{'2-1'}";
		print OUT "\t$TRANSITION{'-10'}\t$TRANSITION{'-11'}\t$TRANSITION{'-12'}\t$TRANSITION{'-1-1'}\n";
	}
	else
	{
		printf " n/a, no known genotype to compare\n";
	
		if (defined($compareListFile))
		{
			printf OUT "$rowElementA\t$rowElementB\tn/a\t%d\t$A_GENO{0}\t$A_GENO{1}\t$A_GENO{2}\t$A_GENO{-1}\t$B_GENO{0}\t$B_GENO{1}\t$B_GENO{2}\t$B_GENO{-1}", $valid;
		}
		else
		{
			printf OUT "$rowElementA\tn/a\t%d\t$A_GENO{0}\t$A_GENO{1}\t$A_GENO{2}\t$A_GENO{-1}\t$B_GENO{0}\t$B_GENO{1}\t$B_GENO{2}\t$B_GENO{-1}", $valid;			
		}
		print OUT "\t$TRANSITION{'00'}\t$TRANSITION{'01'}\t$TRANSITION{'02'}\t$TRANSITION{'0-1'}";
		print OUT "\t$TRANSITION{'10'}\t$TRANSITION{'11'}\t$TRANSITION{'12'}\t$TRANSITION{'1-1'}";
		print OUT "\t$TRANSITION{'20'}\t$TRANSITION{'21'}\t$TRANSITION{'22'}\t$TRANSITION{'2-1'}";
		print OUT "\t$TRANSITION{'-10'}\t$TRANSITION{'-11'}\t$TRANSITION{'-12'}\t$TRANSITION{'-1-1'}\n";
	}
}

close(A);
close(B);
close(OUT);