comparison miRPlant.pl @ 46:ca05d68aca13 draft

Uploaded

45:2cb6add23dfe

my $version=1.00;

use strict;
use Getopt::Long;
use threads;
use threads::shared;
use File::Path;
use File::Basename;
use RNA;
use Term::ANSIColor;

my %opts;
GetOptions(\%opts,"i:s@","tag:s@","format=s","gfa=s","pre=s","mat=s","rfam:s","dis:i","flank:i","mfe:f","idx:s","idx2:s","mis:i","r:i","v:i","e:i","f:i","a:s","M:i","t:i","min:i","max:i","o:s","path:s","D","h");
if (!(defined $opts{i} and defined $opts{format} and defined $opts{gfa} and defined $opts{pre} and defined $opts{mat}) || defined $opts{h}) { #necessary arguments
&usage;

print "miPlant program start:\n   The time is $time!\n";
print "Command line:\n   $0 @ARGV\n";

my $format=$opts{'format'};
if ($format ne "fastq" && $format ne "fq" && $format ne "fasta" && $format ne "fa") { 
	&printErr();
	die "Parameter \"-format\" is error! Parameter is fastq, fq, fasta or fa\n";
}

my $phred_qv=64;

	}
}
sub quantify{
	my $tag=join "\\;" ,@mark;
	system("perl $scipt_path/quantify.pl -p $opts{pre} -m $opts{mat} -r $data -o $dir -time $time -mis $mis -t $t -e $upstream -f $downstream -tag $tag");
#	print "\nquantify.pl -p $opts{pre} -m $opts{mat} -r $data -o $dir -time $time -mis $mis -t $t -e $upstream -f $downstream -tag $tag\n";
}
sub filterbylength{
	my $tmpmark=join ",", @mark;
	system("perl $scipt_path/filterReadsByLength.pl -i $collapsed -o $data -min $min_nt -max $max_nt -mark $tmpmark");
	system("perl $scipt_path/Length_Distibution.pl -i $preprocess/reads_length_distribution.txt  -o $preprocess/length.html");

		push @mark,$tmp[1];
		&check_rawdata($tmp[0]);
	}
	close CON;
	if (@filein != @mark) {
		&printErr();
		die "Maybe config file have some wrong!!!\n";
	}
}
sub check_rawdata{
	my ($fileforcheck)=@_;
	if (!(-s $fileforcheck)) {
		&printErr();
		die "Can not find $fileforcheck, or file is empty!!!\n";
	}
	if ($format eq "fasta" || $format eq "fa") {
		&checkfa($fileforcheck);
	}

	my ($file_reads)=@_;
	open N,"<$file_reads";
	my $line=<N>;
	chomp $line;
    if($line !~ /^>\S+/){
        printErr();
        die "The first line of file $file_reads does not start with '>identifier'
Reads file $file_reads is not a valid fasta file\n\n";
    }
    if(<N> !~ /^[ACGTNacgtn]*$/){
        printErr();
        die "File $file_reads contains not allowed characters in sequences
Allowed characters are ACGTN
Reads file $file_reads is not a fasta file\n\n";
    }
	close N;

		chomp $a;
		chomp $b;
		chomp $c;
		chomp $d;
		if($a!~/^\@/){
			&printErr();
			die "$file_reads is not a fastq file\n\n";
		}
		if($b!~ /^[ACGTNacgtn]*$/){
			&printErr();
			die "File $file_reads contains not allowed characters in sequences
Allowed characters are ACGTN
Reads file $file_reads is not a fasta file\n\n";
		}
		if ($c!~/^\@/ && $c!~/^\+/) {
			&printErr();
			die "$file_reads is not a fastq file\n\n";
		}
		if ((length $b) != (length $d)) {
			&printErr();
			die "$file_reads is not a fastq file\n\n";
		}
		my @qv=split //,$d;
		for (my $j=0;$j<@qv ;$j++) {
			my $q=ord($qv[$j])-64;

			push @ret, $file;
		}
	}
	closedir I;
	if (@ret != 1) {
		&printErr();

		die "Can not find directory or file which name has string: $str !!!\n";
	}
	return $ret[0];
}

sub printErr{
    print STDERR color 'bold red';
    print STDERR "Error: ";
    print STDERR color 'reset';
}
=cut
sub Time{
        my $time=time();
        my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];
        $month++;
        $year+=1900;

sub usage{
print <<"USAGE";
Version $version
Usage:
$0 -i -format -gfa -index -pre -mat -rfam -D -a -M -min -max -mis -e -f -v -t  -o  -path
options:
-i string,  input file#input files information file
		/path/filename	mark
		/path/filename	mark
		...

-format string,#specific input rawdata file format : fastq|fq|fasta|fa

-gfa string,  input file # genome fasta. sequence file
-idx string, genome file index, file-prefix #(must be indexed by bowtie-build) The parameter
                string must be the prefix of the bowtie index. For instance, if
                the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then

46:ca05d68aca13

my $version=1.00;

use strict;
use Getopt::Long;
use threads;
#use threads::shared;
use File::Path;
use File::Basename;
#use RNA;
#use Term::ANSIColor;

my %opts;
GetOptions(\%opts,"i:s@","tag:s@","format=s","gfa=s","pre=s","mat=s","rfam:s","dis:i","flank:i","mfe:f","idx:s","idx2:s","mis:i","r:i","v:i","e:i","f:i","a:s","M:i","t:i","min:i","max:i","o:s","path:s","D","h");
if (!(defined $opts{i} and defined $opts{format} and defined $opts{gfa} and defined $opts{pre} and defined $opts{mat}) || defined $opts{h}) { #necessary arguments
&usage;

print "miPlant program start:\n   The time is $time!\n";
print "Command line:\n   $0 @ARGV\n";

my $format=$opts{'format'};
if ($format ne "fastq" && $format ne "fq" && $format ne "fasta" && $format ne "fa") { 
	#&printErr();
	die "Parameter \"-format\" is error! Parameter is fastq, fq, fasta or fa\n";
}

my $phred_qv=64;

	}
}
sub quantify{
	my $tag=join "\\;" ,@mark;
	system("perl $scipt_path/quantify.pl -p $opts{pre} -m $opts{mat} -r $data -o $dir -time $time -mis $mis -t $t -e $upstream -f $downstream -tag $tag");
	print "\nquantify.pl -p $opts{pre} -m $opts{mat} -r $data -o $dir -time $time -mis $mis -t $t -e $upstream -f $downstream -tag $tag\n";
}
sub filterbylength{
	my $tmpmark=join ",", @mark;
	system("perl $scipt_path/filterReadsByLength.pl -i $collapsed -o $data -min $min_nt -max $max_nt -mark $tmpmark");
	system("perl $scipt_path/Length_Distibution.pl -i $preprocess/reads_length_distribution.txt  -o $preprocess/length.html");

		push @mark,$tmp[1];
		&check_rawdata($tmp[0]);
	}
	close CON;
	if (@filein != @mark) {
		#&printErr();
		die "Maybe config file have some wrong!!!\n";
	}
}
sub check_rawdata{
	my ($fileforcheck)=@_;
	if (!(-s $fileforcheck)) {
		#&printErr();
		die "Can not find $fileforcheck, or file is empty!!!\n";
	}
	if ($format eq "fasta" || $format eq "fa") {
		&checkfa($fileforcheck);
	}

	my ($file_reads)=@_;
	open N,"<$file_reads";
	my $line=<N>;
	chomp $line;
    if($line !~ /^>\S+/){
        #printErr();
        die "The first line of file $file_reads does not start with '>identifier'
Reads file $file_reads is not a valid fasta file\n\n";
    }
    if(<N> !~ /^[ACGTNacgtn]*$/){
        #printErr();
        die "File $file_reads contains not allowed characters in sequences
Allowed characters are ACGTN
Reads file $file_reads is not a fasta file\n\n";
    }
	close N;

		chomp $a;
		chomp $b;
		chomp $c;
		chomp $d;
		if($a!~/^\@/){
			#&printErr();
			die "$file_reads is not a fastq file\n\n";
		}
		if($b!~ /^[ACGTNacgtn]*$/){
			#&printErr();
			die "File $file_reads contains not allowed characters in sequences
Allowed characters are ACGTN
Reads file $file_reads is not a fasta file\n\n";
		}
		if ($c!~/^\@/ && $c!~/^\+/) {
			#&printErr();
			die "$file_reads is not a fastq file\n\n";
		}
		if ((length $b) != (length $d)) {
			#&printErr();
			die "$file_reads is not a fastq file\n\n";
		}
		my @qv=split //,$d;
		for (my $j=0;$j<@qv ;$j++) {
			my $q=ord($qv[$j])-64;

			push @ret, $file;
		}
	}
	closedir I;
	if (@ret != 1) {
		#&printErr();

		die "Can not find directory or file which name has string: $str !!!\n";
	}
	return $ret[0];
}

=cut

sub printErr{
    print STDERR color 'bold red';
    print STDERR "Error: ";
    print STDERR color 'reset';
}
sub Time{
        my $time=time();
        my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];
        $month++;
        $year+=1900;

sub usage{
print <<"USAGE";
Version $version
Usage:

$0 -i -format -gfa -index -pre -mat -rfam -D -a -M -min -max -mis -e -f -v -t  -o  -path
options:
-i input files, # raw data file, can be multipe eg. -i xxx.fq -i xxx .fq ...
-tag string # raw data file names, -tag xxx -tag xxx

-format string,#specific input rawdata file format : fastq|fq|fasta|fa

-path scirpt path

-gfa string,  input file # genome fasta. sequence file
-idx string, genome file index, file-prefix #(must be indexed by bowtie-build) The parameter
                string must be the prefix of the bowtie index. For instance, if
                the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then