Mercurial > repos > clairetn > detrprok_scripts

comparison detrprok_scripts/splitTranscriptGff.pl @ 0:a53eb951b164 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author clairetn
date Mon, 25 Mar 2013 05:39:29 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a53eb951b164
1 #!/usr/bin/perl -w
2 ###
3 # But : croiser
4 #
5 # Entrees : 2 fichiers gff à croiser
6 #
7 # Sortie : gff affiche a l'ecran
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 intervalls by couple of overlapping elements
16
17 =head1 SYNOPSIS
18
19 % intervallsExtractorGff.pl -i file1.gff -j file2.gff -s strand [-h]
20
21 =head1 DESCRIPTION
22 This script will sort 2 gff files and compute distance between 2 successives lines.
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 Compare 2 input gff files: an annotations file (included elements) and transcription units file (extended elements).
31 For each couple of overlapping elements, split the transcription unit in 5'utr, CDS, 3'utr or "operon" in case of successives genes included in the transcription unit.
32
33 =head1 KWON BUG
34 Fisrt and last elements on the genome sequence can be misclassified in the result file.
35
36 =head1 AUTHOR
37 Claire Toffano-Nioche - sep.11
38
39 =cut
40 #-----------------------
41 sub feedPositionTab { my ($val, $pF, $pB, @info) = @_ ;
42 #print "feedPositionTab::$#info, ", ($#info+1)/4," \n";
43 for (my $i=0 ; $i <= $#info ; $i+=4) { # for each extended element
44 #print "....$info[$i+2]\n";
45 if ($info[$i+3] =~ /\+/) {
46 for (my $c = $info[$i+1] ; $c <= $info[$i+2] ; $c++) { @$pF[$c]=$val } ; # Forward sequence
47 } else {
48 for (my $c = $info[$i+1] ; $c <= $info[$i+2] ; $c++) { @$pB[$c]=$val } ; # Backward sequence
49 }
50 }
51 #print "feedPos...:: ", join(".", @$pF[0..100]), "\n";
52 #print "feedPos...:: ", join(".", @$pB[0..100]), "\n";
53 }
54 #-----------------------
55 sub recupInfo { my ($pInfo, @lines) = @_ ;
56 my $i=0 ;
57 while ($i<=$#lines){
58 chomp($lines[$i]);
59 my @line = split("\t",$lines[$i]);
60 if (defined ($line[0])) {
61 if (!($line[0] =~ m/^\s*$|^#/)) { # skip both null and comment lines
62 push(@$pInfo, $line[0], $line[3], $line[4], $line[6]) ; # 0=name, 3=begin, 4=end, 6=strand
63 }
64 }
65 $i=$i+1 ;
66 }
67 #print "recupInfo::end=", $i, "\n" ;
68 }
69 #-----------------------
70 sub tagName { my ($seqN, $posB, $posE, $strand) = @_ ;
71 my $tagN=$seqN.$strand.$posB."..".$posE;
72 #print "tagName::",join("_",@_)," and tagName:$tagN\n";
73 return $tagN;
74 }
75 #-----------------------
76 sub transitionAnalysis {
77 my ($pos, $seq, $s, $pdebAmont, $pfinAmont, $pdebIn, $pfinIn, $pdebAval, $pfinAval, $ptag) = @_ ;
78 my $enCours = @$ptag[$pos] ;
79 my $precedant = ($s =~ /\+/?@$ptag[$pos-1]:@$ptag[$pos+1]) ;
80 if ($enCours ne $precedant) {
81 #print "transi...:: $s, $pos, $precedant, $enCours\n";
82 #print "transition::$$pdebAmont, $$pfinAmont, $$pdebIn, $$pfinIn, $$pdebAval, $$pfinAval\n";
83 SWITCH: for ($precedant.$enCours) {
84 /01/ && do { $$pdebAmont = $pos ; last SWITCH ;};
85 /02/ && do { $$pdebIn = $pos ; last SWITCH ;};
86 /10/ && do { $$pfinAval = ($s =~/\+/?$pos-1:$pos+1) ;
87 if (($s =~ /\+/)and ($$pdebAval!=$$pfinAval)) {
88 printf "%s\tsplit\tutr3\t%s\t%s\t.\t%s\t.\tName=%s;\n",
89 $seq, $$pdebAval, $$pfinAval, $s, &tagName($seq, $$pdebAval, $$pfinAval, $s) ;
90 #if ($$pdebAval==$$pfinAval) { print "transition 10 +\n"};
91 } elsif ($$pfinAval!=$$pdebAval) {
92 printf "%s\tsplit\tutr3\t%s\t%s\t.\t%s\t.\tName=%s;\n",
93 $seq, $$pfinAval, $$pdebAval, $s, &tagName($seq, $$pfinAval, $$pdebAval, $s) ;
94 #if ($$pfinAval==$$pdebAval){ print "transition 10 -\n"};
95 }
96 $$pdebAval = 0 ; $$pfinAval = 0 ;
97 last SWITCH ;
98 };
99 /12/ && do { $$pdebIn = $pos ; $$pfinAmont=($s =~/\+/?$pos-1:$pos+1) ;
100 my $type="utr5";
101 if ($$pdebAmont == 0) { # in case of interOperon : utr5-CDS-interOperon-CDS-utr3
102 $$pdebAmont=($s =~/\+/?$$pfinIn+1:$$pfinIn-1) ;
103 $type="operon" ;
104 }
105 if (($s =~ /\+/) and ($$pdebAmont!=$$pfinAmont)) {
106 printf "%s\tsplit\t%s\t%s\t%s\t.\t%s\t.\tName=%s;\n",
107 $seq, $type, $$pdebAmont, $$pfinAmont, $s, &tagName($seq, $$pdebAmont, $$pfinAmont, $s) ;
108 # if ($$pdebAmont==$$pfinAmont) { print "transition 12 +\n"};
109 } elsif ($$pfinAmont!=$$pdebAmont) {
110 printf "%s\tsplit\t%s\t%s\t%s\t.\t%s\t.\tName=%s;\n",
111 $seq, $type, $$pfinAmont, $$pdebAmont, $s, &tagName($seq, $$pfinAmont, $$pdebAmont, $s) ;
112 #if ($$pfinAmont==$$pdebAmont) { print "transition 12 -\n"} ;
113 }
114 $$pdebAmont = 0 ; $$pfinAmont = 0 ;
115 last SWITCH ;
116 };
117 /20/ && do { $$pfinIn=($s =~/\+/?$pos-1:$pos+1) ;
118 if (($s =~ /\+/) and ($$pdebIn!=$$pfinIn)) {
119 printf "%s\tsplit\tgene\t%s\t%s\t.\t%s\t.\tName=%s;\n",
120 $seq, $$pdebIn, $$pfinIn, $s, &tagName($seq, $$pdebIn, $$pfinIn, $s) ;
121 } elsif ($$pfinIn!=$$pdebIn) {
122 printf "%s\tsplit\tgene\t%s\t%s\t.\t%s\t.\tName=%s;\n",
123 $seq, $$pfinIn, $$pdebIn, $s, &tagName($seq, $$pfinIn, $$pdebIn, $s) ;
124 }
125 $$pdebIn = 0 ; $$pfinIn = 0 ;
126 last SWITCH ;
127 };
128 /21/ && do { $$pdebAval=$pos ; $$pfinIn=($s =~/\+/?$pos-1:$pos+1) ;
129 if (($s =~ /\+/) and ($$pdebIn!=$$pfinIn)) {
130 printf "%s\tsplit\tgene\t%s\t%s\t.\t%s\t.\tName=%s;\n",
131 $seq, $$pdebIn, $$pfinIn, $s, &tagName($seq, $$pdebIn, $$pfinIn, $s) ;
132 } elsif ($$pfinIn!=$$pdebIn) {
133 printf "%s\tsplit\tgene\t%s\t%s\t.\t%s\t.\tName=%s;\n",
134 $seq, $$pfinIn, $$pdebIn, $s, &tagName($seq, $$pfinIn, $$pdebIn, $s) ;
135 }
136 #$$pdebIn = 0 ; $$pfinIn = 0 ;
137 last SWITCH ;
138 };
139 }
140 }
141 }
142 #-----------------------
143 my ($fileNameI, $fileNameE, $strand) = ("", "", 0) ;
144 # command line check
145 foreach my $num (0 .. $#ARGV) {
146 SWITCH: for ($ARGV[$num]) {
147 /--input1|-i/ && do {
148 $fileNameI=$ARGV[$num+1];
149 open ( fichierGffI, "< $fileNameI" ) or die "Can't open gff file: \"$fileNameI\"\n" ;
150 last };
151 /--input2|-j/ && do {
152 $fileNameE=$ARGV[$num+1];
153 open ( fichierGffE, "< $fileNameE" ) or die "Can't open gff file: \"$fileNameE\"\n" ;
154 last };
155 /--strand|-s/ && do {
156 if ($ARGV[$num+1] eq "s") { $strand=1};
157 last };
158 /--help|-h/ && do { exec("pod2text $0\n") ; die };
159 }
160 }
161 # memory declarations:
162 my @infoI ; my @infoE ;
163 my $seqName ;
164 my @tagF ; my @tagB ; # for Forward and Backward sequences
165 # data retrieval:
166 my @linesI = <fichierGffI> ; my @linesE = <fichierGffE> ;
167 close fichierGffI ; close fichierGffE ;
168 #print "gff files read ; number of lines : ",$#linesI+1," + ",$#linesE+1,"\n";
169 # positions management:
170 &recupInfo(\@infoI, @linesI) ;
171 #print "number of treated elements:",($#infoI+1)/4,"\n";
172 &recupInfo(\@infoE, @linesE) ;
173 #print "number of treated elements:",($#infoE+1)/4,"\n";
174 # treatement:
175 # transform gff lines into chromosomal position tags : 0 for nothing, 1 resp. 2 for extended resp. included elements
176 if (($#infoI) and ($#infoE)) {
177 $seqName=$infoI[0] ;
178 #print "end : $infoE[$#infoE-1]\n";
179 for (my $i=0 ; $i <= $infoE[$#infoE-1] ; $i++) { $tagF[$i] = 0 ; $tagB[$i] = 0 ; } ; # "O" tag in all chr. positions
180 #print "seqName : $seqName\n" ;
181 &feedPositionTab(1, \@tagF, \@tagB, @infoE) ; # "1" tag for all extended elements
182 &feedPositionTab(2, \@tagF, \@tagB, @infoI) ; # "2" tag for all included elements
183 #print join("", @tagF), "\n";
184 #print join("", @tagB), "\n";
185 # transition management:
186 my ($beginUpstream, $endUpstream, $beginIncluded, $endIncluded, $beginDownstream, $endDownstream)
187 = (0, 0, 0, 0, 0, 0) ;
188 for (my $i=1 ; $i <= $#tagF-1 ; $i+=1) {
189 &transitionAnalysis($i, $seqName, "+", \$beginUpstream, \$endUpstream, \$beginIncluded, \$endIncluded, \$beginDownstream, \$endDownstream, \@tagF) ;
190 }
191 ($beginUpstream, $endUpstream, $beginIncluded, $endIncluded, $beginDownstream, $endDownstream)
192 = ($infoE[$#infoE-1], $infoE[$#infoE-1], $infoE[$#infoE-1], $infoE[$#infoE-1], $infoE[$#infoE-1], $infoE[$#infoE-1]) ;
193 for (my $i=$#tagB-1 ; $i >= 1 ; $i-=1) {
194 &transitionAnalysis($i, $seqName, "-", \$beginUpstream, \$endUpstream, \$beginIncluded, \$endIncluded, \$beginDownstream, \$endDownstream, \@tagB) ;
195 }
196 }
197 exit(0) ;