<commit_before>#!/usr/bin/perl
use strict;
use Getopt::Long;

my $gap = 5000000;
my $file = shift(@ARGV);
my $help = 0;

GetOptions("gap=i" => \$gap, "help" => \$help);

chomp $file;
my $uniqfile = $file;
$uniqfile =~ s/\.bam/.rmdup.bam/;

if($help || $file eq "")
{
	print "Usage: uniqsample.pl [options] infile\ninfile: Sorted file of reads\nOptions:\n --gap n\tCalculate unique read target counts at multiples of n [default 5000000]\n"; 
	exit 1;
}


#Samtools must be in system path
open(IN, "samtools view -h $file | ") or die "$!";   
open(UNIQ, " | samtools view -bSo $uniqfile -") or die "$!";

my ($chr, $start, $cigar, $strand, $oldline); 

my $header;
my $length = 0;
my %counts;
my %startstrandcounts;
my $lc = 0;
my $uniq = 0;
my $target;
my ($l, $ct, $op);
my @s;
my @numop;
my @cigararray;

while(<IN>)
{
	$l = $_;

	#Store header
	if($l =~ /^@/)
	{
		$header .= $l;
		print UNIQ "$l";
	}


	else
	{
		@s = split "\t", $l;

		$strand = $s[1] & 16;

		@cigararray = split(/(?<=[A-Z])/, $s[5]);

		$length = 0;

		if(!($s[1] & 1024))
		{
			print UNIQ "$l";
		}	

		if($cigararray[0] =~ /S$/)   #Remove soft clip length for read length
		{
			my @nummasked = split(/(?<=[A-Z])/, $cigararray[0]);
			$s[3] -= $nummasked[0];
		}
		

		foreach $op (@cigararray)
		{
			if ($op =~ /[MDNXP]$/)  #Calculate read length from CIGAR string
			{
				my @numop = split(/(?<=[A-Z])/, $op);
				$length += $numop[0];
			}
				
		}
		
		#Count how many reads map to identical chromosome, start, length (end) and strand
		#if($s[3] == $start && $s[2] eq $chr)
		if($s[2] eq $chr)
		{
			$startstrandcounts{"$s[3]\t$length\t$strand"}++;
		}

		else
		{
			foreach $ct (values %startstrandcounts)
			{
				$counts{$ct}++;
			}
			
			%startstrandcounts = ();
			$chr = $s[2];
			#$start = $s[3];
		}
	}
}

#Increment number of reads with $ct duplicates
foreach $ct (values %startstrandcounts)
	{
		$counts{$ct}++;
	}
	
%startstrandcounts = ();



close(IN);
close(UNIQ);

#Create R command for calculating total number of reads 
# $uniq = current number of unique reads
# $target = desired number of unique reads
# Solve $target - $uniq + sum(n * $counts{n}) = 0 

my $Rcomm = "uniroot(function(x){";

foreach my $dupct (sort {$a <=> $b} keys %counts)
{
	if($dupct)
	{
		$uniq += $counts{$dupct};
		$lc += $dupct * $counts{$dupct};
		$Rcomm .= "$counts{$dupct}*x^$dupct+";
	}
}

$Rcomm .= "-$uniq";
my @filehandles;
my @targets;
my @probs;

print "$target\t$gap\t$uniq\t".$gap * int($uniq/$gap)."\n";

#For each target read count, write to a file, run and parse R command, and open filehandles for bams<commit_msg>Add larger protecting stack in R to avoid error<commit_after>#!/usr/bin/perl
use strict;
use Getopt::Long;

my $gap = 5000000;
my $file = shift(@ARGV);
my $help = 0;

GetOptions("gap=i" => \$gap, "help" => \$help);

chomp $file;
my $uniqfile = $file;
$uniqfile =~ s/\.bam/.rmdup.bam/;

if($help || $file eq "")
{
	print "Usage: uniqsample.pl [options] infile\ninfile: Sorted file of reads\nOptions:\n --gap n\tCalculate unique read target counts at multiples of n [default 5000000]\n"; 
	exit 1;
}


#Samtools must be in system path
open(IN, "samtools view -h $file | ") or die "$!";   
open(UNIQ, " | samtools view -bSo $uniqfile -") or die "$!";

my ($chr, $start, $cigar, $strand, $oldline); 

my $header;
my $length = 0;
my %counts;
my %startstrandcounts;
my $lc = 0;
my $uniq = 0;
my $target;
my ($l, $ct, $op);
my @s;
my @numop;
my @cigararray;

while(<IN>)
{
	$l = $_;

	#Store header
	if($l =~ /^@/)
	{
		$header .= $l;
		print UNIQ "$l";
	}


	else
	{
		@s = split "\t", $l;

		$strand = $s[1] & 16;

		@cigararray = split(/(?<=[A-Z])/, $s[5]);

		$length = 0;

		if(!($s[1] & 1024))
		{
			print UNIQ "$l";
		}	

		if($cigararray[0] =~ /S$/)   #Remove soft clip length for read length
		{
			my @nummasked = split(/(?<=[A-Z])/, $cigararray[0]);
			$s[3] -= $nummasked[0];
		}
		

		foreach $op (@cigararray)
		{
			if ($op =~ /[MDNXP]$/)  #Calculate read length from CIGAR string
			{
				my @numop = split(/(?<=[A-Z])/, $op);
				$length += $numop[0];
			}
				
		}
		
		#Count how many reads map to identical chromosome, start, length (end) and strand
		#if($s[3] == $start && $s[2] eq $chr)
		if($s[2] eq $chr)
		{
			$startstrandcounts{"$s[3]\t$length\t$strand"}++;
		}

		else
		{
			foreach $ct (values %startstrandcounts)
			{
				$counts{$ct}++;
			}
			
			%startstrandcounts = ();
			$chr = $s[2];
			#$start = $s[3];
		}
	}
}

#Increment number of reads with $ct duplicates
foreach $ct (values %startstrandcounts)
	{
		$counts{$ct}++;
	}
	
%startstrandcounts = ();



close(IN);
close(UNIQ);

#Create R command for calculating total number of reads 
# $uniq = current number of unique reads
# $target = desired number of unique reads
# Solve $target - $uniq + sum(n * $counts{n}) = 0 

my $Rcomm = "options(expressions=10000);
uniroot(function(x){";

foreach my $dupct (sort {$a <=> $b} keys %counts)
{
	if($dupct)
	{
		$uniq += $counts{$dupct};
		$lc += $dupct * $counts{$dupct};
		$Rcomm .= "$counts{$dupct}*x^$dupct+";
	}
}

$Rcomm .= "-$uniq";
my @filehandles;
my @targets;
my @probs;

print "$target\t$gap\t$uniq\t".$gap * int($uniq/$gap)."\n";

#For each target read count, write to a file, run and parse R command, and open filehandles for bams
{'fraction_non_alphanumeric': 0.20749716016660355, 'fraction_numerical': 0.016471033699356305, 'mean_word_length': 3.309135399673736, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1255250', 'n_tokens_mistral': 2513, 'n_tokens_neox': 2200, 'n_words': 586}