
########   extract fasta sequence files from fna directory based on gene ids of gbk directory. Output file: ./out/...
use DBI;
$dbh = DBI->connect('DBI:mysql:metazoa_1001', 'root', '500068'
	           ) || die "Could not connect to database: $DBI::errstr";




system "rm  -rf ./out_8_taxon_species_remove_unwanted";
system "mkdir ./out_8_taxon_species_remove_unwanted";

$dir = "./gbk";
@files = <$dir/*>;

foreach $gbk_file (@files)
{

	open TXT1 , "<$gbk_file" or die "can't open file1:$!";;


	$ref_seq_name = "";



	while(<TXT1>)
	{
		chomp;	

#		if($_=~ /^LOCUS\s+(\S+)\s+.*/)

		if($_=~ /.*VERSION\s+.*\s+GI\:(\d+)/)
		{         
			$ref_seq_name = $1;
			$sth = $dbh->prepare("SELECT * FROM `cds_rRNA_target_accession_2` WHERE `accession_num` = '$acc_name'");
			$sth->execute();

			while (@row = $sth->fetchrow_array) {


				if($row[2] =~ /(\S+)/)
				{     
					$rna_region = $1;
					@seq_range_array = ();
					$complement = 0;

		#			if($rna_region =~ /(join\()?(complement\()?(\<?(\d+)\.\.\>?(\d+).*)\)?/)
					if($rna_region =~ /(?:join)?(?:complement)?\(?(\<?\d+.*\d+)\)?/)
					{				
						$seq_range_rna = $1;
						@seq_range_array = grep { /\S/ } split(/[\.\.|,|\<|\>]/, $seq_range_rna);								
					}
			
					if($rna_region =~ /^complement.*/)
					{
						$complement = 1;						

					}		
					    	
#					$r = 1;	

					&rna_fasta($row[5],$ref_seq_name);


				}						
			
	
			}
			
		}

		if($_=~/^\s+ACCESSION\s+(\S+).*/)
		{
			$acc_name = $1;




		}


	
	}

	close TXT1;

}


###  extract rRNA sequence
sub rna_fasta{
	my($out_file,$ref_seq_name_file) = @_;

print "$ref_seq_name_file\n";
	open OUT, ">>./out_8_taxon_species_remove_unwanted/$out_file.txt";
#	open TXT3, "<./fna/$ref_seq_name_file\.fna" or die "can't open fna file:$!";
	open TXT3, "<./fna_8_taxon_species_remove_unwanted/$ref_seq_name_file\.fna" or last;
	$s = 0;
	@rna_seq_array = ();
	$rna_seq = "";

	while(<TXT3>)
	{
		chomp;
		### print out fasta file title
		if($s == 0)
		{
			if($_ =~ /(\>\d+.*)/)
			{
				$fasta_title_1 = $1;
#				$fasta_title_2 = $2;
				print OUT "$fasta_title_1\n";
			
				if($complement == 0)
				{
					$c = 0;
					$fasta_range_array = "";
					@range_array = ();
					foreach $range (@seq_range_array)
					{
						if($c == 0)
						{
							$range_1 = $range;
							$c = 1;
							next;

						}

						if($c == 1)
						{
							$range_2  = $range;
							push @range_array, "$range_1\-$range_2";
							$c = 0;
							
						}													

					}
					
					$fasta_range_array = join ",", @range_array;
				
#					print OUT "$fasta_range_array$fasta_title_2\n";


				}


				if($complement == 1)
				{
					$c = 0;
					$fasta_range_array = "";
					@range_array = ();
					foreach $range (@seq_range_array)
					{
						if($c == 0)
						{
							$range_1 = $range;
							$c = 1;
							next;

						}

						if($c == 1)
						{
							$range_2  = $range;
							push @range_array, "c$range_2\-$range_1";
							$c = 0;
							
						}													

					}
					
					$fasta_range_array = join ",", @range_array;
				
#					print OUT "$fasta_range_array$fasta_title_2\n";					

				}

			}	

			
		}

		unless($s == 0)
		{
			$rna_seq .= $_; 
		}
		
		
		$s++;
	}
	close TXT3;
	### whole RNA sequence		
	@rna_seq_array = split //, $rna_seq;
	$u = 0;
	$a = 1;
	$b = 1;
	$comp_seq = "";

	foreach $seq_range_rna (@seq_range_array)
	{
		if($u == 0)
		{
			$seq_range_rna_1 = $seq_range_rna;
			$u = 1;	
			next;
		}

		if($u == 1 )
		{
			$seq_range_rna_2 = $seq_range_rna;
			$u = 2;

		}	


		if($u == 2)
		{
			if($complement == 0)
			{
				for($t = $seq_range_rna_1-1;$t<$seq_range_rna_2;$t++)
				{
					print OUT "$rna_seq_array[$t]";	
					if($a % 70 == 0)
					{
						print OUT "\n";

					}

					$a++;

				}

			}

			if($complement == 1)
			{  
				for($t = $seq_range_rna_1-1;$t<$seq_range_rna_2;$t++)
				{
					$comp_seq .= "$rna_seq_array[$t]";						

				}
				
				$reverse_comp_seq = reverse($comp_seq);
				$reverse_comp_seq =~ tr/ACGTRYKMBVDHacgtrykmbvdh/TGCAYRMKVBHDtgcayrmkvbhd/;
				@comp_seq_array = split //, $reverse_comp_seq;
				$length_comp_array = scalar(@comp_seq_array);
				for($t = 0;$t<$length_comp_array;$t++)
				{
					print OUT "$comp_seq_array[$t]";
					
					if($b % 70 == 0)
					{
						print OUT "\n";

					}

					$b++;
				}					

			}

			$u = 0;
		}


	}
	print OUT "\n";

@seq_range_array = ();
$seq_range_rna_1 = "";		
$seq_range_rna_2 = "";
$acc_name = "";

}



