use List::Util qw(shuffle);
$|=1;

$genome_file="C:/PublicFolder/Software/Bowtie/Hsap.genome";
$phasi_out="phasiRNA_test_dataset.fas";
@background=("C:/PublicFolder/Software/SRAtoolkit/SRR950876.fasta.clipped","C:/PublicFolder/Software/SRAtoolkit/SRR950878.fasta.clipped");
@background_reads=(100000,500000,1000000,2000000,3000000,4000000,5000000);

open(GENOME,$genome_file)||die print $!;
while(<GENOME>)
	{
	$_=~s/\s*$//;
	if($_=~/^>/)
		{
		$head=$_;
		print".";
		}
	else
		{
		$genome{$head}.=$_;
		}
	}
close GENOME;

### CREATE FASTA FILE WITH ARTIFICIAL PHASED RNA SEQUENCES ###
open(TESTSET,">$phasi_out")||die print $!;
foreach$chr(keys%genome)
	{
	next if($chr!~/^>Chr/);
	print"\n$chr";
	
	# Settings
	$phasi_size=21;
	$phasi_n=100; # must be an even integer
	$start=1000001;
	$increment=10000;
	$reads_per_phasiRNA=1;

	$reduce=-1;
	$start=$start-$increment;
	while(1)
		{
		@phasi=();
		$start+=$increment;
		$end=$start+($phasi_size*($phasi_n/2));
		$stretch=substr($genome{$chr},$start,($phasi_size+2)*($phasi_n+1));
		next if($stretch=~/N/);
		print" $start";
		
		# make + phasiRs
		$pos=2-$phasi_size;
		foreach$i(1..$phasi_n/2)
			{
			$pos+=$phasi_size;
			push(@phasi,substr($stretch,$pos,$phasi_size));
			}
		
		# make - phasiRs
		$pos=-$phasi_size;
		foreach$i(1..$phasi_n/2)
			{
			$pos+=$phasi_size;
			$pre=substr($stretch,$pos,$phasi_size);
			$pre=reverse$pre;
			$pre=~tr/ATGC/TACG/;
			push(@phasi,$pre);
			}
		
		# randomly reject phasiRs
		@phasi=shuffle@phasi;
		$reduce++;
		$phasi_count=$phasi_n-$reduce;
		foreach(1..$reduce)
			{
			shift@phasi;
			}
		
		$phasiR_seqid=0;
		foreach$seq(@phasi)
			{
			$phasiR_seqid++;
			$phasiR_readid=0;
			foreach(1..$reads_per_phasiRNA)
				{
				$phasiR_readid++;
				print TESTSET"$chr:$start-$end:$phasiR_seqid:$phasiR_readid\n$seq\n";
				}
			}
		
		last if($phasi_count==10);
		}
	}
close TESTSET;

### CREATE FASTA FILES WITH PHASED RNA SEQUENCES AND BACKGROUND NOISE ###
$background_file=0;
foreach$file(@background)
	{
	$background_file++;
	@noise_seq=();
	open(IN,$file)||die print$!;
	while(<IN>)
		{
		if($_!~/^>/&&$_!~/N/)
			{
			$_=~s/\s*$//;
			push(@noise_seq,$_);
			}
		}
	close IN;
	$n_reads=@noise_seq;
	print"\Loaded $n_reads from $file.";
	@noise_seq=shuffle@noise_seq;
	foreach$n_background(@background_reads)
		{
		print".";
		open(OUT,">$background_file.$n_background")||die print$!;
		
		# add phased sequences
		open(TESTSET,$phasi_out)||die print$!;
		while(<TESTSET>)
			{
			print OUT$_;
			}
		close TESTSET;
		
		# add noise sequences
		foreach(0..$n_background-1)
			{
			print OUT">noise\n$noise_seq[$_]\n";
			}
		close OUT;
		
		# create collapsed fasta file. Needs the Perl script TBr2_collapse.pl from the NGS toolbox
		system("TBr2_collapse.pl -i $background_file.$n_background -o $background_file.$n_background.collapsed");
		}
	}