calc_human_mouse_ratio.pl

#! /usr/bin/env perl 

=head1 NAME

=head1 SYNOPSIS

  perl calc_human_mouse_ratio.pl [-file MatrixMarket file] [-g gene files] [-o prefix output parsed file]  [-mingenes minimum number of gene per cell] [-mintags minimum number of tags per gene] [-h help]

=head1 DESCRIPTION

  This script parses a matrixMarket file generated by indrop pipeline and filter the results according to some params.
  It also make a graph for further clustering suitable for MCL (?)

  
Typical usage is as follows:

  % perl calc_human_mouse_ratio.pl -file input.est.mtx -t -o output parsed  -g input.est.cells.tsv -mingenes 10 -mintags 1
  
=head2 Options

The following options are accepted:

 --file=<mtx name>   	(Mandatory).

 --gene=<tsv file name>   	(Mandatory).

 --mingenes=<symbol>		  			An integer (default 0)

 --mintags=<symbol>		  			An integer (default 0)

 --o=<file>		  		Output prefix (default output)

 --help                   	This documentation.


=head1 AUTHOR

Luca Cozzuto <luca.cozzuto@crg.es> 

=cut
use warnings;
use strict;
use Data::Dumper;
use File::Basename;
use Pod::Usage;
use Getopt::Long;
use List::Util 'sum';
use List::MoreUtils qw(uniq);

my $USAGE = "  perl parseIndropOut.pl [-file] [-o] [-g] [-mingenes] [-mintags] [-h help]
";

my ($input,$output,$genefile,$mingenes,$mintags,$show_help);

&GetOptions(    	
			'file|f=s'		=> \$input,
			'gene|g=s'		=> \$genefile,
			'output|o=s'	=> \$output,
			'mingenes|mg=s'		=> \$mingenes,
			'mintags|mt=s'		=> \$mintags,
			'help|h'        => \$show_help
			)
  or pod2usage(-verbose=>2);
pod2usage(-verbose=>2) if $show_help;

if (!$output) { $output = "output"}
if (!$input) { die ("Please specify input file")}
if (!$genefile) { die ("Please specify tsv gene file")}
if (!$mintags) { $mintags = 0}
if (!$mingenes) { $mingenes = 0};

#READ MTX FILE
my %hash;
my $numrow = 0;
my $totGenes = 0;
my $totCells = 0;
my $outfile	= $output."_parsed.mtx";
my $statfile	= $output."_parsed.tab";
my $resumefile	= $output."_parsed.stats";


open(my $fh, "<", $input)
	or die "Can't open < $input: $!";

open(my $fg, "<", $genefile)
	or die "Can't open < $input: $!";

open(my $of, ">", $outfile)
	or die "Can't open > $outfile: $!";

open(my $orf, ">", $resumefile)
	or die "Can't open > $resumefile: $!";


open(my $osf, ">", $statfile)
	or die "Can't open > $statfile: $!";
	
my @geneDescs;

while (my $row = <$fg>) {
	chomp $row;
	push (@geneDescs, $row);
}


while ( my $row = <$fh>) {
	$numrow++;
	chomp $row;
	if ($numrow==2) {
		my @fields = split(" ", $row);
		$totGenes=$fields[0];
		$totCells=$fields[1];
	}
	elsif ($numrow>2) { 
		my @fields = split(" ", $row);
		my $geneID = $fields[0]-1;
		my $geneNameRaw = $geneDescs[$geneID];
		my $CellID = $fields[1];
		my $tagCount = $fields[2];
		# get genes with tag number higher than mintags
		if ($tagCount >= $mintags) {
			$hash{$CellID}{$geneNameRaw} = $tagCount;
		}		
	}
}
close($fh);

my %statGenes;
my %speciesNames;
	
foreach my $cellID (keys (%hash)) {
	my $numgenes = scalar (keys ($hash{$cellID}));
	# get cells with gene number higher than mingenes
	if ($numgenes >= $mingenes) {
		foreach my $geneNameRaw (keys $hash{$cellID}) {
		    my $tagCount = $hash{$cellID}{$geneNameRaw};
			# write to parsed file
			print $of "cell-".$cellID." ".$geneNameRaw." ".$tagCount."\n";
			my @vals = split("_", $geneNameRaw);
			my $species = $vals[0];
			my $geneName = $vals[1];
			$statGenes{$cellID}{$species}++;
			$speciesNames{$species} = $species;
		}
	}
} 
close($of);
my $header = "cellID";
foreach my $speciesName (sort keys %speciesNames) {
	$header .= "\t".$speciesName;
}
print $osf $header."\n";

my %resume;
$resume{"mixed"} = 0;

foreach my $cellID (keys %statGenes) {
	my $row = "cell-".$cellID;
 	my %species = %{$statGenes{$cellID}};
	my $totCount = sum (values %species);

	$resume{"mixed"}++;

	foreach my $speciesName (sort keys %speciesNames) {
		if (!$resume{$speciesName}) {
			$resume{$speciesName} = 0;
		}
		my $speciesCount = 0;
		if ($species{$speciesName}) {
			$speciesCount = $species{$speciesName};
		}

		$row .= "\t".$speciesCount;

		my $ratio = $speciesCount/$totCount;
		if ($ratio>=0.9) {
			$resume{$speciesName}++;
			$resume{"mixed"}--;
		}
	}
	print $osf $row."\n";
}
close($osf);

foreach my $speciesName (sort keys %resume) {
	print $orf $speciesName."\t".$resume{$speciesName}."\n";

}


close($orf);
exit;