Mercurial > repos > yufei-luo > s_mart

comparison SMART/bacteriaRegulatoryRegion_Detection/splitTranscriptGff.pl @ 31:0ab839023fe4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author m-zytnicki
date Tue, 30 Apr 2013 14:33:21 -0400
parents 94ab73e8a190
children
comparison
equal deleted inserted replaced
30:5677346472b5 31:0ab839023fe4
1 #!/usr/bin/perl -w
2 ###
3 # Main : defining utr and intergenic operonic intervalles from a transcripts file following a referencies file
4 #
5 # Input : 2 gff files to intersect, transcript queries vs referencies
6 #
7 # Output : resulting gff file printing to standard output
8 #
9 ###------------------------------------------------------
10 use vars qw($USAGE);
11 use strict;
12
13 =head1 NAME
14
15 splitTranscriptGff.pl - compare 2 input gff files and define utr and intergenic operonic intervalles by couple of overlapping elements
16
17 =head1 SYNOPSIS
18
19 % intervallsExtractorGff.pl -i referencies.gff -j transcriptQueries.gff -s strand [-h]
20
21 =head1 DESCRIPTION
22 This script will intersect 2 gff files and compute distance between 2 successives lines. Take care both of sorting by positions the input files and of that referencies are included in transcriptQueries.
23
24 -i|--input1 fileName gff input file name: included elements
25 -j|--input2 fileName gff input file name: extended elements
26 [-s|--strand] [s|d] s for single strand (colinear) or d for double strands (antisense) [default d]
27 [-h|--help] help mode then die
28
29 =head1 USECASE
30 Define many fragments for each extended element (transcript): UTR5, gene, UTR3, "inOperon" for intergenomic region between 2 genes
31 intervallsExtractorGff.pl -i CDSannotations.gff -j RNAseqTranscripts.gff > UTRsGenesOperonsLists.gff;
32
33 =head1 KWON BUGS
34 No disjonction of overlapping elements of the included elements (-i file).
35 In usecase, overlapping genes are fused in one long gene.
36
37 =head1 AUTHOR
38 Claire Toffano-Nioche - sep.11
39
40 =cut
41 #-----------------------
42 sub feedPositionTab { my ($val, $pF, $pB, @info) = @_ ;
43 #print "feedPositionTab::$#info, ", ($#info+1)/4," \n";
44 for (my $i=0 ; $i <= $#info ; $i+=4) { # for each extended element
45 #print "....$info[$i+2]\n";
46 if ($info[$i+3] =~ /\+/) {
47 for (my $c = $info[$i+1] ; $c <= $info[$i+2] ; $c++) { @$pF[$c]=$val } ; # sequence Forward
48 } else {
49 for (my $c = $info[$i+1] ; $c <= $info[$i+2] ; $c++) { @$pB[$c]=$val } ; # sequence Backward
50 }
51 }
52 #print "feedPos...:: ", join(".", @$pF[0..100]), "\n";
53 #print "feedPos...:: ", join(".", @$pB[0..100]), "\n";
54 }
55 #-----------------------
56 sub recupInfo { my ($pInfo, @lines) = @_ ;
57 for (my $i=0 ; $i <= ($#lines+1)*4-1 ; $i+=4) {
58 my @line = split("\t",$lines[$i/4]);
59 push(@$pInfo, $line[0], $line[3], $line[4], $line[6]) ; # 0=nom, 3=debut, 4=fin, 6=sens
60 }
61 #print "recupInfo::fin=", ($#lines+1)*4, "\n" ;
62 }
63 #-----------------------
64 sub tagName { my ($seqN, $posB, $posE, $strand) = @_ ;
65 my $tagN=$seqN.$strand.$posB."..".$posE;
66 #print "tagName:",join("_",@_)," et tagName:$tagN\n";
67 return $tagN;
68 }
69 #-----------------------
70 sub transitionAnalysis {
71 my ($pos, $seq, $s, $pdebAmont, $pfinAmont, $pdebIn, $pfinIn, $pdebAval, $pfinAval, $ptag) = @_ ;
72 my $enCours = 0 ; my $precedant = 0 ;
73 $enCours = @$ptag[$pos] ;
74 $precedant = ($s =~ /\+/?@$ptag[$pos-1]:@$ptag[$pos+1]) ;
75 if ($enCours ne $precedant) {
76 #print "transi...:: $s, $pos, $precedant, $enCours\n";
77 #print "transition::$$pdebAmont, $$pfinAmont, $$pdebIn, $$pfinIn, $$pdebAval, $$pfinAval\n";
78 SWITCH: for ($precedant.$enCours) {
79 /01/ && do { $$pdebAmont = $pos ; last SWITCH ;};
80 /02/ && do { $$pdebIn = $pos ; last SWITCH ;};
81 /10/ && do { $$pfinAval = ($s =~/\+/?$pos-1:$pos+1) ;
82 if (($s =~ /\+/)and ($$pdebAval!=$$pfinAval)) {
83 printf "%s\tsplit\tutr3\t%s\t%s\t.\t%s\t.\tName=%s;\n",
84 $seq, $$pdebAval, $$pfinAval, $s, &tagName($seq, $$pdebAval, $$pfinAval, $s) ;
85 #if ($$pdebAval==$$pfinAval) { print "transition 10 +\n"};
86 } elsif ($$pfinAval!=$$pdebAval) {
87 printf "%s\tsplit\tutr3\t%s\t%s\t.\t%s\t.\tName=%s;\n",
88 $seq, $$pfinAval, $$pdebAval, $s, &tagName($seq, $$pfinAval, $$pdebAval, $s) ;
89 #if ($$pfinAval==$$pdebAval){ print "transition 10 -\n"};
90 }
91 $$pdebAval = 0 ; $$pfinAval = 0 ;
92 last SWITCH ;
93 };
94 /12/ && do { $$pdebIn = $pos ; $$pfinAmont=($s =~/\+/?$pos-1:$pos+1) ;
95 my $type="utr5";
96 if ($$pdebAmont == 0) { # in case of interOperon : utr5-CDS-interOperon-CDS-utr3
97 $$pdebAmont=($s =~/\+/?$$pfinIn+1:$$pfinIn-1) ;
98 $type="inOperon" ;
99 }
100 if (($s =~ /\+/) and ($$pdebAmont!=$$pfinAmont)) {
101 printf "%s\tsplit\t%s\t%s\t%s\t.\t%s\t.\tName=%s;\n",
102 $seq, $type, $$pdebAmont, $$pfinAmont, $s, &tagName($seq, $$pdebAmont, $$pfinAmont, $s) ;
103 # if ($$pdebAmont==$$pfinAmont) { print "transition 12 +\n"};
104 } elsif ($$pfinAmont!=$$pdebAmont) {
105 printf "%s\tsplit\t%s\t%s\t%s\t.\t%s\t.\tName=%s;\n",
106 $seq, $type, $$pfinAmont, $$pdebAmont, $s, &tagName($seq, $$pfinAmont, $$pdebAmont, $s) ;
107 #if ($$pfinAmont==$$pdebAmont) { print "transition 12 -\n"} ;
108 }
109 $$pdebAmont = 0 ; $$pfinAmont = 0 ;
110 last SWITCH ;
111 };
112 /20/ && do { $$pfinIn=($s =~/\+/?$pos-1:$pos+1) ;
113 if (($s =~ /\+/) and ($$pdebIn!=$$pfinIn)) {
114 printf "%s\tsplit\tgene\t%s\t%s\t.\t%s\t.\tName=%s;\n",
115 $seq, $$pdebIn, $$pfinIn, $s, &tagName($seq, $$pdebIn, $$pfinIn, $s) ;
116 } elsif ($$pfinIn!=$$pdebIn) {
117 printf "%s\tsplit\tgene\t%s\t%s\t.\t%s\t.\tName=%s;\n",
118 $seq, $$pfinIn, $$pdebIn, $s, &tagName($seq, $$pfinIn, $$pdebIn, $s) ;
119 }
120 $$pdebIn = 0 ; $$pfinIn = 0 ;
121 last SWITCH ;
122 };
123 /21/ && do { $$pdebAval=$pos ; $$pfinIn=($s =~/\+/?$pos-1:$pos+1) ;
124 if (($s =~ /\+/) and ($$pdebIn!=$$pfinIn)) {
125 printf "%s\tsplit\tgene\t%s\t%s\t.\t%s\t.\tName=%s;\n",
126 $seq, $$pdebIn, $$pfinIn, $s, &tagName($seq, $$pdebIn, $$pfinIn, $s) ;
127 } elsif ($$pfinIn!=$$pdebIn) {
128 printf "%s\tsplit\tgene\t%s\t%s\t.\t%s\t.\tName=%s;\n",
129 $seq, $$pfinIn, $$pdebIn, $s, &tagName($seq, $$pfinIn, $$pdebIn, $s) ;
130 }
131 #$$pdebIn = 0 ; $$pfinIn = 0 ;
132 last SWITCH ;
133 };
134 }
135 }
136 }
137 #-----------------------
138 my ($fileNameI, $fileNameE, $strand) = ("", "", 0) ;
139 # command line check
140 foreach my $num (0 .. $#ARGV) {
141 SWITCH: for ($ARGV[$num]) {
142 /--input1|-i/ && do {
143 $fileNameI=$ARGV[$num+1];
144 open ( fichierGffI, "< $fileNameI" ) or die "Can't open gff file: \"$fileNameI\"\n" ;
145 last };
146 /--input2|-j/ && do {
147 $fileNameE=$ARGV[$num+1];
148 open ( fichierGffE, "< $fileNameE" ) or die "Can't open gff file: \"$fileNameE\"\n" ;
149 last };
150 /--strand|-s/ && do {
151 if ($ARGV[$num+1] eq "s") { $strand=1};
152 last };
153 /--help|-h/ && do { exec("pod2text $0\n") ; die };
154 }
155 }
156 # memory declarations:
157 my @infoI ; my @infoE ;
158 my $seqName ;
159 my @tagF ; my @tagB ; # Forward and Backward sequence
160 # data retrieval:
161 my @linesI = <fichierGffI> ; my @linesE = <fichierGffE> ;
162 close fichierGffI ; close fichierGffE ;
163 #print "gff files read ; number of lines : $#lines1 + $#lines2\n";
164 # positions management
165 &recupInfo(\@infoI, @linesI) ;
166 &recupInfo(\@infoE, @linesE) ;
167 # treatement:
168 # transform gff lines into chromosomal position tags : 0 for nothing, 1 resp. 2 for extended resp. included elements
169 if (($#infoI) and ($#infoE)) {
170 $seqName=$infoI[0] ;
171 #print "fin : $infoE[$#infoE-1]\n";
172 for (my $i=0 ; $i <= $infoE[$#infoE-1] ; $i++) { $tagF[$i] = 0 ; $tagB[$i] = 0 ; } ; # "O" tag in all chr. positions
173 #print "seqName : $seqName\n" ;
174 &feedPositionTab(1, \@tagF, \@tagB, @infoE) ; # "1" tag for all extended elements
175 &feedPositionTab(2, \@tagF, \@tagB, @infoI) ; # "2" tag for all included elements
176 #print join("", @tagF), "\n";
177 #print join("", @tagB), "\n";
178 # transition management:
179 my ($beginUpstream, $endUpstream, $beginIncluded, $endIncluded, $beginDownstream, $endDownstream)
180 = (0, 0, 0, 0, 0, 0) ;
181 for (my $i=1 ; $i <= $#tagF-1 ; $i+=1) {
182 &transitionAnalysis($i, $seqName, "+", \$beginUpstream, \$endUpstream, \$beginIncluded, \$endIncluded, \$beginDownstream, \$endDownstream, \@tagF) ;
183 }
184 ($beginUpstream, $endUpstream, $beginIncluded, $endIncluded, $beginDownstream, $endDownstream) = ($infoE[$#infoE-1], $infoE[$#infoE-1], $infoE[$#infoE-1], $infoE[$#infoE-1], $infoE[$#infoE-1], $infoE[$#infoE-1]) ;
185 for (my $i=$#tagB-1 ; $i >= 1 ; $i-=1) {
186 &transitionAnalysis($i, $seqName, "-", \$beginUpstream, \$endUpstream, \$beginIncluded, \$endIncluded, \$beginDownstream, \$endDownstream, \@tagB) ;
187 }
188 }
189 exit(0) ;