Mercurial > repos > amadeo > amadeo

comparison Tools/Second/step1_galaxy.pl @ 0:229d36377838 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author amadeo
date Mon, 05 Sep 2016 05:53:08 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:229d36377838
1 #!/usr/bin/perl -w
2 $|=1;
3 use warnings;
4 use strict;
5
6 #Script that takes a gff format file from MEME suite as input and orders it by genes,
7 #so it will create an output with all the information grouped by genes. Motifs will be mixed.
8
9 #Declaration of variables
10 my $line;
11 my @cols;
12 my %hash1;
13 my %hash2;
14 my @list1;
15 my @list2;
16 my $gene;
17 my $pos1;
18 my $n;
19 my $index;
20 my $position;
21 my $scalar;
22 my $TF;
23 my $counter=0;#it gives you the number of lines of the gff file. It is a good way to check that the information is not lost.
24
25 #Files that I am going to use
26
27 if(@ARGV < 2){
28 print "\nUsage: step1.pl fimo.gff fimo-position-sorted.gff e\n\n";
29 exit(0);
30 }
31
32 #I open both files, FIMO as the input and OUTPUT as the ouput.
33 open(FIMO, "<$ARGV[0]") ||
34 die "File '$ARGV[0]' not found\n";
35 open(OUTPUT, ">$ARGV[1]") ||
36 die "File '>$ARGV[1]' not found\n";
37
38
39 while (<FIMO>) {
40 $line=$_; #assigning line to variable $line | $_ is a special default variable that here holds the line contents
41 chomp $line; #avoid \n on last field
42 @cols=split; #Splits the string EXPR into a list of strings and returns the list in list context, or the size of the list in scalar context.
43 #This is very useful because the data of the gff file can be called using this variable.
44
45 if ($line=~/^#/){ #prints the first line of the gff file that is different from the rest
46 printf OUTPUT "%s\n", $line;
47 $counter++;
48 }
49 else { #considers the other lines of the file
50 $gene=substr $cols[0],0,21; #variable that returns the name of the gene of the line
51 $pos1 = $cols[3]; #variable that returns the motif's start position on the gene
52 $TF= substr $cols[8],5,8; #variable that returns the name of the motif
53
54 #I use two arrays (list1 and list2) list1 returns the name of the genes and list2 the lines with all the information.
55 #Notice that the gene and its line will have the same position in both list.
56 if (not exists $hash1{$gene}{$TF}{$pos1}) {
57 $hash1{$gene}{$TF}{$pos1}=1;
58 push @list1, $gene;
59 push @list2, $line;
60 }
61
62 }
63
64 }
65
66 #In this section I sort the list1 (genes) by the name of the genes, so I will take the position of every gene sorted
67 #and I will use the position to print out the lines in the order that I want. The main function of this script
68 #is to write the gff file but having the genes sorted by blocks.
69 $n= scalar @list1;
70 my @list_pos_sorted= sort { $list1[$a] cmp $list1[$b] } 0..($n - 1);
71 for (my $i=0; $i <(scalar @list_pos_sorted); $i++){
72 $index=$list_pos_sorted[$i];
73 $position = $list1[$index];
74 #print $hash2{$position};
75 printf OUTPUT "%s\n", $list2[$index];
76 $counter++;
77 }
78 print $counter;