#!/usr/bin/perl -w
use strict;
#############################################################
################### #######################
################### DADES #######################
################### #######################
#############################################################
## ##
## ##
## AUTORES: Natàlia Pérez i Laia Pibernat ##
## ##
## ##
#############################################################
#############################################################
##El programa següent prediu les posicions d'interacció de ##
##factors de transcripció en una seqüència promotora així ##
## com mostra el valor de l'score i el p-value de cada ##
## interacció ##
#############################################################
############# DECLARACIÓ DE VARIABLES ########################
my %motiu = ("A" => [0, 0, 0, 0, 0, 0, 0, 0, 0], #emmagatzema els valors de les matrius dels factors de transcripcio proporcionades
"C" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"G" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"T" => [0, 0, 0, 0, 0, 0, 0, 0, 0]);
my $i = 0; #
my @v; #
my $posicio = 0; #
my $score = 0; #
open(FITXERMEU,"<$ARGV[0]");
while (){ #anem llegint les línies del document, d'una en una.
chomp($_);
if (m/(FA)/) { #expressió regular per a que trobi la primera línia amb l'identificador del factor de transcripció
print "$_\n"; #mostra la primera línia d'informació de cada FT
}
if (m/\d+\t(\d+)\t(\d+)\t(\d+)\t(\d+)/) { #expressió regular per trobar la fila de números
$motiu{"A"}[$i] = $1; # assignem els valors de les variables anònimes de les files a la posició corresponent dels vectors del hash
$motiu{"C"}[$i] = $2; # primera fila que llegirà, posició dels vectors.
$motiu{"G"}[$i] = $3;
$motiu{"T"}[$i] = $4;
$i = $i +1;
}
if ((m/\/\//) && ($i > 2)){ #condició un cop tenim la hash plena
###conversió en matriu de pesos
### càlcul de les proporcions de nucleòtids en la seqüència promotora (pB(x))
open (SEQPROM, "<$ARGV[1]");
my $x = ;
chomp($x);
my @vector = split(//,$x);
my $recorre = 0;
my $pa = 0;
my $pc = 0 ;
my $pt = 0;
my $pg = 0;
my $logpa = 0;
my $logpc = 0;
my $logpt = 0;
my $logpg = 0;
while ($recorre < scalar(@vector)){
if ($vector[$recorre] eq "A" || $vector[$recorre] eq "a") {
$pa = $pa + 1;
}
if( $vector[$recorre] eq "C" || $vector[$recorre] eq "c") {
$pc = $pc + 1;
}
if ($vector[$recorre] eq "T" || $vector[$recorre] eq "t") {
$pt = $pt + 1;
}
if ($vector[$recorre] eq "G" || $vector[$recorre] eq "g") {
$pg = $pg + 1;
}
$recorre = $recorre + 1;
}
$logpa = log($pa/scalar(@vector));
$logpc = log($pc/scalar(@vector));
$logpt = log($pt/scalar(@vector));
$logpg = log($pg/scalar(@vector));
close (SEQPROM);
### càlcul de piM a partir de la hash actual--> obtenim hash de log de freqüències relatives
my %frequencia = ("A" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"C" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"G" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"T" => [0, 0, 0, 0, 0, 0, 0, 0, 0]);
my $h =0;
while ($h < $i){ #recorrem taula hash
if ($motiu{"A"}[$h] == 0){
$frequencia{"A"}[$h]= -999;
}
else {
$frequencia{"A"}[$h] = log ($motiu{"A"}[$h] / ($motiu{"A"}[$h] + $motiu{"C"}[$h] + $motiu{"G"}[$h] + $motiu{"T"}[$h]));
}
if ($motiu{"C"}[$h] == 0){
$frequencia{"C"}[$h]= -999;
}
else{
$frequencia{"C"}[$h] = log ($motiu{"C"}[$h] / ($motiu{"A"}[$h] + $motiu{"C"}[$h] + $motiu{"G"}[$h] + $motiu{"T"}[$h]));
}
if($motiu{"T"}[$h] == 0){
$frequencia{"T"}[$h]= -999;
}
else{
$frequencia{"T"}[$h] = log ($motiu{"T"}[$h] / ($motiu{"A"}[$h] + $motiu{"C"}[$h] + $motiu{"G"}[$h] + $motiu{"T"}[$h]));
}
if($motiu{"G"}[$h] == 0){
$frequencia{"G"}[$h]= -999;
}
else{
$frequencia{"G"}[$h] = log ($motiu{"G"}[$h] / ($motiu{"A"}[$h] + $motiu{"C"}[$h] + $motiu{"G"}[$h] + $motiu{"T"}[$h]));
}
$h = $h +1;
}
#### resta de log, construcció de la matriu de pesos
my %pesos = ("A" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"C" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"G" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"T" => [0, 0, 0, 0, 0, 0, 0, 0, 0]);
my $w =0;
while ($w < $i){ #recorrem taula hash
$pesos{"A"}[$w]= $frequencia{"A"}[$w]- $logpa;
$pesos{"C"}[$w]= $frequencia{"C"}[$w]- $logpc;
$pesos{"T"}[$w]= $frequencia{"T"}[$w]- $logpt;
$pesos{"G"}[$w]= $frequencia{"G"}[$w]- $logpg;
$w = $w +1;
}
## printar la taula hash
my @k = ("A", "C", "G", "T");
my $p = 0;
while ($p < scalar(@k)){
my $j = 0;
print "$k[$p] ";
while ($j < $i ) {
print "\t$pesos{$k[$p]}[$j]";
$j = $j + 1;
}
print "\n";
$p = $p + 1;
}
### anàlisi seq promotora
open (PROMOTER, "<$ARGV[1]");
my $y = ;
@v=split(//,$y);
my $pv = 0;
my @van;
my $millor = -999999;
my $posicio = 0;
while ($pv < (scalar(@v)-$i+1)) {
my $pva = 0;
my $voltes= 0;
while ($voltes < $i ){ ## omplim vector d'anàlisi, serà tan llarg com la longitud del motiu ($i)
$van[$pva] = $v[$pv];
$pv = $pv +1;
$pva = $pva +1;
$voltes = $voltes +1;
}
## càlcul d'score del vector anàlisi
$score = 0;
$pva = 0;
while ($pva < scalar(@van)) {
if ($van[$pva] eq 'A' || $van[$pva] eq 'a' ){
$score = $score + $pesos{"A"}[$pva];
}
if ($van[$pva] eq 'C' || $van[$pva] eq 'c'){
$score = $score + $pesos{"C"}[$pva];
}
if ($van[$pva] eq 'T'|| $van[$pva] eq 't'){
$score = $score + $pesos{"T"}[$pva];
}
if ($van[$pva] eq 'G' || $van[$pva] eq 'g'){
$score = $score + $pesos{"G"}[$pva];
}
$pva = $pva +1 ;
}
if ($score > $millor){
$millor = $score;
$posicio = $pv - $i;
}
$pv = $pv - $i +1;
}
print "posicio: $posicio\n";
print "score: $millor\n";
close (PROMOTER);
###### PART 4: estimar la freqüència d'observar a l'atzar aquest score:
##Algoritme per agafar una seqüència random:
my $finsa100 = 0;
my $c = 0;
while ($finsa100 < 100){
my $n = scalar(@v);
my $ii = $n - 1;
while ($ii >= 0) {
my $jj = int(rand($ii+1));
if ($ii != $jj) {
my $tmp = $v[$ii];
$v[$ii] = $v[$jj];
$v[$jj] = $tmp;
}
$ii = $ii - 1;
}
##Realitzem el score de la seqüència random actual (es farà per les 100 vegades!):
my $pvrandom = 0;
my @vanrandom;
my $millorrandom = -999999;
while ($pvrandom < scalar(@v)-$i+1) {
my $pvarandom=0;
my $voltesrandom= 0;
while ($voltesrandom < $i ){ ## omplim vector d'anàlisi, serà tan llarg com la longitud del motiu ($i)
$vanrandom[$pvarandom] = $v[$pvrandom];
$pvrandom = $pvrandom +1;
$pvarandom = $pvarandom +1;
$voltesrandom = $voltesrandom +1;
}
## càlcul d'score del vector anàlisi
my $scorerandom = 0;
$pvarandom = 0;
while ($pvarandom < scalar(@vanrandom)) {
if ($vanrandom[$pvarandom] eq 'A' || $vanrandom[$pvarandom] eq 'a' ){
$scorerandom = $scorerandom + $pesos{"A"}[$pvarandom];
}
if ($vanrandom[$pvarandom] eq 'C' || $vanrandom[$pvarandom] eq 'c'){
$scorerandom = $scorerandom + $pesos{"C"}[$pvarandom];
}
if ($vanrandom[$pvarandom] eq 'T'|| $vanrandom[$pvarandom] eq 't'){
$scorerandom = $scorerandom + $pesos{"T"}[$pvarandom];
}
if ($vanrandom[$pvarandom] eq 'G' || $vanrandom[$pvarandom] eq 'g'){
$scorerandom = $scorerandom + $pesos{"G"}[$pvarandom];
}
$pvarandom = $pvarandom +1 ;
}
if ($scorerandom > $millorrandom){
$millorrandom = $scorerandom;
}
$pvrandom = $pvrandom - $i +1;
}
if ($millorrandom > $millor){
$c = $c + 1;
}
$finsa100 = $finsa100 + 1;
}
my $pvalue = $c / 100;
print "p-value $pvalue\n";
print "\n";
# deixar la hash a 0.
%motiu = ("A" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"C" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"G" => [0, 0, 0, 0, 0, 0, 0, 0, 0],
"T" => [0, 0, 0, 0, 0, 0, 0, 0, 0]);
$i = 0; #comencem a omplir la nova hash
}
}
close (FITXERMEU);