N stat fasta.pl

From RAJ INFO
Jump to navigation Jump to search


#! /usr/bin/perl -w
# *************************************************************
# This Programme does the following in sequence
# -> Parses a multi fasta
# -> Makes the separate file with completely 'N', >80 Ns and <80 Ns
# -> Written as individual files.
# Author : Rajkumar (itc@rajkumar.in)
# Release: JUL 2009
# *************************************************************

# Libraries
	use strict;
	use Bio::Perl;
	use Bio::SeqIO;

# Variables
	my $path        = "/home/raj/bio";
	my $path1       = "$path/user/raja";
	my $path2       = "$path1/blast/BES-NR";
	my $multifasta  = "$path2/BES.fa";
	
	my $outfile1    = "$path2/BES_All_N.fa";
	my $outfile2    = "$path2/BES_MoreThan_80N.fa";
	my $outfile3    = "$path2/BES_LessThan_81N.fa";
	my $outfile4    = "$path2/BES_withNo_N.fa";
	my $outfile5    = "$path2/BES_Atlest_1N.fa";
	my $stat        = "$path2/stat.txt";
	
# Getting the Grip of folders
	if (-e $outfile1)   {`rm $outfile1`;}
	if (-e $outfile2)   {`rm $outfile2`;}
	if (-e $outfile3)   {`rm $outfile3`;}
	if (-e $outfile4)   {`rm $outfile4`;}
	if (-e $outfile5)   {`rm $outfile5`;}

# making several Fasta files of GDNA
	my $in  = Bio::SeqIO->new(-file => "$multifasta", -format => 'Fasta');
	my $identifier;
	my $sequence;
	
	while ( my $seq = $in->next_seq() ) {
		$identifier = $seq->id;
  		$sequence = $seq->seq;
  		
  		#Trim the identifier
  		#$identifier =~ s/_Size:[0-9]+bp//g; 
		
		my ($fasta) = $_;
		$fasta = ">$identifier\n$sequence\n";

		my $count_seq = length($sequence);
		my $count_N = ($sequence =~ tr/N//);
  		
  		
		if ($count_seq eq $count_N) {
		
			open (FASTA, ">>".$outfile1);
			print FASTA $fasta;
			close FASTA;		
		}  		
		
		
		if ($count_N >80) {
		
			open (FASTA, ">>".$outfile2);
			print FASTA $fasta;
			close FASTA;		
		}  		
  		
		if ($count_N <=80) {
		
			open (FASTA, ">>".$outfile3);
			print FASTA $fasta;
			close FASTA;		
		}  		

		if ($count_N <1) {
		
			open (FASTA, ">>".$outfile4);
			print FASTA $fasta;
			close FASTA;		
		}

		if ($count_N >=1) {
		
			open (FASTA, ">>".$outfile5);
			print FASTA $fasta;
			close FASTA;		
		}  				  		

	}



# Generating STAT
	
	my $c_outfile1 = `grep -c ">" $outfile1`; chomp $c_outfile1;
	my $c_outfile2 = `grep -c ">" $outfile2`; chomp $c_outfile2;	
	my $c_outfile3 = `grep -c ">" $outfile3`; chomp $c_outfile3;
	my $c_outfile4 = `grep -c ">" $outfile4`; chomp $c_outfile4;
	my $c_outfile5 = `grep -c ">" $outfile5`; chomp $c_outfile5;		
	my $total      = `grep -c ">" $multifasta`; chomp $total;	
	

	my $stat_content = "Total sequences:\t$total\nBES with No 'N':\t$c_outfile4\nBES with atleast one 'N':\t$c_outfile5\nBES with all 'N':\t$c_outfile1\nBES with more than 80 'N's:\t$c_outfile2\nBES with less than 81 'N's:\t$c_outfile3\n";
	open (F, ">$stat");
	print F $stat_content;
	close F;		

# END