Mercurial > repos > big-tiandm > mirplant2

changeset 22:c0aecae5eb03 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Deleted selected files
author big-tiandm
date Fri, 25 Jul 2014 05:57:20 -0400
parents b963c4a161d2
children 45de5e1ff487
files miRPlant.pl
diffstat 1 files changed, 0 insertions(+), 503 deletions(-) [+]
line wrap: on
line diff
--- a/miRPlant.pl	Fri Jul 25 05:57:15 2014 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,503 +0,0 @@
-#!/usr/bin/perl -w
-#Filename:
-#Author: Tian Dongmei
-#Email: tiandm@big.ac.cn
-#Date: 2014-4-22
-#Modified:
-#Description: plant microRNA prediction
-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;
-}
-
-my $time=&Time();
-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;
-
-
-my @inputfiles=@{$opts{'i'}};
-my @inputtags=@{$opts{'tag'}};
-
-my  $mypath=`pwd`;
-chomp $mypath;
-
-my $dir=defined $opts{'o'} ? $opts{'o'} : "$mypath/miRPlant_out/";
-
-
-unless ($dir=~/\/$/) {$dir.="/";}
-if (not -d $dir) {
-	mkdir $dir;
-}
-my $config=$dir."/input_config";
-open CONFIG,">$config";
-	for (my $i=0;$i<@inputfiles;$i++) {
-		print CONFIG $inputfiles[$i],"\t",$inputtags[$i],"\n";
-	}
-close CONFIG;
-
-my $scipt_path=defined $opts{'path'} ? $opts{'path'} : "/Users/big/galaxy-dist/tools/myTools/";
-
-my $a="ATCTCGTATG";  #adapter
-if (defined $opts{'a'}) {$a=$opts{'a'};}
-
-my $m=6;  #adapter minimum mapped nt
-if (defined $opts{'M'}) {$m=$opts{'M'};}
-
-my $t=1; #threads number
-if (defined $opts{'t'}) {$t=$opts{'t'};}
-
-my $min_nt=19; # minimum reads length
-if (defined $opts{'min'}) {$min_nt=$opts{'min'};}
-
-my $max_nt=28; #maximum reads length
-if (defined $opts{'max'}) {$max_nt=$opts{'max'};}
-
-my $mis=0; #mismatch number for microRNA
-if (defined $opts{'mis'}) {$mis=$opts{'mis'};}
-
-my $mis_rfam=0;# mismatch number for rfam
-if (defined $opts{'v'}) {$mis_rfam=$opts{'v'};}
-
-my $hit=25; # maximum reads mapping hits in genome
-if (defined $opts{'r'}) {$hit=$opts{'r'};}
-
-my $upstream = 2; # microRNA 5' extension
-$upstream = $opts{'e'} if(defined $opts{'e'});
-
-my $downstream = 5;# microRNA 3' extension
-$downstream = $opts{'f'} if(defined $opts{'f'});
-
-my $maxd=defined $opts{'dis'} ? $opts{'dis'} : 200;
-my $flank=defined $opts{'flank'} ? $opts{'flank'} :10;
-my $mfe=defined $opts{'mfe'} ? $opts{'mfe'} : -20;
-
-$time=&Time();
-print "$time, Checking input file!\n";
-
-my (@filein,@mark,@clean);
-#&read_config();
-@filein=@inputfiles;
-@mark=@inputtags;
-
-&checkfa($opts{pre});
-&checkfa($opts{mat});
-&checkfa($opts{gfa});
-
-
-##### clip adpter --> clean data  start
-$time=&Time();
-print "$time, Preprocess:\n   trim adapter, reads collapse and filter reads by length.\n";
-
-$time=~s/:/-/g;
-$time=~s/ /-/g;
-my $preprocess=$dir."preProcess_${time}/";
-mkdir $preprocess;
-my $can_use_threads = eval 'use threads; 1';
-if ($can_use_threads) {
-# Do processing using threads
-	print "Do processing using threads\n";
-	my @filein1=@filein; my @mark1=@mark;
-	while (@filein1>0) {
-		my @thrs; my @res;
-		for (my $i=0;$i<$t ;$i++) {
-			last if(@filein1==0);
-			my $in=shift @filein1;
-			my $out=shift @mark1;
-			push @clean,$preprocess.$out."_clips_adapter.fq";
-			$thrs[$i]=threads->create(\&clips,$in,$out);
-		}
-		for (my $i=0;$i<@thrs;$i++) {
-			$res[$i]=$thrs[$i]->join();
-		}
-	}
-} else {
-# Do not processing using threads
-	print "Do not processing using threads\n";
-	for (my $i=0;$i<@filein ;$i++) {
-		my $in=$filein[$i];
-		my $out=$mark[$i];
-		push @clean,$preprocess.$out."_clips_adapter.fq";
-		&clips($in,$out);
-	}
-}
-
-##### clip adpter --> clean data  end
-
-my $collapsed=$preprocess."collapse_reads.fa";
-my $data=$preprocess."collapse_reads_${min_nt}_${max_nt}.fa"; ## raw clean data
-my $data2;   ### mirbase not mapped reads
-my $data3;   ### rfam not mapped reads
-&collapse(\@clean,$collapsed);   #collapse reads to tags
-
-&filterbylength();  # filter <$min_nt  && >$max_nt
-
-print "The final clean data file is $data, only contains reads which length is among $min_nt\~$max_nt\n\n";
-
-$time=Time();
-print "$time: known microRNA quantify!\n\n";
-
-chdir $dir;
-
-$time=~s/:/-/g;
-$time=~s/ /-/g;
-my $known_result=$dir."miRNA_Express_${time}/";
-&quantify(); ### known microRAN quantify
-
-
-#my $miR_exp_dir=&search($known_result,"miRNA_Express_");
-$data2=$known_result."/mirbase_not_mapped.fa";
-
-my $pathfile="$dir/path.txt";
-open PA,">$pathfile";
-print PA "$config\n";
-print PA "$preprocess\n";
-print PA "$known_result\n";
-
-if (defined $opts{'rfam'}) {              #rfam mapping and analysis
-	$time=Time();
-	print "$time: RNA annotate!\n\n";
-	$time=~s/:/-/g;
-	$time=~s/ /-/g;
-	my $rfam_exp_dir=$dir."rfam_match_${time}";
-	&rfam();
-	#my $rfam_exp_dir=&search($dir,"rfam_match_");
-	$data3=$rfam_exp_dir."/rfam_not_mapped.fa";
-print PA "$rfam_exp_dir\n";
-
-	my $tag=join "\\;" ,@mark;
-	system("perl $scipt_path/count_rfam_express.pl -i $rfam_exp_dir/rfam_mapped.bwt -tag $tag -o rfam_non-miRNA_annotation.txt");
-}
-
-my $data4=$data;
-if (defined $opts{'D'}) {                 #genome mapping 
-	$data4=$data3;
-}else{
-	$data4=$data2;
-}
-
-$time=Time();
-print "$time: Genome alignment!\n\n";
-$time=~s/:/-/g;
-$time=~s/ /-/g;
-my $genome_map=$dir."genome_match_${time}";
-&genome($data4);
-print PA "$genome_map\n";
-#my $genome_map=&search($dir,"genome_match_");
-my $mapfile=$genome_map."/genome_mapped.bwt";
-my $mapfa=$genome_map."/genome_mapped.fa";
-my $unmap=$genome_map."/genome_not_mapped.fa";
-
-#$time=Time();
-#print "$time: Novel microRNA prediction!\n\n";
-
-&predict($mapfa);
-
-close PA;
-system("perl $scipt_path/html.pl -i $pathfile -format $format -o $dir/result.html");
-
-$time=Time();
-print "$time: Program end!!\n";
-
-############################## sub programs ###################################
-sub predict{
-	my ($file)=@_;
-	$time=&Time();
-	print "$time: Novel microRNA prediction!\n\n";
-	$time=~s/:/-/g;
-	$time=~s/ /-/g;
-	my $predict=$dir."miRNA_predict_${time}";
-print PA "$predict\n";
-	mkdir $predict;
-	chdir $predict;
-	system("perl $scipt_path/precursors.pl -map $mapfile -g $opts{gfa} -d $maxd -f $flank -o $predict/excised_precursor.fa -s $predict/excised_precursor_struc.txt -e $mfe");
-#	print "\nprecursors.pl -map $mapfile -g $opts{gfa} -d $maxd -f $flank -o $predict/excised_precursor.fa -s $predict/excised_precursor_struc.txt -e $mfe\n";
-	
-	system("bowtie-build -f excised_precursor.fa excised_precursor");
-#	print "\nbowtie-build -f excised_precursor.fa excised_precursor\n";
-	
-	system("bowtie -v $mis -f -p $t -m $hit -a --best --strata excised_precursor $file > precursor_mapped.bwt");
-#	print "\nbowtie -v $mis -f -p $t -m $hit -a --best --strata excised_precursor $file > precursor_mapped.bwt\n";
-	
-	system("perl $scipt_path/convert_bowtie_to_blast.pl  precursor_mapped.bwt $file excised_precursor.fa > precursor_mapped.bst");
-#	print "\nconvert_bowtie_to_blast.pl  precursor_mapped.bwt $file excised_precursor.fa > precursor_mapped.bst\n";
-
-	system("sort +3 -25  precursor_mapped.bst  > signatures.bst");
-#	print "\nsort +3 -25  precursor_mapped.bst  > ../signatures.bst\n";
-
-	chdir $dir;
-	system("perl $scipt_path/miRDeep_plant.pl $predict/signatures.bst $predict/excised_precursor_struc.txt novel_tmp_dir -y > microRNA_prediction.mrd");
-#	print "\nmiRDeep_plant.pl $dir/signatures.bst $predict/excised_precursor_struc.txt tmp_dir -y > microRNA_prediction.txt\n";
-	system("rm novel_tmp_dir -rf");
-	my $tag=join "," ,@mark;
-	system("perl $scipt_path/miRNA_Express_and_sequence.pl -i microRNA_prediction.mrd -list novel_microRNA_express.txt -fa novel_microRNA_mature.fa -pre novel_microRNA_precursor.fa -tag $tag");
-}
-
-sub genome{
-	my ($file)=@_;
-	if(defined $opts{'idx'}){
-		system("perl $scipt_path/matching.pl -i $file -g $opts{gfa} -v $mis -p $t -r $hit -o $dir -index $opts{idx} -time $time") ;
-#		print "\nmatching.pl -i $file -g $opts{gfa} -v $mis -p $t -r $hit -o $dir -index $opts{idx} -time $time\n";
-	}else{
-		system("perl $scipt_path/matching.pl -i $file -g $opts{gfa} -v $mis -p $t -r $hit -o $dir -time $time") ;
-#		print "\nmatching.pl -i $file -g $opts{gfa} -v $mis -p $t -r $hit -o $dir -time $time\n";
-	}
-}
-sub rfam{
-	if (defined $opts{'idx2'}) {
-		system("perl $scipt_path/rfam.pl -i $data2 -ref $opts{rfam} -v $mis_rfam -p $t -o $dir -index $opts{idx2} -time $time");
-#		print "\nrfam.pl -i $data2 -ref $opts{rfam} -v $mis_rfam -p $t -o $dir -index $opts{idx2} -time $time\n";
-	}else{
-		system("perl $scipt_path/rfam.pl -i $data2 -ref $opts{rfam} -v $mis_rfam -p $t -o $dir -time $time");
-#		print "\nrfam.pl -i $data2 -ref $opts{rfam} -v $mis_rfam -p $t -o $dir -time $time\n";
-	}
-}
-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");
-#	print "\nfilterReadsByLength.pl -i $collapsed -o $data -min $min_nt -max $max_nt -mark $tmpmark\n";
-
-}
-sub collapse{
-	my ($ins,$data)=@_;
-		my $str="";
-	for (my $i=0;$i<@{$ins};$i++) {
-		$str .="-i $$ins[$i] ";
-	}
-	system ("perl $scipt_path/collapseReads2Tags.pl $str -mark seq -o $data -format $format");
-#	print "\ncollapseReads2Tags.pl $str -mark seq -o $data -format $format\n";
-}
-
-sub clips{
-	my ($in,$out)=@_;
-	my $adapter=$preprocess.$out."_clips_adapter.fq";
-	if($format eq "fq" || $format eq "fastq"){
-		system("fastx_clipper -a $a -M $m -Q $phred_qv -i $in -o $adapter") ;
-		print "\nfastx_clipper -a $a -M $m -Q $phred_qv -i $in -o $adapter\n";
-	}
-	if($format eq "fa" || $format eq "fasta"){
-		system("fastx_clipper -a $a -M $m -i $in -o $adapter") ;
-        #		print "\nfastx_clipper -a $a -M $m -i $in -o $adapter\n";
-	}
-	#my $clean=$preprocess.$out."_clean.fq";
-	#system("filterReadsByLength.pl -i $adapter -o $clean -min $min_nt -max $max_nt ");
-
-	return;
-}
-
-sub read_config{
-	open CON,"<$config";
-	while (my $aline=<CON>) {
-		chomp $aline;
-		my @tmp=split/\t/,$aline;
-		push @filein,$tmp[0];
-		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);
-	}
-	if ($format eq "fastq" || $format eq "fq") {
-		&checkfq($fileforcheck);
-	}
-}
-sub checkfa{
-	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;
-}
-sub checkfq{
-	my ($file_reads)=@_;
-
-	open N,"<$file_reads";
-	for (my $i=0;$i<10;$i++) {
-		my $a=<N>;
-		my $b=<N>;
-		my $c=<N>;
-		my $d=<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;
-			if($q<0){$phred_qv=33;}
-		}
-	}
-	close N;
-}
-
-sub search{
-	my ($dir,$str)=@_;
-	opendir I,$dir;
-	my @ret;
-	while (my $file=readdir I) {
-		if ($file=~/$str/) {
-			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;
-        if (length($sec) == 1) {$sec = "0"."$sec";}
-        if (length($min) == 1) {$min = "0"."$min";}
-        if (length($hour) == 1) {$hour = "0"."$hour";}
-        if (length($day) == 1) {$day = "0"."$day";}
-        if (length($month) == 1) {$month = "0"."$month";}
-        #print "$year-$month-$day $hour:$min:$sec\n";
-        return("$year-$month-$day-$hour-$min-$sec");
-}
-=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;
-        if (length($sec) == 1) {$sec = "0"."$sec";}
-        if (length($min) == 1) {$min = "0"."$min";}
-        if (length($hour) == 1) {$hour = "0"."$hour";}
-        if (length($day) == 1) {$day = "0"."$day";}
-        if (length($month) == 1) {$month = "0"."$month";}
-        #print "$year-$month-$day $hour:$min:$sec\n";
-        return("$year-$month-$day $hour:$min:$sec");
-}
-
-
-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
-                the prefix is 'h_sapiens_37_asm'.##can be null
-
--pre string,  input file #species specific microRNA precursor sequences
--mat string,  input file #species specific microRNA mature sequences
-
--rfam string,  input file# rfam database file, microRNAs must not be contained in this file## if not define, rfam small RNA will not be count.
--idx2 string,  rfam 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
-                the prefix is 'h_sapiens_37_asm'.##can be null
-
--D      If [-D] is specified,will discard rfam mapped reads(nead -rfam).
-
--a string,  ADAPTER string. default is ATCTCGTATG.
--M int,    require minimum adapter alignment length of N. If less than N nucleotides aligned with the adapter - don't clip it. 
--min int,  reads min length,default is 19.
--max int,  reads max length,default is 28.
-
--mis [int]     number of allowed mismatches when mapping reads to precursors, default 0
--e [int]     number of nucleotides upstream of the mature sequence to consider, default 2
--f [int]     number of nucleotides downstream of the mature sequence to consider, default 5
--v <int>     report end-to-end hits w/ <=v mismatches; ignore qualities,default 0; used in rfam alignment
--r int       a read is allowed to map up to this number of positions in the genome,default is 25 
-
--dis <int>   Maximal space between miRNA and miRNA* (200)
--flank <int>   Flank sequence length of miRNA precursor (10)
--mfe <folat> Maximal free energy allowed for a miRNA precursor (-20)
-
--t int,    number of threads [1]
-
--o output directory# absolute path
--h help
-USAGE
-exit(1);
-}
-