$|=1;
$species="Homo sapiens";
$miRBase_mature_FASTA="miRBase_mature";
$miRBase_hairpin_FASTA="miRBase_hairpin";

$unmodified=1000;
$abundance_1=300;
$abundance_2=200;
$abundance_3=100;

$freq_A=0.15;
$freq_AA=0.1;
$freq_U=0.15;
$freq_UU=0.1;

$internal_mod_freq=0.1;

%mod=
	(
	"A"=>"T",
	"T"=>"A",
	"G"=>"C",
	"C"=>"G",
	);

open(MATURE,$miRBase_mature_FASTA)||die print$!;
while(<MATURE>)
	{
	$_=~s/\s*$//;
	if($_=~s/^>//)
		{
		if($_=~/$species/)
			{
			$hsa=1;
			$_=~s/ .+$//;
			$head=$_;
			}
		else
			{
			$hsa=0;
			}
		}
	elsif($hsa==1)
		{
		$_=~s/U/T/g;
		$mature{$_}=$head;
		}
	}
close MATURE;

open(HAIRPIN,$miRBase_hairpin_FASTA)||die print$!;
while(<HAIRPIN>)
	{
	$_=~s/\s*$//;
	if($_=~s/^>//)
		{
		if($_=~/$species/)
			{
			$hsa=1;
			$_=~s/ .+$//;
			$head=$_;
			}
		else
			{
			$hsa=0;
			}
		}
	elsif($hsa==1)
		{
		$_=~s/U/T/g;
		$hairpin{$head}.=$_;
		}
	}
close HAIRPIN;

$species_output=$species;
$species_output=~s/ /_/;
open(FASTA,">test_miR_dataset-$species_output.fas");
open(SUMMARY,">test_miR_dataset.$species_output.txt");
foreach$miRseq(keys%mature)
	{
	foreach$prec_name(keys%hairpin)
		{
		if($hairpin{$prec_name}=~/$miRseq/) 
			{
			$freq=$unmodified;$seq=$&;$info="no_offset";make_seqs(); # original, mature miR
			$freq=$abundance_1;$seq=substr($&,0,(length$&)-1);$info="-1nt_3p";make_seqs();
			$freq=$abundance_2;$seq=substr($&,0,(length$&)-2);$info="-2nt_3p";make_seqs();
			$freq=$abundance_3;$seq=substr($&,0,(length$&)-3);$info="-3nt_3p";make_seqs();
			$freq=$abundance_1;$seq=substr($&,1);$info="-1nt_5p";make_seqs();
			$freq=$abundance_2;$seq=substr($&,2);$info="-2nt_5p";make_seqs();
			$freq=$abundance_3;$seq=substr($&,3);$info="-3nt_5p";make_seqs();
			$freq=$abundance_2;$seq=substr($&,1,(length$&)-2);$info="-1nt_5p, -1nt_3p";make_seqs();
			$freq=$abundance_3;$seq=substr($&,2,(length$&)-4);$info="-2nt_5p, -2nt_3p";make_seqs();
			}
		if($hairpin{$prec_name}=~/$miRseq./)
			{
			$freq=$abundance_1;$seq=$&;$info="+1nt_3p";make_seqs();
			$freq=$abundance_2;$seq=substr($&,1);$info="-1nt_5p, +1nt_3p";make_seqs();
			}
		if($hairpin{$prec_name}=~/$miRseq../)
			{
			$freq=$abundance_2;$seq=$&;$info="+2nt_3p";make_seqs();
			$freq=$abundance_3;$seq=substr($&,2);$info="-2nt_5p, +2nt_3p";make_seqs();
			}
		if($hairpin{$prec_name}=~/$miRseq.../)
			{
			$freq=$abundance_3;$seq=$&;$info="+3nt_3p";make_seqs();
			}
		if($hairpin{$prec_name}=~/.$miRseq/)
			{
			$freq=$abundance_1;$seq=$&;$info="+1nt_5p";make_seqs();
			$freq=$abundance_2;$seq=substr($&,0,(length$&)-1);$info="+1nt_5p, -1nt_3p";make_seqs();
			}
		if($hairpin{$prec_name}=~/..$miRseq/)
			{
			$freq=$abundance_2;$seq=$&;$info="+2nt_5p";make_seqs();
			$freq=$abundance_3;$seq=substr($&,0,(length$&)-2);$info="+2nt_5p, -2nt_3p";make_seqs();
			}
		if($hairpin{$prec_name}=~/...$miRseq/)
			{
			$freq=$abundance_3;$seq=$&;$info="+3nt_5p";make_seqs();
			}
		if($hairpin{$prec_name}=~/.$miRseq./)
			{
			$freq=$abundance_2;$seq=$&;$info="+1nt_5p, +1nt_3p";make_seqs();
			}
		if($hairpin{$prec_name}=~/..$miRseq../)
			{
			$freq=$abundance_3;$seq=$&;$info="+2nt_5p, +2nt_3p";make_seqs();
			}
		}
	}
close FASTA;

foreach$miR_name(keys%summary)
	{
	$total_miRs=0;
	print SUMMARY"$miR_name";
	foreach$mod(sort{$a cmp $b}keys%{$summary{$miR_name}})
		{
		print SUMMARY"\n$mod\t$summary{$miR_name}{$mod}";
		$total_miRs+=$summary{$miR_name}{$mod};
		}
	print SUMMARY"\nTOTAL\t$total_miRs\n\n";
	}
close SUMMARY;

sub make_seqs
	{
	$add_A=0;
	$add_AA=0;
	$add_U=0;
	$add_UU=0;
	$modified=0;
	foreach$i(1..$freq)
		{
		if($add_A<$freq*$freq_A)
			{
			$output_seq=$seq."A";
			$output_info=$info.", +A";
			if($add_A<$freq*$freq_A*$internal_mod_freq)
				{
				@seq=split("",$output_seq);
				$seq[9]=$mod{$seq[9]};
				$output_seq=join("",@seq);
				$output_info=$output_info.", SNP-10";
				}
			$add_A++;
			}
		elsif($add_AA<$freq*$freq_AA)
			{
			$output_seq=$seq."AA";
			$output_info=$info.", +AA";
			if($add_AA<$freq*$freq_AA*$internal_mod_freq)
				{
				@seq=split("",$output_seq);
				$seq[9]=$mod{$seq[9]};
				$output_seq=join("",@seq);
				$output_info=$output_info.", SNP-10";
				}
			$add_AA++;
			}
		elsif($add_U<$freq*$freq_U)
			{
			$output_seq=$seq."T";
			$output_info=$info.", +T";
			if($add_U<$freq*$freq_U*$internal_mod_freq)
				{
				@seq=split("",$output_seq);
				$seq[9]=$mod{$seq[9]};
				$output_seq=join("",@seq);
				$output_info=$output_info.", SNP-10";
				}
			$add_U++;
			}
		elsif($add_UU<$freq*$freq_UU)
			{
			$output_seq=$seq."TT";
			$output_info=$info.", +TT";
			if($add_UU<$freq*$freq_UU*$internal_mod_freq)
				{
				@seq=split("",$output_seq);
				$seq[9]=$mod{$seq[9]};
				$output_seq=join("",@seq);
				$output_info=$output_info.", SNP-10";
				}
			$add_UU++;
			}
		elsif($modified<$freq*$internal_mod_freq)
			{
			@seq=split("",$seq);
			$seq[9]=$mod{$seq[9]};
			$output_seq=join("",@seq);
			$output_info=$info.", SNP-10";
			$modified++;
			}
		else
			{
			$output_info=$info;
			$output_seq=$seq;
			}
		print FASTA">$mature{$miRseq}|$prec_name|$output_info\n$output_seq\n";
		$summary{$mature{$miRseq}}{$output_info}++;
		}
	}