|
18
|
1 #!/usr/bin/perl -w
|
|
|
2 ###
|
|
|
3 # But : protocol permettant la detection d'RNA non codant potentiel
|
|
|
4 #
|
|
|
5 # Entrees : fichier de mapping Smart gff3
|
|
|
6 # fichier gff des gènes
|
|
|
7 # fichier gff des clusters Cis regulateur potentiel
|
|
|
8 #
|
|
|
9 # Sortie : fichier gff des clusters ARN nc
|
|
|
10 #
|
|
|
11 ###------------------------------------------------------
|
|
|
12 use vars qw($USAGE);
|
|
|
13 use strict;
|
|
|
14
|
|
|
15 =head1 NAME
|
|
|
16
|
|
|
17 interElementGff.pl - creation of a new Gff corresponding to the region of two successive Elements
|
|
|
18
|
|
|
19 =head1 SYNOPSIS
|
|
|
20
|
|
|
21 % interElementGff.pl -i inputFile.gff3 -o outputFile.gff3 [-s 50] [-a 20] [-n seqName] [-h]
|
|
|
22
|
|
|
23 =head1 DESCRIPTION
|
|
|
24 This script will determine cluster ok ncRNA.
|
|
|
25
|
|
|
26 -i|--input fileName gff input file name
|
|
|
27 -o|--output fileName gff output file name
|
|
|
28 -n|--name seqName sequence name
|
|
|
29 -p|--print print parameters used
|
|
|
30
|
|
|
31 -f5ff n number of nt to exclude from 5' seed when gene before is Forward, seed is Forward and next gene is Forward [default 0]
|
|
|
32 -ff3f n number... " ...[default 0]
|
|
|
33
|
|
|
34 -f5fr n number... " ...[default 0]
|
|
|
35 -ff3r n number... " ...[default 0]
|
|
|
36
|
|
|
37 -fr3f n number... " ...[default 0]
|
|
|
38 -fr5f n number... " ...[default 0]
|
|
|
39
|
|
|
40 -f3rr n number... " ...[default 0]
|
|
|
41 -fr5r n number... " ...[default 0]
|
|
|
42
|
|
|
43 -r5ff n number... " ...[default 0]
|
|
|
44 -rf3f n number... " ...[default 0]
|
|
|
45
|
|
|
46 -r5fr n number... " ...[default 0]
|
|
|
47 -rf3r n number... " ...[default 0]
|
|
|
48
|
|
|
49 -r3rf n number... " ...[default 0]
|
|
|
50 -rr5f n number... " ...[default 0]
|
|
|
51
|
|
|
52 -r3rr n number... " ...[default 0]
|
|
|
53 -rr5r n number... " ...[default 0]
|
|
|
54
|
|
|
55 [-h|--help] help mode then die
|
|
|
56
|
|
|
57
|
|
|
58 USAGE_CASE
|
|
|
59
|
|
|
60 % interElementGff.pl -i inputFile.gff3 -o outputFile.gff3 -ff 53 -rr 23 -n NC_011744
|
|
|
61
|
|
|
62 BUG
|
|
|
63
|
|
|
64 Caution : input file needs to be sorted on positions
|
|
|
65
|
|
|
66 Caution : for -f/r options add +3 bp to include stop codon if not in input file
|
|
|
67
|
|
|
68 =head1 AUTHOR - CTN - apr.11
|
|
|
69 (from RNA-Vibrio/protocol_NC_V2.pl - Claire KUCHLY)
|
|
|
70
|
|
|
71 =cut
|
|
|
72 #----------------------------------------------------------------------------
|
|
|
73 # check command line :
|
|
|
74 my ($IDfile, $OutputFileName, $f5ff, $ff3f, $f5fr, $ff3r, $f3rf, $fr5f, $f3rr,$fr5r, $r5ff, $rf3f, $r5fr, $rf3r, $r3rf, $rr5f, $r3rr, $rr5r, $seqName, $printParameters) =
|
|
|
75 (undef, undef , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, "", 0) ;
|
|
|
76 if ($#ARGV==0) {
|
|
|
77 die (exec("pod2text $0\n"));
|
|
|
78 } else {
|
|
|
79 foreach my $num (0 .. $#ARGV) {
|
|
|
80 SWITCH: for ($ARGV[$num]) {
|
|
|
81 /--input|-i/ && do { $IDfile=$ARGV[$num+1];
|
|
|
82 open(F,"<$IDfile") or die "Error: Can't open \"$IDfile\", $!";
|
|
|
83 last; };
|
|
|
84 /-f5ff/ && do { $f5ff=$ARGV[$num+1]+1; last; }; # need +1 for intervall computations
|
|
|
85 /-ff3f/ && do { $ff3f=$ARGV[$num+1]+1; last; };
|
|
|
86
|
|
|
87 /-f5fr/ && do { $f5fr=$ARGV[$num+1]+1; last; };
|
|
|
88 /-ff3r/ && do { $ff3r=$ARGV[$num+1]+1; last; };
|
|
|
89
|
|
|
90 /-f3rf/ && do { $f3rf=$ARGV[$num+1]+1; last; };
|
|
|
91 /-fr5f/ && do { $fr5f=$ARGV[$num+1]+1; last; };
|
|
|
92
|
|
|
93 /-f3rr/ && do { $f3rr=$ARGV[$num+1]+1; last; };
|
|
|
94 /-fr5r/ && do { $fr5r=$ARGV[$num+1]+1; last; };
|
|
|
95
|
|
|
96 /-r5ff/ && do { $r5ff=$ARGV[$num+1]+1; last; };
|
|
|
97 /-rf3f/ && do { $rf3f=$ARGV[$num+1]+1; last; };
|
|
|
98
|
|
|
99 /-r5fr/ && do { $r5fr=$ARGV[$num+1]+1; last; };
|
|
|
100 /-rf3r/ && do { $rf3r=$ARGV[$num+1]+1; last; };
|
|
|
101
|
|
|
102 /-r3rf/ && do { $r3rf=$ARGV[$num+1]+1; last; };
|
|
|
103 /-rr5f/ && do { $rr5f=$ARGV[$num+1]+1; last; };
|
|
|
104
|
|
|
105 /-r3rr/ && do { $r3rr=$ARGV[$num+1]+1; last; };
|
|
|
106 /-rr5r/ && do { $rr5r=$ARGV[$num+1]+1; last; };
|
|
|
107
|
|
|
108 # /--name|-n/ && do { $seqName=$ARGV[$num+1]; last; };
|
|
|
109 /--print|-p/ && do { $printParameters=1; last; };
|
|
|
110 /--output|-o/ && do { $OutputFileName=$ARGV[$num+1];
|
|
|
111 open(S,">$OutputFileName") or die "Error : Can't open result file \"$OutputFileName\", $!";
|
|
|
112 last; };
|
|
|
113 /--help|-h/ && do { exec("pod2text $0\n") ; die };
|
|
|
114 }
|
|
|
115 }
|
|
|
116 if ($printParameters) {
|
|
|
117 print "
|
|
|
118 --> f5ff ",$f5ff-1," --> ff3f ",$ff3f-1," --> ;
|
|
|
119 --> f5fr ",$f5fr-1," --> ff3r ",$ff3r-1," <-- ;
|
|
|
120 --> f3rf ",$f3rf-1," <-- fr5f ",$fr5f-1," --> ;
|
|
|
121 --> f3rr ",$f3rr-1," <-- fr5r ",$fr5r-1," <-- ;
|
|
|
122 <-- r5ff ",$r5ff-1," --> rf3f ",$rf3f-1," --> ;
|
|
|
123 <-- r5fr ",$r5fr-1," --> rf3r ",$rf3r-1," <-- ;
|
|
|
124 <-- r3rf ",$r3rf-1," <-- rr5f ",$rr5f-1," --> ;
|
|
|
125 <-- r3rr ",$r3rr-1," <-- rr5r ",$rr5r-1," <-- ;\n";
|
|
|
126 }
|
|
|
127 ##NC_011753.2 RefSeq gene 367 834 . - . locus_tag=VS_0001;db_xref=GeneID:7162789
|
|
|
128 my $finSeedSens;
|
|
|
129 my $finSeedAntisens;
|
|
|
130 my $debSeedSens;
|
|
|
131 my $debSeedAntisens;
|
|
|
132 my $info_gene="";
|
|
|
133 my $sensGeneAvant = "+" ; # 1rst seed definition : geneAvant (gene[i-1]) doesn't exist
|
|
|
134 my @chromList;
|
|
|
135 while(my $ligne = <F>){
|
|
|
136 chomp($ligne);
|
|
|
137 my @list = split(/\t/,$ligne);
|
|
|
138 if ((scalar(@chromList) == 0) or ($chromList[$#chromList] ne $list[0])){
|
|
|
139 push(@chromList, $list[0]);
|
|
|
140 my $finSeedSens;
|
|
|
141 my $finSeedAntisens;
|
|
|
142 my $debSeedSens;
|
|
|
143 my $debSeedAntisens;
|
|
|
144 my $info_gene="";
|
|
|
145 my $sensGeneAvant = "+" ; # 1rst seed definition : geneAvant (gene[i-1]) doesn't exist
|
|
|
146 }
|
|
|
147 if (($sensGeneAvant eq "+") and ($list[6] eq "+")) { #CTN ie geneavant == f, geneapres == f
|
|
|
148 $debSeedSens += $f5ff;
|
|
|
149 $finSeedSens = $list[3]- $ff3f;
|
|
|
150 $debSeedAntisens += $f3rf;
|
|
|
151 $finSeedAntisens = $list[3]- $fr5f;
|
|
|
152 } elsif (($sensGeneAvant eq "+") and ($list[6] eq "-")) { #CTN ie geneaavant == f, geneapres == r
|
|
|
153 $debSeedSens += $f5fr;
|
|
|
154 $finSeedSens = $list[3]- $ff3r;
|
|
|
155 $debSeedAntisens += $f3rr;
|
|
|
156 $finSeedAntisens = $list[3]- $fr5r;
|
|
|
157 } elsif (($sensGeneAvant eq "-") and ($list[6] eq "+")) { #CTN ie geneaavant == r, geneapres == f
|
|
|
158 $debSeedSens += $r5ff;
|
|
|
159 $finSeedSens = $list[3]- $rf3f;
|
|
|
160 $debSeedAntisens += $r3rf;
|
|
|
161 $finSeedAntisens = $list[3]- $rr5f;
|
|
|
162 } else { #CTN ie geneaavant == r, geneapres == r
|
|
|
163 $debSeedSens += $r5fr;
|
|
|
164 $finSeedSens = $list[3]- $rf3r;
|
|
|
165 $debSeedAntisens += $r3rr;
|
|
|
166 $finSeedAntisens = $list[3]- $rr5r;
|
|
|
167 }
|
|
|
168 if ($debSeedSens <= 0) { $debSeedSens=1 ; } # 1srt
|
|
|
169 if ($debSeedAntisens <= 0) { $debSeedAntisens=1 ; }
|
|
|
170 if($debSeedSens < $finSeedSens){ # only "real" seed
|
|
|
171 #print "$gene_avant\nNC_011753\tperso\tseed\t$deb_seed\t$fin_seed\t.\t+\t.\tgeneavant=$info_gene;geneapres=$list[@list-1]\n$ligne\n\n";
|
|
|
172 #
|
|
|
173
|
|
|
174 print S "$list[0]\tperso\tseedIR\t$debSeedSens\t$finSeedSens\t.\t+\t.\tgeneavant=$info_gene;geneapres=$list[@list-1]\n";
|
|
|
175 }
|
|
|
176 if ($debSeedAntisens < $finSeedAntisens){
|
|
|
177 print S "$list[0]\tperso\tseedIR\t$debSeedAntisens\t$finSeedAntisens\t.\t-\t.\tgeneavant=$info_gene;geneapres=$list[@list-1]\n";
|
|
|
178 }
|
|
|
179 $sensGeneAvant = $list[6] ; # GFF : column 6 gives strand
|
|
|
180 $debSeedSens = $list[4];
|
|
|
181 $debSeedAntisens = $list[4];
|
|
|
182 $info_gene = $list[@list-1];
|
|
|
183 }
|
|
|
184 close F;
|
|
|
185 close S;
|
|
|
186 exit(0);
|
|
|
187 }
|
