From 2e9b5aa5d418bdf6911b63ab247bfb4a21ceba51 Mon Sep 17 00:00:00 2001
From: Vincent Danjean <Vincent.Danjean@ens-lyon.org>
Date: Sat, 1 Apr 2006 16:37:43 +0000
Subject: [PATCH] =?UTF-8?q?Qualitatif:=20correction=20bugs,=20d=C3=A9place?=
=?UTF-8?q?ment=20de=20code=20et=20suite=20du=20d=C3=A9veloppement?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
git-svn-id: svn+ssh://imag/users/huron/danjean/svnroot/claire/altree/trunk@153 cf695345-040a-0410-a956-b889e835fe2e
---
ALTree/Input.pm | 3 +-
ALTree/Tree.pm | 57 +++
altree | 443 +-----------------
.../ancestor_absent/association/run_altree | 2 +-
.../ancestor_absent/association/test.res.log | 2 +-
.../ancestor_absent/association/test.tree | 6 +-
to_rewrite.pm | 407 ++++++++++++++++
7 files changed, 489 insertions(+), 431 deletions(-)
create mode 100644 to_rewrite.pm
diff --git a/ALTree/Input.pm b/ALTree/Input.pm
index e785732..9901687 100644
--- a/ALTree/Input.pm
+++ b/ALTree/Input.pm
@@ -32,7 +32,7 @@ use Data::Dumper;
#controle in a hash: $correspondance{$HaploID}->{"case"} and
#$correspondance{$HaploID}->{"controle"}
-sub ReadCorrespond
+sub ReadCorrespondQuali
{
my($name_correspond) =shift;
my($ligne, @tableau);
@@ -98,6 +98,7 @@ sub ReadCorrespond
#}
return(\%correspondance);
}
+#Read the file correspond.txt and put haplotype ID... (to be continued)
sub ReadCorrespondQuanti
{
my($name_correspond) =shift;
diff --git a/ALTree/Tree.pm b/ALTree/Tree.pm
index 331bf68..ced0739 100644
--- a/ALTree/Tree.pm
+++ b/ALTree/Tree.pm
@@ -146,4 +146,61 @@ sub ChangeRoot {
$self->{"root"}=$newroot;
}
+# Niveau 0 : noeud le plus profond sans frère ni oncle avec plusieurs fils
+sub FillLevel
+{
+ my $self=shift;
+
+ my $set_level;
+ $set_level=sub {
+ my($present_node)=shift;
+ my($level)=shift;
+ my($child);
+ my($childlevel)=$level;
+
+ if ($level>0 || $present_node->NbChildren() > 1) {
+ $childlevel++;
+ }
+ my($realchildlevel)=$childlevel;
+ foreach $child (@{$present_node->{"children"}}) {
+ $realchildlevel=$set_level->($child, $childlevel);
+ }
+ $level=$realchildlevel-1;
+ $present_node->{"level"}=$level;
+ if ($level == 0) {
+ $self->{"level0node"}=$present_node;
+ }
+ return $level;
+ };
+
+ $set_level->($self->GetRoot(), 0);
+}
+
+sub FillHeight
+{
+ my $self=shift;
+
+ my $set_height;
+ $set_height=sub {
+ my($present_node)=shift;
+ my($height)=shift;
+ my($child);
+ $height+=1;
+ foreach $child (@{$present_node->{"children"}}) {
+ $set_height->($child, $height);
+ }
+ $present_node->{"height"}=$height;
+ };
+ $set_height->($self->GetRoot(), 0);
+}
+
+sub GetLevel0 {
+ my $self=shift;
+
+ if (! exists($self->{"level0node"})) {
+ $self->FillLevel();
+ }
+ return $self->{"level0node"};
+}
+
1;
diff --git a/altree b/altree
index fa0f092..b9e4aea 100755
--- a/altree
+++ b/altree
@@ -74,6 +74,7 @@ package DataQual;
use constant QUALI => 0;
use constant QUANTI => 1;
+use to_rewrite;
package main;
@@ -235,165 +236,9 @@ sub ClassicalChi2
return ($chi2, $chi2invalid, $error, $sum_control, $sum_case);
}
-sub CalculChi2
-{
- my($tabnodes_a_traiter)=shift;
- my($ddl)=shift;
- my($test_results)=shift;
- my($sign_util)=shift;
- my($chi2, $chi2invalid, $error, $sum_control, $sum_case);
- my($significatif);
- my($p_value);
-
- ($chi2, $chi2invalid, $error, $sum_control, $sum_case)=
- ClassicalChi2($tabnodes_a_traiter);
- if ($error != 0) {
- # TODO: A vérifier : est-ce OK de mettre $significatif à 0
- # la valeur est utilisée au retour de cette fonction
-# $significatif=0;
- if ($error == 1) {
- $test_results->{"texte"}=
- "No cases, ($sum_control controls)";
- if ($sign_util==SignUtil::YES) {
- $significatif=ALTree::Chi2::NON_SIGNIFICATIF;
- }
- } elsif ($error == 2) {
- $test_results->{"texte"}="No controls: only $sum_case cases";
- if ($sum_case>=Seuil::ONLY_CASE) {
- if ($sign_util==SignUtil::YES) {
- $significatif=1;
- $test_results->{"sign"}=ALTree::Chi2::SIGNIFICATIF;
- }
- } else {
- if ($sign_util==SignUtil::YES) {
- $significatif=0;
- $test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
- }
- }
- #$test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
- } elsif ($error == 4) {
- $test_results->{"texte"}="Only one clade";
- if ($sign_util==SignUtil::YES) {
- $significatif=0;
- $test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
- }
- # Manque plein de trucs par rapport à la fonction dans chi2tree...
- } else {
- die "invalid error $error\n";
- }
- } else {
- if ($chi2invalid !=0) {
- $test_results->{"warning"}="Small sample size correction used";
- # J'ai pas compté dans combien de branches...
- if ($ddl == 1) {
- $p_value=ALTree::CUtils::bilateral($tabnodes_a_traiter->[0]->{"case"},
- $tabnodes_a_traiter->[0]->{"control"},
- $tabnodes_a_traiter->[1]->{"case"},
- $tabnodes_a_traiter->[1]->{"control"});
- if ($sign_util==SignUtil::YES) {
- $significatif=ALTree::Chi2::chi2_fisher_significatif($p_value);
- }
- } else {
- my(@clades, $node);
- foreach $node (@{$tabnodes_a_traiter}) {
- push @clades, ($node->{"case"} + $node->{"control"});
- # remplit un tableau contenant les effectifs
- # totaux des différents clades à utiliser dans le
- # réechantillonnage
- }
- ($p_value)= ALTree::Chi2::reech_chi2($sum_case, $sum_control,
- $ddl+1, $chi2, \@clades);
- $test_results->{"warning"}.=" ($p_value)";
- if ($sign_util==SignUtil::YES) {
- $significatif= ALTree::Chi2::reech_significatif ($p_value);
- if ($significatif !=
- ALTree::Chi2::chi2_significatif($ddl, $chi2)) {
- $test_results->{"warning"}.=" Result has changed !";
- }
- }
- }
- } else {
- if ($sign_util==SignUtil::YES) {
- $significatif=ALTree::Chi2::chi2_significatif($ddl, $chi2);
- }
- #my $p=`pochisq $chi2 $ddl`+0; # Verif que les 2 appellent
- #bien la même chose!
- $p_value=ALTree::CUtils::pochisq($chi2,$ddl);
- #if ($p != $p_value) {
- #print STDERR "pochisq: $p != $p_value !\n";
- #}
- }
- if ($sign_util==SignUtil::YES) {
- if ($significatif) {
- $test_results->{"sign"}=ALTree::Chi2::SIGNIFICATIF;
- #$test_results->{"texte"}.="significatif";
- } else {
- $test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
- # $test_results->{"texte"}.="non significatif";
- }
- }
- $test_results->{"chi2"}=$chi2;
- $test_results->{"p_val"}=$p_value;
- #$test_results->{"texte"}.=" [p_value_chi2=$p_value]";
- }
- if ($sign_util==SignUtil::YES) {
- return ($significatif, $p_value);
- } else {
- return ($p_value);
- }
-}
+#REWRITE: sub CalculChi2
-sub parcours_nosplit_chi2split
-{
- my($tabnodes_a_traiter)=shift;
- my($prolonge)=shift;
- my($splitmode)=shift;
- my($node_ecriture)=shift;
- my($sign_util)=shift; # vaut 1 si on a besoin de la significativité, 0 sinon
- my($node, $child, @tab_noeuds_suivants);
- my($val)=0;
- my($test, $p_val);
- my($test_results);
-
- $test_results->{"ddl"}=scalar(@{$tabnodes_a_traiter})-1; # Nb branches -1
- if ($sign_util==SignUtil::YES) {
- ($test, $p_val)=CalculChi2($tabnodes_a_traiter, $test_results->{"ddl"}, $test_results, SignUtil::YES );
- } elsif ($sign_util==SignUtil::NO) {
- ($p_val)=CalculChi2($tabnodes_a_traiter, $test_results->{"ddl"}, $test_results, SignUtil::NO);
- }
- $test_results->{"node_teste"}=$node_ecriture;
- push (@{$node_ecriture->{"res_test"}}, $test_results);
- $test_results->{"level"}=scalar(@{$node_ecriture->{"res_test"}})-1;
-
- if ($sign_util== SignUtil::YES && $test==1 && $splitmode == SplitMode::CHI2SPLIT) { # sign et que on on est en chi2split
- foreach $node (@{$tabnodes_a_traiter}) {
- if (NbFils($node) != 0) {
- my @children=$node->GetChildrenList();
- parcours_nosplit_chi2split(\@children,
- $prolonge, $splitmode, $node);
- }
- }
- } elsif ($sign_util== SignUtil::NO || $test==0 || $splitmode == SplitMode::NOSPLIT) { # ou alors on est en nosplit
- foreach $node (@{$tabnodes_a_traiter}) {
- if (NbFils($node) != 0) {
- $val=1;
- foreach $child ($node->GetChildrenList()) {
- push (@tab_noeuds_suivants, $child);
- }
- } else {
- if ($prolonge == 1) {
- push (@tab_noeuds_suivants, $node);
- }
- }
- }
- if ($val==1) {
- parcours_nosplit_chi2split(\@tab_noeuds_suivants,
- $prolonge, $splitmode, $node_ecriture, $sign_util);
- } else {
- return;
- }
- }
-}
+#REWRITE: sub parcours_nosplit_chi2split
sub FillCaseControl
{
@@ -472,159 +317,9 @@ sub FillQuanti
}
-sub ParcoursQuanti
-{
- my($tabnodes_a_traiter)=shift;
- my($prolonge)=shift;
- my($splitmode)=shift;
- my($node_ecriture)=shift;
- my($sign_util)=shift; # vaut 1 si on a besoin de la significativité, 0 sinon
- my($node, $child, @tab_noeuds_suivants);
- my($val)=0;
- my($test, $res_anova);
- my($test_results);
-
-# $test_results->{"ddl"}=scalar(@{$tabnodes_a_traiter})-1; # Nb branches -1
- my @valeurs;
- my @facteurs;
- my $i=0;
-#DEBUG print STDERR "TTTT ", scalar (@{$tabnodes_a_traiter}), "\n";
- foreach $node (@{$tabnodes_a_traiter}) {
- $i++;
- foreach my $case (@{$node->{"quanti"}}) {
- push (@valeurs, $case->[0]);
- push (@facteurs, $i);
- }
- }
- my $nb_factors=$i;
- $test_results->{"nb_facteurs"}=$nb_factors;
-# DEBUG print STDERR "node ";
-# for (my $i=0; $i<=$#valeurs; $i++) {
-# print STDERR " $valeurs[$i]";
-# print STDERR " ($facteurs[$i])";
-# }
-# print STDERR "\n";
-
-
- if ($sign_util==SignUtil::YES) {
- ($test, $res_anova)=CalculAnovaOneWay($tabnodes_a_traiter, \@valeurs, \@facteurs, $test_results, SignUtil::YES, $nb_factors );
- } elsif ($sign_util==SignUtil::NO) {
- ($res_anova)=CalculAnovaOneWay($tabnodes_a_traiter, \@valeurs, \@facteurs, $test_results, SignUtil::NO, $nb_factors);
- }
- $test_results->{"node_teste"}=$node_ecriture;
- push (@{$node_ecriture->{"res_test"}}, $test_results);
- $test_results->{"level"}=scalar(@{$node_ecriture->{"res_test"}})-1;
-
- if ($sign_util== SignUtil::YES && $test==1 && $splitmode == SplitMode::CHI2SPLIT) { # sign et que on on est en chi2split
- foreach $node (@{$tabnodes_a_traiter}) {
- if (NbFils($node) != 0) {
- my @children=$node->GetChildrenList();
- ParcoursQuanti(\@children,
- $prolonge, $splitmode, $node);
- }
- }
- } elsif ($sign_util== SignUtil::NO || $test==0 || $splitmode == SplitMode::NOSPLIT) { # ou alors on est en nosplit
- foreach $node (@{$tabnodes_a_traiter}) {
- if (NbFils($node) != 0) {
- $val=1;
- foreach $child ($node->GetChildrenList()) {
- push (@tab_noeuds_suivants, $child);
- }
- } else {
- if ($prolonge == 1) {
- push (@tab_noeuds_suivants, $node);
- }
- }
- }
- if ($val==1) {
- ParcoursQuanti(\@tab_noeuds_suivants,
- $prolonge, $splitmode, $node_ecriture, $sign_util);
- } else {
- return;
- }
- }
-}
-
-sub CalculAnovaOneWay
-{
- my $tabnodes_a_traiter=shift;
- my $valeurs=shift;
- my $facteurs=shift;
- my $test_results =shift;
- my $sign_util=shift;
- my $nb_factors=shift;
-
- my $significatif;
- my $res_anova;
- if (scalar (@{$tabnodes_a_traiter}) < 2) {
- $test_results->{"texte"}=
- "Only one category";
- if ($sign_util==SignUtil::YES) {
- $significatif=ALTree::Chi2::NON_SIGNIFICATIF;
- }
- } else {
- if (scalar (@{$valeurs}) != scalar (@{$facteurs})) {
- erreur("Error in the anova data: the number of values ", scalar @${valeurs}, " and the number of factors ", scalar @{$facteurs}, " should be the same\n");
- } else {
- $res_anova=TamuAnova::anova($valeurs, $facteurs, $nb_factors);
- #DEBUG print STDERR $nb_factors, " ", scalar(@{$valeurs}), "\n";
- $test_results->{"F"}=$res_anova->{"F"};
- $test_results->{"p_val"}= $res_anova->{"p"};
- # $test_results->{"warning"}.=" ($p_value)";
- if ($sign_util==SignUtil::YES) {
- $significatif = ALTree::Chi2::chi2_fisher_significatif($res_anova->{"p"});
- }
-
- }
- }
- if ($sign_util==SignUtil::YES) {
- if ($significatif) {
- $test_results->{"sign"}=ALTree::Chi2::SIGNIFICATIF;
- #$test_results->{"texte"}.="significatif";
- } else {
- $test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
- # $test_results->{"texte"}.="non significatif";
- }
- }
- if ($sign_util==SignUtil::YES) {
- return ($significatif, $res_anova);
- } else {
- return ($res_anova);
- }
-
-#DEBUG if (defined ($res_anova)) {
-#DEBUG print STDERR "RESANOVA ", $res_anova->{"F"}, " ", $res_anova->{"p"}, "\n";
-}
-
-
-
-########## GENERAL ##########
-
-sub FillLevel
-{
- my($present_node)=shift;
- my($level)=shift;
- my($child);
-
- $present_node->{"level"}=$level;
- $level+=1;
- foreach $child (@{$present_node->{"children"}}) {
- FillLevel($child, $level);
- }
-}
-
-sub FillHeight
-{
- my($present_node)=shift;
- my($height)=shift;
- my($child);
- $height+=1;
- foreach $child (@{$present_node->{"children"}}) {
- FillHeight($child, $height);
- }
- $present_node->{"height"}=$height;
-}
+#REWRITE: sub ParcoursQuanti
+#REWRITE: sub CalculAnovaOneWay
##################################################
######## AFFICHAGE DE L' ARBRE ##################
@@ -771,113 +466,9 @@ sub InfosQuantiNoperm
return InfosAffichees($node, 4);
}
#Return ddl, level, pvalues and chi2
-sub InfosAffichees
-{
- my($node)=shift;
- my($mode)=shift;
- my($chaine)=Name($node);
- my($lbl_test)=0;
- my $test;
-
- if ($mode==1 || $mode == 2) { # Affiche ou pas les case/control
- $chaine.=" case/control:".$node->{"case"}."/".$node->{"control"};
- }
- if ($mode==3 || $mode == 4) {
- $chaine.= sprintf " mean:%.2f",$node->{"moyenne"};
- }
- if (1) { # affiche les apomorphies
- $chaine.="\n";
- foreach my $apo ($node->GetApoList()) {
- $chaine.= (" Site: ".$apo->GetSiteNb." Sens: ".$apo->GetSensLabel()."\n");
- }
- }
- $chaine.="\n";
- if ($mode==1 || $mode == 2) { # affiche ou pas les ddl
- if (defined $node->{"res_test"}) {
- for $test (@{$node->{"res_test"}}) {
- $chaine.= sprintf "[%d] ddl=%d",
- $test->{"level"}, $test->{"ddl"};
- if ($test->{"ddl"} > 0) {
- $chaine.= sprintf " chi2=%.2f p_value_chi2=%.3g",
- $test->{"chi2"}, $test->{"p_val"};
- # TODO : ça arrive quand on a que des malades ou témoins
- # dans les clades...
- if (not defined($test->{"chi2"})) {
- print "chi2 for ", Name($node),
- "(", $test->{"ddl"}, ")", "\n";
- }
- if (not defined($test->{"p_val"})) {
- print "p_val for ", Name($node),
- "(", $test->{"ddl"}, ")", "\n";
- }
-
- if ($mode ==2) {
- if (defined($test->{"sign"})) {
- if ($test->{"sign"} == ALTree::Chi2::NON_SIGNIFICATIF) {
- $chaine .= " (non significatif)";
- } elsif ($test->{"sign"} == ALTree::Chi2::SIGNIFICATIF) {
- $chaine .= " (significatif)";
- } else {
- ALTree::Utils::internal_error("unknown value ".
- $test->{"sign"});
- }
- }
- if (defined($test->{"texte"})) {
- $chaine .= "\n".$test->{"texte"};
- }
- if (defined($test->{"warning"})) {
- $chaine .= "\n".$test->{"warning"};
- }
- }
- }
- $chaine.="\n";
- }
- }
- } elsif ($mode == 3 || $mode ==4) {
- if (defined $node->{"res_test"}) {
- for $test (@{$node->{"res_test"}}) {
- $chaine.= sprintf "[%d] nb_fact=%d",
- $test->{"level"}, $test->{"nb_facteurs"};
- if ($test->{"nb_facteurs"} > 1) {
- $chaine.= sprintf " F=%.2f p_value=%.3g",
- $test->{"F"}, $test->{"p_val"};
- # TODO : ça arrive quand on a que des malades ou témoins
- # dans les clades...
- if (not defined($test->{"F"})) {
- print "F for ", Name($node),
- "(", $test->{"nb_facteurs"}, ")", "\n";
- }
- if (not defined($test->{"p_val"})) {
- print "p_val for ", Name($node),
- "(", $test->{"nb_facteurs"}, ")", "\n";
- }
-
- if ($mode == 4) {
- if (defined($test->{"sign"})) {
- if ($test->{"sign"} == ALTree::Chi2::NON_SIGNIFICATIF) {
- $chaine .= " (non significatif)";
- } elsif ($test->{"sign"} == ALTree::Chi2::SIGNIFICATIF) {
- $chaine .= " (significatif)";
- } else {
- ALTree::Utils::internal_error("unknown value ".
- $test->{"sign"});
- }
- }
- if (defined($test->{"texte"})) {
- $chaine .= "\n".$test->{"texte"};
- }
- if (defined($test->{"warning"})) {
- $chaine .= "\n".$test->{"warning"};
- }
- }
- }
- $chaine.="\n";
- }
- }
- }
- return($chaine);
-
-}
+
+#REWRITE: sub InfosAffichees
+
##########################################################
######## MODIFICATIONS/CALCULS SUR L'ARBRE ###############
##########################################################
@@ -1042,7 +633,7 @@ sub Association
SplitMode::NOSPLIT, $racine, $sign_util);
}
-sub RepeatAssociation
+sub RepeatAssociationQuali
{
my($tree)=shift;
my($correspondance)=shift;
@@ -1820,7 +1411,7 @@ sub main
# or to a tab containing the different quantitative values
if ($dataqual == DataQual::QUALI) {
- $correspondance=ALTree::Input::ReadCorrespond($name_corres);
+ $correspondance=ALTree::Input::ReadCorrespondQuali($name_corres);
} elsif ($dataqual == DataQual::QUANTI) {
$correspondance=ALTree::Input::ReadCorrespondQuanti($name_corres);
} else {
@@ -1850,8 +1441,8 @@ sub main
#print "\n";
FusionBrNulles($tree->GetRoot());
# Structure changée, on recalcul...
- FillHeight($tree->GetRoot(), 0);
- FillLevel($tree->GetRoot(), 0);
+ $tree->FillHeight();
+ $tree->FillLevel();
if ($dataqual == DataQual::QUALI) {
FillCaseControl($tree->GetRoot(),$correspondance);
@@ -1883,9 +1474,11 @@ sub main
die "invalid value for the number of permutation: $permutation\n";
}
if ($dataqual == DataQual::QUALI) {
- parcours_nosplit_chi2split(\@children, $prolonge, $splitmode, $racine, $sign_util );
+ parcours_nosplit_chi2split(\@children, $prolonge, $splitmode,
+ $racine, $sign_util);
} else {
- ParcoursQuanti(\@children, $prolonge, $splitmode, $racine, $sign_util );
+ ParcoursQuanti(\@children, $prolonge, $splitmode,
+ $racine, $sign_util);
}
{
@@ -1899,13 +1492,13 @@ sub main
my($value_per_line, $ligne_test);
if ($dataqual == DataQual::QUALI) {
AffichageArbre($racine, \&AssociationInfos);
- ($value_per_line, $ligne_test)=RepeatAssociation
+ ($value_per_line, $ligne_test)=RepeatAssociationQuali
($tree, $correspondance, $prolonge,$permutation, $sign_util);
} else {
AffichageArbre($racine, \&InfosQuanti);
($value_per_line, $ligne_test)=RepeatAssociationQuanti
($tree, $correspondance, $prolonge,$permutation, $sign_util);
- }
+ }
my($corrected_values);
$corrected_values=ALTree::CUtils::double_permutation
($permutation+1, $value_per_line, $ligne_test);
diff --git a/test/paup/ancestor_absent/association/run_altree b/test/paup/ancestor_absent/association/run_altree
index 1549d33..6002ab9 100755
--- a/test/paup/ancestor_absent/association/run_altree
+++ b/test/paup/ancestor_absent/association/run_altree
@@ -5,4 +5,4 @@ paup caco.paup
#To perform the association test
../../../../altree -i test.res.log -j nb_cas_control.txt -a -t SNP \
- -p paup -r 1 --tree-to-analyse 1 -o 1_caco.asso
+ -p paup -r 1 --tree-to-analyse 1 -o 1_caco.asso --data-qual qualitative
diff --git a/test/paup/ancestor_absent/association/test.res.log b/test/paup/ancestor_absent/association/test.res.log
index 48fd392..fec2757 100644
--- a/test/paup/ancestor_absent/association/test.res.log
+++ b/test/paup/ancestor_absent/association/test.res.log
@@ -1,7 +1,7 @@
P A U P *
Portable version 4.0b10 for Unix
-Thu Jan 12 21:56:20 2006
+Sat Apr 1 16:01:41 2006
-----------------------------NOTICE-----------------------------
This is a beta-test version. Please report any crashes,
diff --git a/test/paup/ancestor_absent/association/test.tree b/test/paup/ancestor_absent/association/test.tree
index 2398457..2da5b28 100644
--- a/test/paup/ancestor_absent/association/test.tree
+++ b/test/paup/ancestor_absent/association/test.tree
@@ -1,8 +1,8 @@
#NEXUS
-Begin trees; [Treefile saved Thu Jan 12 21:56:20 2006]
+Begin trees; [Treefile saved Sat Apr 1 16:01:41 2006]
[!
->Data file = /home/cbardel/recherche/logiciel/altree/test/paup/ancestor_absent/association/caco.paup
+>Data file = /home/vdanjean/travail/perso/Claire/svn/altree/test/paup/ancestor_absent/association/caco.paup
>Heuristic search settings:
> Optimality criterion = parsimony
> Character-status summary:
@@ -26,7 +26,7 @@ Begin trees; [Treefile saved Thu Jan 12 21:56:20 2006]
> Total number of rearrangements tried = 69978
> Score of best tree(s) found = 20
> Number of trees retained = 89
-> Time used = 1 sec (CPU time = 0.09 sec)
+> Time used = 1 sec (CPU time = 0.04 sec)
]
tree PAUP_1 = [&R] (((((H002,(H008,H014),H009,H006),H007,H012,H011),H013),H001,(H003,H005)),H010);
tree PAUP_2 = [&R] (((((((H002,(H007,H012,H011),H009,H006),H008),H014),H013),H001),(H003,H005)),H010);
diff --git a/to_rewrite.pm b/to_rewrite.pm
new file mode 100644
index 0000000..5d1a89f
--- /dev/null
+++ b/to_rewrite.pm
@@ -0,0 +1,407 @@
+
+package main;
+use strict;
+use diagnostics;
+use warnings;
+
+use ALTree::Chi2 ();
+use ALTree::Import;
+use ALTree::Utils qw(erreur);
+use ALTree::Input qw(PrepareTree);
+#use Newchi2treeUtils;
+use TamuAnova;
+
+sub parcours_nosplit_chi2split
+{
+ my($tabnodes_a_traiter)=shift;
+ my($prolonge)=shift;
+ my($splitmode)=shift;
+ my($node_ecriture)=shift;
+ my($sign_util)=shift; # vaut 1 si on a besoin de la significativité, 0 sinon
+ my($node, $child, @tab_noeuds_suivants);
+ my($val)=0;
+ my($test, $p_val);
+ my($test_results);
+
+ $test_results->{"ddl"}=scalar(@{$tabnodes_a_traiter})-1; # Nb branches -1
+ # $test n'est valable que si $sign_util est à YES
+ ($p_val, $test)=CalculChi2($tabnodes_a_traiter, $test_results->{"ddl"},
+ $test_results, $sign_util );
+ $test_results->{"node_teste"}=$node_ecriture;
+ push (@{$node_ecriture->{"res_test"}}, $test_results);
+ $test_results->{"level"}=scalar(@{$node_ecriture->{"res_test"}})-1;
+
+ if ($sign_util== SignUtil::YES && $test==1 && $splitmode == SplitMode::CHI2SPLIT) { # sign et que on on est en chi2split
+ foreach $node (@{$tabnodes_a_traiter}) {
+ if (NbFils($node) != 0) {
+ my @children=$node->GetChildrenList();
+ parcours_nosplit_chi2split(\@children,
+ $prolonge, $splitmode, $node);
+ }
+ }
+ } elsif ($sign_util== SignUtil::NO || $test==0 || $splitmode == SplitMode::NOSPLIT) { # ou alors on est en nosplit
+ foreach $node (@{$tabnodes_a_traiter}) {
+ if (NbFils($node) != 0) {
+ $val=1;
+ foreach $child ($node->GetChildrenList()) {
+ push (@tab_noeuds_suivants, $child);
+ }
+ } else {
+ if ($prolonge == 1) {
+ push (@tab_noeuds_suivants, $node);
+ }
+ }
+ }
+ if ($val==1) {
+ parcours_nosplit_chi2split(\@tab_noeuds_suivants,
+ $prolonge, $splitmode, $node_ecriture, $sign_util);
+ } else {
+ return;
+ }
+ }
+}
+
+sub ParcoursQuanti
+{
+ my($tabnodes_a_traiter)=shift;
+ my($prolonge)=shift;
+ my($splitmode)=shift;
+ my($node_ecriture)=shift;
+ my($sign_util)=shift; # vaut 1 si on a besoin de la significativité, 0 sinon
+ my($node, $child, @tab_noeuds_suivants);
+ my($val)=0;
+ my($test, $res_anova);
+ my($test_results);
+
+# $test_results->{"ddl"}=scalar(@{$tabnodes_a_traiter})-1; # Nb branches -1
+ my @valeurs;
+ my @facteurs;
+ my $i=0;
+#DEBUG print STDERR "TTTT ", scalar (@{$tabnodes_a_traiter}), "\n";
+ foreach $node (@{$tabnodes_a_traiter}) {
+ $i++;
+ foreach my $case (@{$node->{"quanti"}}) {
+ push (@valeurs, $case->[0]);
+ push (@facteurs, $i);
+ }
+ }
+ my $nb_factors=$i;
+ $test_results->{"nb_facteurs"}=$nb_factors;
+# DEBUG print STDERR "node ";
+# for (my $i=0; $i<=$#valeurs; $i++) {
+# print STDERR " $valeurs[$i]";
+# print STDERR " ($facteurs[$i])";
+# }
+# print STDERR "\n";
+
+
+# if ($sign_util==SignUtil::YES) {
+ ($res_anova, $test)=CalculAnovaOneWay($tabnodes_a_traiter, \@valeurs, \@facteurs, $test_results, $sign_util, $nb_factors );
+# } elsif ($sign_util==SignUtil::NO) {
+# ($res_anova)=CalculAnovaOneWay($tabnodes_a_traiter, \@valeurs, \@facteurs, $test_results, $sign_util, $nb_factors);
+# }
+ $test_results->{"node_teste"}=$node_ecriture;
+ push (@{$node_ecriture->{"res_test"}}, $test_results);
+ $test_results->{"level"}=scalar(@{$node_ecriture->{"res_test"}})-1;
+
+ if ($sign_util== SignUtil::YES && $test==ALTree::Chi2::SIGNIFICATIF
+ && $splitmode == SplitMode::CHI2SPLIT) {
+ # sign et que on on est en chi2split
+ foreach $node (@{$tabnodes_a_traiter}) {
+ if (NbFils($node) != 0) {
+ my @children=$node->GetChildrenList();
+ ParcoursQuanti(\@children,
+ $prolonge, $splitmode, $node);
+ }
+ }
+ } elsif ($sign_util== SignUtil::NO || $test==ALTree::Chi2::NON_SIGNIFICATIF
+ || $splitmode == SplitMode::NOSPLIT) {
+ # ou alors on est en nosplit
+ foreach $node (@{$tabnodes_a_traiter}) {
+ if (NbFils($node) != 0) {
+ $val=1;
+ foreach $child ($node->GetChildrenList()) {
+ push (@tab_noeuds_suivants, $child);
+ }
+ } else {
+ if ($prolonge == 1) {
+ push (@tab_noeuds_suivants, $node);
+ }
+ }
+ }
+ if ($val==1) {
+ ParcoursQuanti(\@tab_noeuds_suivants,
+ $prolonge, $splitmode, $node_ecriture, $sign_util);
+ } else {
+ return;
+ }
+ } else {
+ die("Arghhh");
+ }
+}
+
+sub CalculChi2
+{
+ my($tabnodes_a_traiter)=shift;
+ my($ddl)=shift;
+ my($test_results)=shift;
+ my($sign_util)=shift;
+ my($chi2, $chi2invalid, $error, $sum_control, $sum_case);
+ my($significatif)=ALTree::Chi2::NON_SIGNIFICATIF;
+ my($p_value);
+
+ ($chi2, $chi2invalid, $error, $sum_control, $sum_case)=
+ ClassicalChi2($tabnodes_a_traiter);
+ if ($error != 0) {
+ # TODO: A vérifier : est-ce OK de mettre $significatif à 0
+ # la valeur est utilisée au retour de cette fonction
+# $significatif=0;
+ if ($error == 1) {
+ $test_results->{"texte"}=
+ "No cases, ($sum_control controls)";
+ if ($sign_util==SignUtil::YES) {
+ $significatif=ALTree::Chi2::NON_SIGNIFICATIF;
+ }
+ } elsif ($error == 2) {
+ $test_results->{"texte"}="No controls: only $sum_case cases";
+ if ($sum_case>=Seuil::ONLY_CASE) {
+ if ($sign_util==SignUtil::YES) {
+ $significatif=1;
+ $test_results->{"sign"}=ALTree::Chi2::SIGNIFICATIF;
+ }
+ } else {
+ if ($sign_util==SignUtil::YES) {
+ $significatif=0;
+ $test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
+ }
+ }
+ #$test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
+ } elsif ($error == 4) {
+ $test_results->{"texte"}="Only one clade";
+ if ($sign_util==SignUtil::YES) {
+ $significatif=0;
+ $test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
+ }
+ # Manque plein de trucs par rapport à la fonction dans chi2tree...
+ } else {
+ die "invalid error $error\n";
+ }
+ } else {
+ if ($chi2invalid !=0) {
+ $test_results->{"warning"}="Small sample size correction used";
+ # J'ai pas compté dans combien de branches...
+ if ($ddl == 1) {
+ $p_value=ALTree::CUtils::bilateral($tabnodes_a_traiter->[0]->{"case"},
+ $tabnodes_a_traiter->[0]->{"control"},
+ $tabnodes_a_traiter->[1]->{"case"},
+ $tabnodes_a_traiter->[1]->{"control"});
+ if ($sign_util==SignUtil::YES) {
+ $significatif=ALTree::Chi2::chi2_fisher_significatif($p_value);
+ }
+ } else {
+ my(@clades, $node);
+ foreach $node (@{$tabnodes_a_traiter}) {
+ push @clades, ($node->{"case"} + $node->{"control"});
+ # remplit un tableau contenant les effectifs
+ # totaux des différents clades à utiliser dans le
+ # réechantillonnage
+ }
+ ($p_value)= ALTree::Chi2::reech_chi2($sum_case, $sum_control,
+ $ddl+1, $chi2, \@clades);
+ $test_results->{"warning"}.=" ($p_value)";
+ if ($sign_util==SignUtil::YES) {
+ $significatif= ALTree::Chi2::reech_significatif ($p_value);
+ if ($significatif !=
+ ALTree::Chi2::chi2_significatif($ddl, $chi2)) {
+ $test_results->{"warning"}.=" Result has changed !";
+ }
+ }
+ }
+ } else {
+ if ($sign_util==SignUtil::YES) {
+ $significatif=ALTree::Chi2::chi2_significatif($ddl, $chi2);
+ }
+ #my $p=`pochisq $chi2 $ddl`+0; # Verif que les 2 appellent
+ #bien la même chose!
+ $p_value=ALTree::CUtils::pochisq($chi2,$ddl);
+ #if ($p != $p_value) {
+ #print STDERR "pochisq: $p != $p_value !\n";
+ #}
+ }
+ if ($sign_util==SignUtil::YES) {
+ if ($significatif) {
+ $test_results->{"sign"}=ALTree::Chi2::SIGNIFICATIF;
+ #$test_results->{"texte"}.="significatif";
+ } else {
+ $test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
+ # $test_results->{"texte"}.="non significatif";
+ }
+ }
+ $test_results->{"chi2"}=$chi2;
+ $test_results->{"p_val"}=$p_value;
+ #$test_results->{"texte"}.=" [p_value_chi2=$p_value]";
+ }
+# if ($sign_util==SignUtil::YES) {
+ return ($p_value, $significatif);
+# } else {
+# return ($p_value);
+# }
+}
+
+sub CalculAnovaOneWay
+{
+ my $tabnodes_a_traiter=shift;
+ my $valeurs=shift;
+ my $facteurs=shift;
+ my $test_results =shift;
+ my $sign_util=shift;
+ my $nb_factors=shift;
+
+ my $significatif=ALTree::Chi2::NON_SIGNIFICATIF;
+ my $res_anova;
+ if (scalar (@{$tabnodes_a_traiter}) < 2) {
+ $test_results->{"texte"}=
+ "Only one category";
+ if ($sign_util==SignUtil::YES) {
+ $significatif=ALTree::Chi2::NON_SIGNIFICATIF;
+ }
+ } else {
+ if (scalar (@{$valeurs}) != scalar (@{$facteurs})) {
+ erreur("Error in the anova data: the number of values ", scalar @${valeurs}, " and the number of factors ", scalar @{$facteurs}, " should be the same\n");
+ } else {
+ $res_anova=TamuAnova::anova($valeurs, $facteurs, $nb_factors);
+ #DEBUG print STDERR $nb_factors, " ", scalar(@{$valeurs}), "\n";
+ $test_results->{"F"}=$res_anova->{"F"};
+ $test_results->{"p_val"}= $res_anova->{"p"};
+ # $test_results->{"warning"}.=" ($p_value)";
+ if ($sign_util==SignUtil::YES) {
+ $significatif = ALTree::Chi2::chi2_fisher_significatif($res_anova->{"p"});
+ }
+
+ }
+ }
+ if ($sign_util==SignUtil::YES) {
+ if ($significatif) {
+ $test_results->{"sign"}=ALTree::Chi2::SIGNIFICATIF;
+ #$test_results->{"texte"}.="significatif";
+ } else {
+ $test_results->{"sign"}=ALTree::Chi2::NON_SIGNIFICATIF;
+ # $test_results->{"texte"}.="non significatif";
+ }
+ }
+ return ($res_anova, $significatif);
+
+#DEBUG if (defined ($res_anova)) {
+#DEBUG print STDERR "RESANOVA ", $res_anova->{"F"}, " ", $res_anova->{"p"}, "\n";
+}
+
+sub InfosAffichees
+{
+ my($node)=shift;
+ my($mode)=shift;
+ my($chaine)=Name($node);
+ my($lbl_test)=0;
+ my $test;
+
+ $chaine.=" (LEVEL: ".$node->{"level"}.")";
+ if ($mode==1 || $mode == 2) { # Affiche ou pas les case/control
+ $chaine.=" case/control:".$node->{"case"}."/".$node->{"control"};
+ }
+ if ($mode==3 || $mode == 4) {
+ $chaine.= sprintf " mean:%.2f",$node->{"moyenne"};
+ }
+ if (1) { # affiche les apomorphies
+ $chaine.="\n";
+ foreach my $apo ($node->GetApoList()) {
+ $chaine.= (" Site: ".$apo->GetSiteNb." Sens: ".$apo->GetSensLabel()."\n");
+ }
+ }
+ $chaine.="\n";
+ if ($mode==1 || $mode == 2) { # affiche ou pas les ddl
+ if (defined $node->{"res_test"}) {
+ for $test (@{$node->{"res_test"}}) {
+ $chaine.= sprintf "[%d] ddl=%d",
+ $test->{"level"}, $test->{"ddl"};
+ if ($test->{"ddl"} > 0) {
+ $chaine.= sprintf " chi2=%.2f p_value_chi2=%.3g",
+ $test->{"chi2"}, $test->{"p_val"};
+ # TODO : ça arrive quand on a que des malades ou témoins
+ # dans les clades...
+ if (not defined($test->{"chi2"})) {
+ print "chi2 for ", Name($node),
+ "(", $test->{"ddl"}, ")", "\n";
+ }
+ if (not defined($test->{"p_val"})) {
+ print "p_val for ", Name($node),
+ "(", $test->{"ddl"}, ")", "\n";
+ }
+
+ if ($mode ==2) {
+ if (defined($test->{"sign"})) {
+ if ($test->{"sign"} == ALTree::Chi2::NON_SIGNIFICATIF) {
+ $chaine .= " (non significatif)";
+ } elsif ($test->{"sign"} == ALTree::Chi2::SIGNIFICATIF) {
+ $chaine .= " (significatif)";
+ } else {
+ ALTree::Utils::internal_error("unknown value ".
+ $test->{"sign"});
+ }
+ }
+ if (defined($test->{"texte"})) {
+ $chaine .= "\n".$test->{"texte"};
+ }
+ if (defined($test->{"warning"})) {
+ $chaine .= "\n".$test->{"warning"};
+ }
+ }
+ }
+ $chaine.="\n";
+ }
+ }
+ } elsif ($mode == 3 || $mode ==4) {
+ if (defined $node->{"res_test"}) {
+ for $test (@{$node->{"res_test"}}) {
+ $chaine.= sprintf "[%d] nb_fact=%d",
+ $test->{"level"}, $test->{"nb_facteurs"};
+ if ($test->{"nb_facteurs"} > 1) {
+ $chaine.= sprintf " F=%.2f p_value=%.3g",
+ $test->{"F"}, $test->{"p_val"};
+ # TODO : ça arrive quand on a que des malades ou témoins
+ # dans les clades...
+ if (not defined($test->{"F"})) {
+ print "F for ", Name($node),
+ "(", $test->{"nb_facteurs"}, ")", "\n";
+ }
+ if (not defined($test->{"p_val"})) {
+ print "p_val for ", Name($node),
+ "(", $test->{"nb_facteurs"}, ")", "\n";
+ }
+
+ if ($mode == 4) {
+ if (defined($test->{"sign"})) {
+ if ($test->{"sign"} == ALTree::Chi2::NON_SIGNIFICATIF) {
+ $chaine .= " (non significatif)";
+ } elsif ($test->{"sign"} == ALTree::Chi2::SIGNIFICATIF) {
+ $chaine .= " (significatif)";
+ } else {
+ ALTree::Utils::internal_error("unknown value ".
+ $test->{"sign"});
+ }
+ }
+ if (defined($test->{"texte"})) {
+ $chaine .= "\n".$test->{"texte"};
+ }
+ if (defined($test->{"warning"})) {
+ $chaine .= "\n".$test->{"warning"};
+ }
+ }
+ }
+ $chaine.="\n";
+ }
+ }
+ }
+ return($chaine);
+
+}
+
+1;
--
GitLab