| Previous changeset 11:861eb4fbfbd0 (2014-07-25) Next changeset 13:d1cc2e6ecf90 (2014-07-25) |
|
Commit message:
Uploaded |
|
added:
miRDeep_plant.pl |
| b |
| diff -r 861eb4fbfbd0 -r dc5a29826c7d miRDeep_plant.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/miRDeep_plant.pl Fri Jul 25 05:20:56 2014 -0400 |
| b |
| b'@@ -0,0 +1,1544 @@\n+#!/usr/bin/perl\n+\n+use warnings;\n+use strict;\n+use Getopt::Std;\n+\n+\n+\n+################################# MIRDEEP #################################################\n+\n+################################## USAGE ##################################################\n+\n+\n+my $usage=\n+"$0 file_signature file_structure temp_out_directory\n+\n+This is the core algorithm of miRDeep. It takes as input a file in blastparsed format with\n+information on the positions of reads aligned to potential precursor sequences (signature).\n+It also takes as input an RNAfold output file, giving information on the sequence, structure\n+and mimimum free energy of the potential precursor sequences.\n+\n+Extra arguments can be given. -s specifies a fastafile containing the known mature miRNA\n+sequences that should be considered for conservation purposes. -t prints out the potential\n+precursor sequences that do _not_ exceed the cut-off (default prints out the sequences that\n+exceeds the cut-off). -u gives limited output, that is only the ids of the potential precursors\n+that exceed the cut-off. -v varies the cut-off. -x is a sensitive option for Sanger sequences\n+obtained through conventional cloning. -z consider the number of base pairings in the lower\n+stems (this option is not well tested).\n+\n+-h print this usage\n+-s fasta file with known miRNAs\n+#-o temp directory ,maked befor running the program.\n+-t print filtered\n+-u limited output (only ids)\n+-v cut-off (default 1)\n+-x sensitive option for Sanger sequences\n+-y use Randfold\n+-z consider Drosha processing\n+";\n+\n+\n+\n+\n+\n+############################################################################################\n+\n+################################### INPUT ##################################################\n+\n+\n+#signature file in blast_parsed format\n+my $file_blast_parsed=shift or die $usage;\n+\n+#structure file outputted from RNAfold\n+my $file_struct=shift or die $usage;\n+\n+my $tmpdir=shift or die $usage;\n+#options\n+my %options=();\n+getopts("hs:tuv:xyz",\\%options);\n+\n+\n+\n+\n+\n+\n+#############################################################################################\n+\n+############################# GLOBAL VARIABLES ##############################################\n+\n+\n+#parameters\n+my $nucleus_lng=11;\n+\n+my $score_star=3.9;\n+my $score_star_not=-1.3;\n+my $score_nucleus=7.63;\n+my $score_nucleus_not=-1.17;\n+my $score_randfold=1.37;\n+my $score_randfold_not=-3.624;\n+my $score_intercept=0.3;\n+my @scores_stem=(-3.1,-2.3,-2.2,-1.6,-1.5,0.1,0.6,0.8,0.9,0.9,0);\n+my $score_min=1;\n+if($options{v}){$score_min=$options{v};}\n+if($options{x}){$score_min=-5;}\n+\n+my $e=2.718281828;\n+\n+#hashes\n+my %hash_desc;\n+my %hash_seq;\n+my %hash_struct;\n+my %hash_mfe;\n+my %hash_nuclei;\n+my %hash_mirs;\n+my %hash_query;\n+my %hash_comp;\n+my %hash_bp;\n+\n+#other variables\n+my $subject_old;\n+my $message_filter;\n+my $message_score;\n+my $lines;\n+my $out_of_bound;\n+\n+\n+\n+##############################################################################################\n+\n+################################ MAIN ###################################################### \n+\n+\n+#print help if that option is used\n+if($options{h}){die $usage;}\n+unless ($tmpdir=~/\\/$/) {$tmpdir .="/";}\n+if(!(-s $tmpdir)){mkdir $tmpdir;}\n+$tmpdir .="TMP_DIR/";\n+mkdir $tmpdir;\n+\n+#parse structure file outputted from RNAfold\n+parse_file_struct($file_struct);\n+\n+#if conservation is scored, the fasta file of known miRNA sequences is parsed\n+if($options{s}){create_hash_nuclei($options{s})};\n+\n+#parse signature file in blast_parsed format and resolve each potential precursor\n+parse_file_blast_parsed($file_blast_parsed);\n+`rm -rf $tmpdir`;\n+exit;\n+\n+\n+\n+\n+##############################################################################################\n+\n+############################## SUBROUTINES ###################################################\n+\n+\n+\n+sub parse_file_blast_parsed{\n+\n+# read through the signature blastparsed file, fills up a hash with information on '..b'with\n+# the n\'t nucleotide on the plus strand\n+\n+# This subroutine reverses the begin and end positions of positions of the minus strand so that they\n+# are relative to the 5\' end of the plus strand\t\n+ \n+ my($beg,$end,$lng)=@_;\n+ \n+ my $new_end=$lng-$beg+1;\n+ my $new_beg=$lng-$end+1;\n+ \n+ return($new_beg,$new_end);\n+}\n+\n+sub round {\n+\n+ #rounds to nearest integer\n+ \n+ my($number) = shift;\n+ return int($number + .5);\n+ \n+}\n+\n+\n+sub rev{\n+\n+ #reverses the order of nucleotides in a sequence\n+\n+ my($sequence)=@_;\n+\n+ my $rev=reverse $sequence; \n+\n+ return $rev;\n+}\n+\n+sub com{\n+\n+ #the complementary of a sequence\n+\n+ my($sequence)=@_;\n+\n+ $sequence=~tr/acgtuACGTU/TGCAATGCAA/; \n+ \n+ return $sequence;\n+}\n+\n+sub revcom{\n+ \n+ #reverse complement\n+\n+ my($sequence)=@_;\n+\n+ my $revcom=rev(com($sequence));\n+\n+ return $revcom;\n+}\n+\n+\n+sub max2 {\n+\n+ #max of two numbers\n+ \n+ my($a, $b) = @_;\n+ return ($a>$b ? $a : $b);\n+}\n+\n+sub min2 {\n+\n+ #min of two numbers\n+\n+ my($a, $b) = @_;\n+ return ($a<$b ? $a : $b);\n+}\n+\n+\n+\n+sub score_freq{\n+\n+# scores the count of reads that map to the potential precursor\n+# Assumes geometric distribution as described in methods section of manuscript\n+\n+ my $freq=shift;\n+\n+ #parameters of known precursors and background hairpins\n+ my $parameter_test=0.999;\n+ my $parameter_control=0.6;\n+\n+ #log_odds calculated directly to avoid underflow\n+ my $intercept=log((1-$parameter_test)/(1-$parameter_control));\n+ my $slope=log($parameter_test/$parameter_control);\n+ my $log_odds=$slope*$freq+$intercept;\n+\n+ #if no strong evidence for 3\' overhangs, limit the score contribution to 0\n+ unless($options{x} or $hash_comp{"star_read"}){$log_odds=min2($log_odds,0);}\n+ \n+ return $log_odds;\n+}\n+\n+\n+\n+##sub score_mfe{\n+\n+# scores the minimum free energy in kCal/mol of the potential precursor\n+# Assumes Gumbel distribution as described in methods section of manuscript\n+\n+## my $mfe=shift;\n+\n+ #numerical value, minimum 1\n+## my $mfe_adj=max2(1,-$mfe);\n+\n+ #parameters of known precursors and background hairpins, scale and location\n+## my $prob_test=prob_gumbel_discretized($mfe_adj,5.5,32);\n+## my $prob_background=prob_gumbel_discretized($mfe_adj,4.8,23);\n+\n+## my $odds=$prob_test/$prob_background;\n+## my $log_odds=log($odds);\n+\n+## return $log_odds;\n+##}\n+\n+sub score_mfe{\n+#\tuse bignum;\n+\n+# scores the minimum free energy in kCal/mol of the potential precursor\n+# Assumes Gumbel distribution as described in methods section of manuscript\n+\n+ my ($mfe,$mlng)=@_;\n+\n+ #numerical value, minimum 1\n+ my $mfe_adj=max2(1,-$mfe);\n+my $mfe_adj1=$mfe/$mlng;\n+ #parameters of known precursors and background hairpins, scale and location\n+\tmy $a=1.339e-12;my $b=2.778e-13;my $c=45.834;\n+\tmy $ev=$e**($mfe_adj1*$c);\n+\tprint STDERR "\\n***",$ev,"**\\t",$ev+$b,"\\t";\n+\tmy $log_odds=($a/($b+$ev));\n+\n+\n+\tmy $prob_test=prob_gumbel_discretized($mfe_adj,5.5,32);\n+ my $prob_background=prob_gumbel_discretized($mfe_adj,4.8,23);\n+\n+ my $odds=$prob_test/$prob_background;\n+ my $log_odds_2=log($odds);\n+\tprint STDERR "log_odds :",$log_odds,"\\t",$log_odds_2,"\\n";\n+ return $log_odds;\n+}\n+\n+\n+\n+sub prob_gumbel_discretized{\n+\n+# discretized Gumbel distribution, probabilities within windows of 1 kCal/mol\n+# uses the subroutine that calculates the cdf to find the probabilities\n+\n+ my ($var,$scale,$location)=@_;\n+\n+ my $bound_lower=$var-0.5;\n+ my $bound_upper=$var+0.5;\n+\n+ my $cdf_lower=cdf_gumbel($bound_lower,$scale,$location);\n+ my $cdf_upper=cdf_gumbel($bound_upper,$scale,$location);\n+\n+ my $prob=$cdf_upper-$cdf_lower;\n+\n+ return $prob;\n+}\n+\n+\n+sub cdf_gumbel{\n+\n+# calculates the cumulative distribution function of the Gumbel distribution\n+\n+ my ($var,$scale,$location)=@_;\n+\n+ my $cdf=$e**(-($e**(-($var-$location)/$scale)));\n+\n+ return $cdf;\n+}\n+\n' |