| Previous changeset 22:c0aecae5eb03 (2014-07-25) Next changeset 24:5691802f074b (2014-07-25) |
|
Commit message:
Deleted selected files |
|
removed:
miRDeep_plant.pl |
| b |
| diff -r c0aecae5eb03 -r 45de5e1ff487 miRDeep_plant.pl --- a/miRDeep_plant.pl Fri Jul 25 05:57:20 2014 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| b'@@ -1,1544 +0,0 @@\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' |