| Previous changeset 34:915b62d31f9d (2014-07-31) Next changeset 36:1a6d88b42514 (2014-07-31) |
|
Commit message:
Uploaded |
|
added:
precursors.pl |
| b |
| diff -r 915b62d31f9d -r 0979f7a58686 precursors.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/precursors.pl Thu Jul 31 03:08:13 2014 -0400 |
| [ |
| b'@@ -0,0 +1,789 @@\n+#!/usr/bin/perl -w\n+#Filename:\n+#Author: Tian Dongmei\n+#Email: tiandm@big.ac.cn\n+#Date: 2013/7/19\n+#Modified:\n+#Description: \n+my $version=1.00;\n+\n+use strict;\n+use Getopt::Long;\n+use RNA;\n+\n+my %opts;\n+GetOptions(\\%opts,"map=s","g=s","d:i","f:i","o=s","e:f","s=s","h");\n+if (!(defined $opts{map} and defined $opts{g} and defined $opts{o} and defined $opts{s} ) || defined $opts{h}) { #necessary arguments\n+&usage;\n+}\n+\n+my $filein=$opts{\'map\'};\n+my $faout=$opts{\'o\'};\n+my $strout=$opts{\'s\'};\n+my $genome= $opts{\'g\'};\n+\n+my $maxd=defined $opts{\'d\'} ? $opts{\'d\'} : 200;\n+my $flank=defined $opts{\'f\'}? $opts{\'f\'} : 10;\n+\n+my $MAX_ENERGY=-18;\n+if (defined $opts{\'e\'}) {$MAX_ENERGY=$opts{\'e\'};}\n+my $MAX_UNPAIR=5;\n+my $MIN_PAIR=15;\n+my $MAX_SIZEDIFF=4;\n+my $MAX_BULGE=2;\n+my $ASYMMETRY=5;\n+my $MIN_UNPAIR=0;\n+my $MIN_SPACE=5;\n+my $MAX_SPACE=$maxd;\n+my $FLANK=$flank;\n+\n+######### load in genome sequences start ########\n+my %genome;\n+my %lng;\n+my $name;\n+open IN,"<$genome";\n+while (my $aline=<IN>) {\n+\tchomp $aline;\n+\tnext if($aline=~/^\\#/);\n+\tif ($aline=~/^>(\\S+)/) {\n+\t\t$name=$1;\n+\t\tnext;\n+\t}\n+\t$genome{$name} .=$aline;\n+}\n+close IN;\n+foreach my $key (keys %genome) {\n+\t$lng{$key}=length($genome{$key});\n+}\n+####### load in genome sequences end ##########\n+\n+my %breaks; ### reads number bigger than 3\n+open IN,"<$filein"; #input file \n+while (my $aline=<IN>) {\n+\tchomp $aline;\n+\tmy @tmp=split/\\t/,$aline;\n+\t$tmp[0]=~/_x(\\d+)$/;\n+\tmy $no=$1;\n+\tnext if($no<3);\n+\t#my $trand=&find_strand($tmp[9]);\n+\t#my @pos=split/\\.\\./,$tmp[5];\n+\tmy $end=$tmp[3]+length($tmp[4])-1;\n+\tif($tmp[1] eq "-"){$tmp[4]=revcom($tmp[4]);}\n+\tpush @{$breaks{$tmp[2]}{$tmp[1]}},[$tmp[3],$end,$no,$tmp[4]]; ### 0 base\n+}\n+close IN;\n+\n+my %cites; ### peaks\n+foreach my $chr (keys %breaks) {\n+\tforeach my $strand (keys %{$breaks{$chr}}) {\n+\t\tmy @array=@{$breaks{$chr}{$strand}};\n+\t\t@array=sort{$a->[0]<=>$b->[0]} @array;\n+\t\tfor (my $i=0;$i<@array;$i++) {\n+\t\t\tmy $start=$array[$i][0];my $end=$array[$i][1];\n+\t\t\tmy @subarray=();\n+\t\t\tpush @subarray,$array[$i];\n+\n+\t\t\tfor (my $j=$i+1;$j<@array;$j++) {\n+\t\t\t\tif ($start<$array[$j][1] && $end>$array[$j][0]) {\n+\t\t\t\t\tpush @subarray,$array[$j];\n+\t\t\t\t\t($start,$end)=&newpos($start,$end,$array[$j][0],$array[$j][1]);\n+\t\t\t\t}\n+\t\t\t\telse{\n+\t\t\t\t\t$i=$j;\n+\t\t\t\t\t&find_cites(\\@subarray,$chr,$strand);\n+\t\t\t\t\tlast;\n+\t\t\t\t}\n+\t\t\t}\n+\t\t}\n+\t}\n+}\n+\n+open FA,">$faout"; #output file \n+open STR,">$strout";\n+foreach my $chr (keys %cites) {\n+\tforeach my $strand (keys %{$cites{$chr}}) {\n+\n+\t\tmy @array2=@{$cites{$chr}{$strand}};\n+\t\t@array2=sort{$a->[0]<=>$b->[0]} @array2;\n+\t\t&excise(\\@array2,$chr,$strand);\n+\t}\n+}\n+close FA;\n+close STR;\n+sub oneCiteDn{\n+\tmy ($array,$a,$chr,$strand)=@_;\n+\n+\tmy $ss=$$array[$a][0]-$flank;\n+\t$ss=0 if($ss<0);\n+\tmy $ee=$$array[$a][1]+$maxd+$flank;\n+\t$ee=$lng{$chr} if($ee>$lng{$chr});\n+\t\n+\tmy $seq=substr($genome{$chr},$ss,$ee-$ss+1);\n+\tif($strand eq "-"){$seq=revcom($seq);}\n+\n+\tmy $val=&ffw1($seq,$$array[$a][3],$chr,$strand,$ss,$ee);\n+\treturn $val;\n+}\n+sub oneCiteUp{\n+\tmy ($array,$a,$chr,$strand)=@_;\n+\n+\tmy $ss=$$array[$a][0]-$maxd-$flank;\n+\t$ss=0 if($ss<0);\n+\tmy $ee=$$array[$a][1]+$flank;\n+\t$ee=$lng{$chr} if($ee>$lng{$chr});\n+\t\n+\tmy $seq=substr($genome{$chr},$ss,$ee-$ss+1);\n+\tif($strand eq "-"){$seq=revcom($seq);}\n+\n+\tmy $val=&ffw1($seq,$$array[$a][3],$chr,$strand,$ss,$ee);\n+\treturn $val;\n+\n+}\n+\n+sub twoCites{\n+\tmy ($array,$a,$b,$chr,$strand)=@_;\n+\n+\tmy $ss=$$array[$a][0]-$flank;\n+\t$ss=0 if($ss<0);\n+\tmy $ee=$$array[$b][1]+$flank;\n+\t$ee=$lng{$chr} if($ee>$lng{$chr});\n+\t\n+\tmy $seq=substr($genome{$chr},$ss,$ee-$ss+1);\n+\tif($strand eq "-"){$seq=revcom($seq);}\n+\n+#\tmy( $str,$mfe)=RNA::fold($seq);\n+#\treturn 0 if($mfe>$MAX_ENERGY); ### minimum mfe\n+\tmy $val=&ffw2($seq,$$array[$a][3],$$array[$b][3],$chr,$strand,$ss,$ee);\n+\t\n+\treturn $val;\n+\n+}\n+sub excise{\n+\tmy ($cluster,$chr,$strand)=@_;\n+\n+\tmy $last_pos=0;\n+\tfor (my $i=0;$i<@{$cluster};$i++) {\n+\t\tnext if($$cluster[$i][0]<$last_pos);\n+\t\tmy $ok=0;\n+\t\tfor (my $j=$i+1;$j<@{$cluster} ;$j++) {\n+\t\t\tif($'..b' s2\n+\n+\tmy $pair_num=0;\n+\tmy $overhang1=0;\n+\tmy $overhang2=0;\n+\tmy ($s1,$e1,$s2,$e2);\n+\tforeach my $i ($beg1..$end1) {\n+\t\tif (defined $table->{$i}) {\n+\t\t\tmy $j=$table->{$i};\n+\t\t\tif ($j <= $end2 && $j >= $beg2) {\n+\t\t\t\t$s1=$i;\n+\t\t\t\t$e2=$j;\n+\t\t\t\tlast;\n+\t\t\t}\n+\t\t}\n+\t}\n+\tforeach my $i (reverse ($beg1..$end1)) {\n+\t\tif (defined $table->{$i}) {\n+\t\t\tmy $j=$table->{$i};\n+\t\t\tif ($j <= $end2 && $j >= $beg2) {\n+\t\t\t\t$e1=$i;\n+\t\t\t\t$s2=$j;\n+\t\t\t\tlast;\n+\t\t\t}\n+\t\t}\n+\t}\n+\n+#\tprint "$beg1,$end1 $s1,$e1\\n";\n+#\tprint "$beg2,$end2 $s2,$e2\\n";\n+\n+\tforeach my $i ($beg1..$end1) {\n+\t\tif (defined $table->{$i}) {\n+\t\t\tmy $j=$table->{$i};\n+\t\t\tif ($j <= $end2 && $j >= $beg2) {\n+\t\t\t\t++$pair_num;\n+\t\t\t}\n+\t\t}\n+\t}\n+\tif (defined $s1 && defined $e2) {\n+\t\t$overhang1=($end2-$e2)-($s1-$beg1);\n+\t}\n+\tif (defined $e1 && defined $s2) {\n+\t\t$overhang2=($end1-$e1)-($s2-$beg2);\n+\t}\n+\t\n+\tif ($pair_num < 13) {\n+\t\t$like_mir_duplex=0;\n+\t}\n+\tif ($overhang1 < 0 && $overhang2 < 0) {\n+\t\t$like_mir_duplex=0;\n+\t}\n+\treturn ($like_mir_duplex,$pair_num,$overhang1,$overhang2);\n+}\n+sub parse_struct {\n+\tmy $struct=shift;\n+\tmy $table=shift;\n+\n+\tmy @t=split(\'\',$struct);\n+\tmy @lbs; # left brackets\n+\tforeach my $k (0..$#t) {\n+\t\tif ($t[$k] eq "(") {\n+\t\t\tpush @lbs, $k+1;\n+\t\t}\n+\t\telsif ($t[$k] eq ")") {\n+\t\t\tmy $lb=pop @lbs;\n+\t\t\tmy $rb=$k+1;\n+\t\t\t$table->{$lb}=$rb;\n+\t\t\t$table->{$rb}=$lb;\n+\t\t}\n+\t}\n+\tif (@lbs) {\n+\t\twarn "unbalanced RNA struct.\\n";\n+\t}\n+}\n+sub which_arm {\n+\tmy $substruct=shift;\n+\tmy $arm;\n+\tif ($substruct=~/\\(/ && $substruct=~/\\)/) {\n+\t\t$arm="-";\n+\t}\n+\telsif ($substruct=~/\\(/) {\n+\t\t$arm="5p";\n+\t}\n+\telse {\n+\t\t$arm="3p";\n+\t}\n+\treturn $arm;\n+}\n+sub biggest_bulge {\n+\tmy $struct=shift;\n+\tmy $bulge_size=0;\n+\tmy $max_bulge=0;\n+\twhile ($struct=~/(\\.+)/g) {\n+\t\t$bulge_size=length $1;\n+\t\tif ($bulge_size > $max_bulge) {\n+\t\t\t$max_bulge=$bulge_size;\n+\t\t}\n+\t}\n+\treturn $max_bulge;\n+}\n+sub get_asy {\n+\tmy($table,$a1,$a2)=@_;\n+\tmy ($pre_i,$pre_j);\n+\tmy $asymmetry=0;\n+\tforeach my $i ($a1..$a2) {\n+\t\tif (defined $table->{$i}) {\n+\t\t\tmy $j=$table->{$i};\n+\t\t\tif (defined $pre_i && defined $pre_j) {\n+\t\t\t\tmy $diff=($i-$pre_i)+($j-$pre_j);\n+\t\t\t\t$asymmetry += abs($diff);\n+\t\t\t}\n+\t\t\t$pre_i=$i;\n+\t\t\t$pre_j=$j;\n+\t\t}\n+\t}\n+\treturn $asymmetry;\n+}\n+\n+sub peaks{\n+\tmy @cluster=@{$_[0]};\n+\n+\treturn if(@cluster<1);\n+\n+\tmy $max=0; my $index=-1;\n+\tfor (my $i=0;$i<@cluster;$i++) {\n+\t\tif($cluster[$i][2]>$max){\n+\t\t\t$max=$cluster[$i][2];\n+\t\t\t$index=$i;\n+\t\t}\n+\t}\n+#\t&excise(\\@cluster,$index,$_[1],$_[2]);\n+\treturn($index);\n+}\n+\n+sub find_cites{\n+\tmy @tmp=@{$_[0]};\n+\tmy $i=&peaks(\\@tmp);\n+\t\n+\tmy $start=$tmp[$i][0];\n+\tmy $total=0; my $node5=0;\n+\tfor (my $j=0;$j<@tmp ;$j++) {\n+\t\t$total+=$tmp[$j][2];\n+\t\t$node5 +=$tmp[$j][2] if($tmp[$j][0]-$start<=2 && $tmp[$j][0]-$start>=-2);\n+\t}\n+\tpush @{$cites{$_[1]}{$_[2]}},$tmp[$i] if($node5/$total>0.80 && $tmp[$i][2]/$node5>0.5);\n+}\n+\n+sub newpos{\n+\tmy ($a,$b,$c,$d)=@_;\n+\tmy $s= $a>$c ? $c : $a;\n+\tmy $e=$b>$d ? $b : $d;\n+\treturn($s,$e);\n+}\n+\n+sub rev{\n+\n+ my($sequence)=@_;\n+\n+ my $rev=reverse $sequence; \n+\n+ return $rev;\n+}\n+\n+sub com{\n+\n+ my($sequence)=@_;\n+\n+ $sequence=~tr/acgtuACGTU/TGCAATGCAA/; \n+ \n+ return $sequence;\n+}\n+\n+sub revcom{\n+\n+ my($sequence)=@_;\n+\n+ my $revcom=rev(com($sequence));\n+\n+ return $revcom;\n+}\n+\n+sub find_strand{\n+\n+ #A subroutine to find the strand, parsing different blast formats\n+ my($other)=@_;\n+\n+ my $strand="+";\n+\n+ if($other=~/-/){\n+\t$strand="-";\n+ }\n+\n+ if($other=~/minus/i){\n+\t$strand="-";\n+ }\n+\n+ return($strand);\n+}\n+sub usage{\n+print <<"USAGE";\n+Version $version\n+Usage:\n+$0 -map -g -d -f -o -s -e\n+options:\n+ -map input file# align result # bst. format\n+ -g input file # genome sequence fasta format\n+ -d <int> Maximal space between miRNA and miRNA* (200)\n+ -f <int> Flank sequence length of miRNA precursor (10)\n+ -o output file# percursor fasta file\n+ -s output file# precursor structure file\n+ -e <folat> Maximal free energy allowed for a miRNA precursor (-18 kcal/mol)\n+\n+ -h help\n+USAGE\n+exit(1);\n+}\n+\n' |