| Previous changeset 13:d1cc2e6ecf90 (2014-07-25) Next changeset 15:28e9e08d840c (2014-07-25) |
|
Commit message:
Uploaded |
|
added:
miRPlant.pl |
| b |
| diff -r d1cc2e6ecf90 -r 554fbaf5f451 miRPlant.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/miRPlant.pl Fri Jul 25 05:21:31 2014 -0400 |
| [ |
| b'@@ -0,0 +1,503 @@\n+#!/usr/bin/perl -w\n+#Filename:\n+#Author: Tian Dongmei\n+#Email: tiandm@big.ac.cn\n+#Date: 2014-4-22\n+#Modified:\n+#Description: plant microRNA prediction\n+my $version=1.00;\n+\n+use strict;\n+use Getopt::Long;\n+use threads;\n+use threads::shared;\n+use File::Path;\n+use File::Basename;\n+#use RNA;\n+use Term::ANSIColor;\n+\n+my %opts;\n+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");\n+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\n+&usage;\n+}\n+\n+my $time=&Time();\n+print "miPlant program start:\\n The time is $time!\\n";\n+print "Command line:\\n $0 @ARGV\\n";\n+\n+my $format=$opts{\'format\'};\n+if ($format ne "fastq" && $format ne "fq" && $format ne "fasta" && $format ne "fa") { \n+\t&printErr();\n+\tdie "Parameter \\"-format\\" is error! Parameter is fastq, fq, fasta or fa\\n";\n+}\n+\n+my $phred_qv=64;\n+\n+\n+my @inputfiles=@{$opts{\'i\'}};\n+my @inputtags=@{$opts{\'tag\'}};\n+\n+my $mypath=`pwd`;\n+chomp $mypath;\n+\n+my $dir=defined $opts{\'o\'} ? $opts{\'o\'} : "$mypath/miRPlant_out/";\n+\n+\n+unless ($dir=~/\\/$/) {$dir.="/";}\n+if (not -d $dir) {\n+\tmkdir $dir;\n+}\n+my $config=$dir."/input_config";\n+open CONFIG,">$config";\n+\tfor (my $i=0;$i<@inputfiles;$i++) {\n+\t\tprint CONFIG $inputfiles[$i],"\\t",$inputtags[$i],"\\n";\n+\t}\n+close CONFIG;\n+\n+my $scipt_path=defined $opts{\'path\'} ? $opts{\'path\'} : "/Users/big/galaxy-dist/tools/myTools/";\n+\n+my $a="ATCTCGTATG"; #adapter\n+if (defined $opts{\'a\'}) {$a=$opts{\'a\'};}\n+\n+my $m=6; #adapter minimum mapped nt\n+if (defined $opts{\'M\'}) {$m=$opts{\'M\'};}\n+\n+my $t=1; #threads number\n+if (defined $opts{\'t\'}) {$t=$opts{\'t\'};}\n+\n+my $min_nt=19; # minimum reads length\n+if (defined $opts{\'min\'}) {$min_nt=$opts{\'min\'};}\n+\n+my $max_nt=28; #maximum reads length\n+if (defined $opts{\'max\'}) {$max_nt=$opts{\'max\'};}\n+\n+my $mis=0; #mismatch number for microRNA\n+if (defined $opts{\'mis\'}) {$mis=$opts{\'mis\'};}\n+\n+my $mis_rfam=0;# mismatch number for rfam\n+if (defined $opts{\'v\'}) {$mis_rfam=$opts{\'v\'};}\n+\n+my $hit=25; # maximum reads mapping hits in genome\n+if (defined $opts{\'r\'}) {$hit=$opts{\'r\'};}\n+\n+my $upstream = 2; # microRNA 5\' extension\n+$upstream = $opts{\'e\'} if(defined $opts{\'e\'});\n+\n+my $downstream = 5;# microRNA 3\' extension\n+$downstream = $opts{\'f\'} if(defined $opts{\'f\'});\n+\n+my $maxd=defined $opts{\'dis\'} ? $opts{\'dis\'} : 200;\n+my $flank=defined $opts{\'flank\'} ? $opts{\'flank\'} :10;\n+my $mfe=defined $opts{\'mfe\'} ? $opts{\'mfe\'} : -20;\n+\n+$time=&Time();\n+print "$time, Checking input file!\\n";\n+\n+my (@filein,@mark,@clean);\n+#&read_config();\n+@filein=@inputfiles;\n+@mark=@inputtags;\n+\n+&checkfa($opts{pre});\n+&checkfa($opts{mat});\n+&checkfa($opts{gfa});\n+\n+\n+##### clip adpter --> clean data start\n+$time=&Time();\n+print "$time, Preprocess:\\n trim adapter, reads collapse and filter reads by length.\\n";\n+\n+$time=~s/:/-/g;\n+$time=~s/ /-/g;\n+my $preprocess=$dir."preProcess_${time}/";\n+mkdir $preprocess;\n+my $can_use_threads = eval \'use threads; 1\';\n+if ($can_use_threads) {\n+# Do processing using threads\n+\tprint "Do processing using threads\\n";\n+\tmy @filein1=@filein; my @mark1=@mark;\n+\twhile (@filein1>0) {\n+\t\tmy @thrs; my @res;\n+\t\tfor (my $i=0;$i<$t ;$i++) {\n+\t\t\tlast if(@filein1==0);\n+\t\t\tmy $in=shift @filein1;\n+\t\t\tmy $out=shift @mark1;\n+\t\t\tpush @clean,$preprocess.$out."_clips_adapter.fq";\n+\t\t\t$thrs[$i]=threads->create(\\&clips,$in,$out);\n+\t\t}\n+\t\tfor (my $i=0;$i<@thrs;$i++) {\n+\t\t\t$res[$i]=$thrs[$i]->join();\n+\t\t}\n+\t}\n+} else {\n+# Do not processing using threads\n+\tprint "Do not processing using threads\\n";\n+\tfor (my $i=0;$i<@filein ;$i++) {\n+\t\tmy $in=$filein[$i];\n+\t\tmy $out=$mark[$i];\n+\t\tpush @clean,$preprocess.$out."_clips_adapter.fq";\n+\t\t&clips($in,$out);\n+\t}\n+}\n+\n+##### clip adpter --> clean data end\n+\n+my $collapsed=$preprocess."collapse_reads.fa";\n+my $data='..b'+\tclose N;\n+}\n+\n+sub search{\n+\tmy ($dir,$str)=@_;\n+\topendir I,$dir;\n+\tmy @ret;\n+\twhile (my $file=readdir I) {\n+\t\tif ($file=~/$str/) {\n+\t\t\tpush @ret, $file;\n+\t\t}\n+\t}\n+\tclosedir I;\n+\tif (@ret != 1) {\n+\t\t&printErr();\n+\n+\t\tdie "Can not find directory or file which name has string: $str !!!\\n";\n+\t}\n+\treturn $ret[0];\n+}\n+\n+sub printErr{\n+ print STDERR color \'bold red\';\n+ print STDERR "Error: ";\n+ print STDERR color \'reset\';\n+}\n+=cut\n+sub Time{\n+ my $time=time();\n+ my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];\n+ $month++;\n+ $year+=1900;\n+ if (length($sec) == 1) {$sec = "0"."$sec";}\n+ if (length($min) == 1) {$min = "0"."$min";}\n+ if (length($hour) == 1) {$hour = "0"."$hour";}\n+ if (length($day) == 1) {$day = "0"."$day";}\n+ if (length($month) == 1) {$month = "0"."$month";}\n+ #print "$year-$month-$day $hour:$min:$sec\\n";\n+ return("$year-$month-$day-$hour-$min-$sec");\n+}\n+=cut\n+sub Time{\n+ my $time=time();\n+ my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];\n+ $month++;\n+ $year+=1900;\n+ if (length($sec) == 1) {$sec = "0"."$sec";}\n+ if (length($min) == 1) {$min = "0"."$min";}\n+ if (length($hour) == 1) {$hour = "0"."$hour";}\n+ if (length($day) == 1) {$day = "0"."$day";}\n+ if (length($month) == 1) {$month = "0"."$month";}\n+ #print "$year-$month-$day $hour:$min:$sec\\n";\n+ return("$year-$month-$day $hour:$min:$sec");\n+}\n+\n+\n+sub usage{\n+print <<"USAGE";\n+Version $version\n+Usage:\n+$0 -i -format -gfa -index -pre -mat -rfam -D -a -M -min -max -mis -e -f -v -t -o -path\n+options:\n+-i string, input file#input files information file\n+\t\t/path/filename\tmark\n+\t\t/path/filename\tmark\n+\t\t...\n+\n+-format string,#specific input rawdata file format : fastq|fq|fasta|fa\n+\n+-gfa string, input file # genome fasta. sequence file\n+-idx string, genome file index, file-prefix #(must be indexed by bowtie-build) The parameter\n+ string must be the prefix of the bowtie index. For instance, if\n+ the first indexed file is called \'h_sapiens_37_asm.1.ebwt\' then\n+ the prefix is \'h_sapiens_37_asm\'.##can be null\n+\n+-pre string, input file #species specific microRNA precursor sequences\n+-mat string, input file #species specific microRNA mature sequences\n+\n+-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.\n+-idx2 string, rfam file index, file-prefix #(must be indexed by bowtie-build) The parameter\n+ string must be the prefix of the bowtie index. For instance, if\n+ the first indexed file is called \'h_sapiens_37_asm.1.ebwt\' then\n+ the prefix is \'h_sapiens_37_asm\'.##can be null\n+\n+-D If [-D] is specified,will discard rfam mapped reads(nead -rfam).\n+\n+-a string, ADAPTER string. default is ATCTCGTATG.\n+-M int, require minimum adapter alignment length of N. If less than N nucleotides aligned with the adapter - don\'t clip it. \n+-min int, reads min length,default is 19.\n+-max int, reads max length,default is 28.\n+\n+-mis [int] number of allowed mismatches when mapping reads to precursors, default 0\n+-e [int] number of nucleotides upstream of the mature sequence to consider, default 2\n+-f [int] number of nucleotides downstream of the mature sequence to consider, default 5\n+-v <int> report end-to-end hits w/ <=v mismatches; ignore qualities,default 0; used in rfam alignment\n+-r int a read is allowed to map up to this number of positions in the genome,default is 25 \n+\n+-dis <int> Maximal space between miRNA and miRNA* (200)\n+-flank <int> Flank sequence length of miRNA precursor (10)\n+-mfe <folat> Maximal free energy allowed for a miRNA precursor (-20)\n+\n+-t int, number of threads [1]\n+\n+-o output directory# absolute path\n+-h help\n+USAGE\n+exit(1);\n+}\n+\n' |