SEQUEL_stats.pl

	use strict;
	use Getopt::Long;
	
	
	my %PROGRAM;
	$PROGRAM{'name'}	= "SEQUEL_stats.pl";
	$PROGRAM{'version'}	= "260717A";
	
	my $PROGRAM_INFO = "### File generated by: \"".$PROGRAM{'name'}."\" (".$PROGRAM{'version'}.")\n";
	
	
	my %PARAMS;
	
	GetOptions(\%PARAMS, 'Hn=s', 'HpSn=s', 'HpSp=s');
	
	if(!exists($PARAMS{'Hn'})   || !defined($PARAMS{'Hn'}))		{ die("\n'Hn' file not specified!\n\n"); }
	if(!exists($PARAMS{'HpSn'}) || !defined($PARAMS{'HpSn'}))	{ die("\n'HpSn' file not specified!\n\n"); }
	if(!exists($PARAMS{'HpSp'}) || !defined($PARAMS{'HpSp'}))	{ die("\n'HpSp' file not specified!\n\n"); }
	
	
	my @g_hist	= ();
	
	for(my $i=1;$i<=3000;++$i) { push(@g_hist, [0,0]); }	# Histogram has 3000 bins of 100nt width. 300kb should be enough to capture all PacBio SEQUEL reads ...
	
	my %g_LQdata	= ('id' => "", 'count' => 0, 'zr_seq' => 0);
	my %g_MQdata	= ('id' => "", 'count' => 0, 'zr_seq' => 0);
	my %g_HQdata	= ('id' => "", 'count' => 0, 'zr_seq' => 0);
	
	$g_HQdata{'pr'}	= {'rds' => [], 'seq' => 0, 'mean' => 0, 'median' => 0, 'n50' => 0, 'zor' => 0, 'psr' => 0};
	$g_HQdata{'sr'}	= {'rds' => [], 'seq' => 0, 'mean' => 0, 'median' => 0, 'n50' => 0, 'a_cnt' => [], 'l_flg' => []};
	
	
	parse_LQ_data($PARAMS{'Hn'}, \%g_LQdata);				# a.k.a. "Hn"   = "High quality region: negative" (i.e. the basecaller couldn't find a region in the ZMW read it would consider 'proper' sequence)
	parse_MQ_data($PARAMS{'HpSn'}, \%g_MQdata);				# a.k.a. "HpSn" = "High quality region: positive; Subread: negative" (i.e. a good ZMW read interval was found but for some reason the sequence is considered 'bad')
	parse_HQ_data($PARAMS{'HpSp'}, \%g_HQdata, \@g_hist);	# a.k.a. "HpSp" = "High quality region: positive; Subread: positive" (i.e. that's the good stuff ;)
	
	print_data($PARAMS{'Hn'}, $PARAMS{'HpSn'}, $PARAMS{'HpSp'}, \%g_LQdata, \%g_MQdata, \%g_HQdata, \@g_hist);
	
	
	
	
	sub parse_LQ_data
	{
		my $path	= shift;
		my $data	= shift;
		
		
		open(LQ, "<", $path) or die("\nCan't open LQ file:\t\"".$path."\"!\n\n");
		
		while(my $line = <LQ>) {
			chomp($line);
			
			if($data->{'id'} eq "") { $data->{'id'} = (split(/\t/, $line))[0]; }
			
			$data->{'zr_seq'} += (split(/\t/, $line))[2];
			
			++$data->{'count'};
		}
		
		close LQ;
		
		
		return 0;
	}
	
	
	sub parse_MQ_data
	{
		my $path	= shift;
		my $data	= shift;
		
		
		open(MQ, "<", $path) or die("\nCan't open MQ file:\t\"".$path."\"!\n\n");
		
		while(my $line = <MQ>) {
			chomp($line);
			
			if($data->{'id'} eq "") { $data->{'id'} = (split(/\t/, $line))[0]; }
			
			$data->{'zr_seq'} += (split(/\t/, $line))[2];
			
			++$data->{'count'};
		}
		
		close MQ;
		
		
		return 0;
	}
	
	
	sub parse_HQ_data
	{
		my $path	= shift;
		my $data	= shift;
		my $hist	= shift;
		
		
		my @tmp	= ();
		
		
		open(HQ, "<", $path) or die("\nCan't open HQ file:\t\"".$path."\"!\n\n");
		
		while(my $line = <HQ>) {
			chomp($line);
			
			++$data->{'count'};
			
			@tmp	= split(/\t/, $line);
			
			if($data->{'id'} eq "") { $data->{'id'} = $tmp[0]; }
			
			$data->{'zr_seq'}		+= $tmp[2];
			$data->{'pr'}->{'seq'}	+= $tmp[4];
			$data->{'sr'}->{'seq'}	+= $tmp[6];
			
			push(@{$data->{'pr'}->{'rds'}}, $tmp[4]);
			push(@{$data->{'sr'}->{'rds'}}, $tmp[6]);
			
			if($tmp[9]  == 1) { push(@{$data->{'sr'}->{'a_cnt'}}, $tmp[6]); }	# Column 10 contains the number of adapters found in the HQ region. If it is "1" it has to be at one side of the longest subread. Hence, the insert must be at least that long
			if($tmp[12] == 1) { push(@{$data->{'sr'}->{'l_flg'}}, $tmp[6]); }	# Column 13 is a binary flag. If it is "1", the longest subread is flanked by two adapters and hence can be used to estimate the library size distribution
			
			#-- Fill histgram bins for 
			++$hist->[ceilling($tmp[4] / 100)]->[0];	# Polymerase reads
			++$hist->[ceilling($tmp[6] / 100)]->[1];	# Subreads
			#--
		}
		
		close HQ;
		
		
		@{$data->{'pr'}->{'rds'}} = sort { $b <=> $a } @{$data->{'pr'}->{'rds'}};
		@{$data->{'sr'}->{'rds'}} = sort { $b <=> $a } @{$data->{'sr'}->{'rds'}};
		
		
		return 0;
	}
	
	
	sub print_data
	{
		my $lqFile	= shift;
		my $mqFile	= shift;
		my $hqFile	= shift;
		
		my $lq		= shift;
		my $mq		= shift;
		my $hq		= shift;
		
		my $hist	= shift;
		
		
		my %files	= ();
		$files{'lq'}->{'stats'}	= $lqFile.".stats";
		$files{'mq'}->{'stats'}	= $mqFile.".stats";
		$files{'hq'}->{'stats'}	= $hqFile.".stats";
		$files{'hq'}->{'a_cnt'}	= $hqFile.".aCnt";
		$files{'hq'}->{'l_flg'}	= $hqFile.".lFlg";
		$files{'hq'}->{'hist'}	= $hqFile.".hist";
		
		
		$hq->{'pr'}->{'zor'}	= $hq->{'count'} / ($lq->{'count'} + $mq->{'count'} + $hq->{'count'});
		$hq->{'pr'}->{'psr'}	= $hq->{'sr'}->{'seq'} / $hq->{'pr'}->{'seq'};
		
		$hq->{'pr'}->{'mean'}	= calc_mean($hq->{'pr'}->{'rds'});
		$hq->{'pr'}->{'median'}	= calc_median($hq->{'pr'}->{'rds'});
		$hq->{'pr'}->{'n50'}	= calc_n50($hq->{'pr'}->{'rds'}, $hq->{'pr'}->{'seq'});
		
		$hq->{'sr'}->{'mean'}	= calc_mean($hq->{'sr'}->{'rds'});
		$hq->{'sr'}->{'median'}	= calc_median($hq->{'sr'}->{'rds'});
		$hq->{'sr'}->{'n50'}	= calc_n50($hq->{'sr'}->{'rds'}, $hq->{'sr'}->{'seq'});
		
		
		open(LQ, ">", $files{'lq'}->{'stats'}) or die("\nUnable to generate LQ stats file:\t\"".$files{'lq'}->{'stats'}."\"!\n\n");
		
		print LQ $lq->{'id'}."\tZR\tCOUNT\t".$lq->{'count'}."\n";
		print LQ $lq->{'id'}."\tZR\tSEQ\t".$lq->{'zr_seq'}."\n";
		
		close LQ;
		
		
		open(MQ, ">",$files{'mq'}->{'stats'}) or die("\nUnable to generate MQ stats file:\t\"".$files{'mq'}->{'stats'}."\"!\n\n");
		
		print MQ $mq->{'id'}."\tZR\tCOUNT\t".$mq->{'count'}."\n";
		print MQ $mq->{'id'}."\tZR\tSEQ\t".$mq->{'zr_seq'}."\n";
		
		close MQ;
		
		
		open(HQ, ">", $files{'hq'}->{'stats'}) or die("\nUnable to generate HQ stats file:\t\"".$files{'hq'}->{'stats'}."\"!\n\n");
		
		print HQ $hq->{'id'}."\tZR\tCOUNT\t".$hq->{'count'}."\n";
		print HQ $hq->{'id'}."\tZR\tSEQ\t".$hq->{'zr_seq'}."\n";
		print HQ $hq->{'id'}."\tPR\tSEQ\t".$hq->{'pr'}->{'seq'}."\n";
		print HQ $hq->{'id'}."\tSR\tSEQ\t".$hq->{'sr'}->{'seq'}."\n";
		print HQ $hq->{'id'}."\tPR\tMEAN\t".$hq->{'pr'}->{'mean'}."\n";
		print HQ $hq->{'id'}."\tSR\tMEAN\t".$hq->{'sr'}->{'mean'}."\n";
		print HQ $hq->{'id'}."\tPR\tMEDIAN\t".$hq->{'pr'}->{'median'}."\n";
		print HQ $hq->{'id'}."\tSR\tMEDIAN\t".$hq->{'sr'}->{'median'}."\n";
		print HQ $hq->{'id'}."\tPR\tN50\t".$hq->{'pr'}->{'n50'}."\n";
		print HQ $hq->{'id'}."\tSR\tN50\t".$hq->{'sr'}->{'n50'}."\n";
		print HQ $hq->{'id'}."\tPR\tZOR\t".$hq->{'pr'}->{'zor'}."\n";
		print HQ $hq->{'id'}."\tPR\tPSR\t".$hq->{'pr'}->{'psr'}."\n";
		
		close HQ;
		
		
		open(HQRLH, ">", $files{'hq'}->{'hist'}) or die("\nUnable to generate HQ read length histogram file:\t\"".$files{'hq'}->{'hist'}."\"!\n\n");
		for(my $i=0;$i<scalar(@$hist);++$i) { print HQRLH $hq->{'id'}."\t".($i*100)."\t".$hist->[$i]->[0]."\t".$hist->[$i]->[1]."\n"; }
		close HQRLH;
		
		
		open(HQA, ">", $files{'hq'}->{'a_cnt'}) or die("\nUnable to generate HQ 'aCnt' file:\t\"".$files{'hq'}->{'a_cnt'}."\"!\n\n");
		foreach my $entry (@{$hq->{'sr'}->{'a_cnt'}}) { print HQA $hq->{'id'}."\t".$entry."\n"; }
		close HQA;
		
		open(HQF, ">", $files{'hq'}->{'l_flg'}) or die("\nUnable to generate HQ 'lFlg' file:\t\"".$files{'hq'}->{'l_flg'}."\"!\n\n");
		foreach my $entry (@{$hq->{'sr'}->{'l_flg'}}) { print HQF $hq->{'id'}."\t".$entry."\n"; }
		close HQF;
		
		
		return 0;
	}
	
	
	sub ceilling
	{
		my $value	= shift;
		
		
		if($value == int($value)) { return $value; }
		else { return (int($value) + 1); }
	}
	
	
	sub calc_mean
	{
		my $array	= shift;
		
		
		my $mean = 0;
		
		foreach my $entry (@$array) { $mean += $entry; }
		$mean /= scalar(@$array);
		
		
		return $mean;
	}
	
	
	sub calc_median
	{
		my $array	= shift;
		
		
		my $median = 0;
		
		if(scalar(@$array) % 2)	{ $median = $array->[(scalar(@$array) - 1) / 2]; }
		else					{ $median = ($array->[scalar(@$array) / 2] + $array->[(scalar(@$array) / 2) - 1]) / 2; }
		
		
		return $median;
	}
	
	
	sub calc_n50
	{
		my $array	= shift;
		my $total	= shift;
		
		
		my $idx	= 0;
		my $n50 = 0;
		
		while(($n50 + $array->[$idx]) < ($total / 2)) { $n50 += $array->[$idx]; ++$idx; }
		
		$n50	= $array->[$idx];
		
		
		return $n50;
	}
	
	
	1;