altree 1.3.1-1 / usr / bin / altree

#!/usr/bin/perl

eval 'exec /usr/bin/perl  -S $0 ${1+"$@"}'
    if 0; # not running under some shell

use strict;
use diagnostics;
use warnings;
use Getopt::Long qw(:config permute);
use Pod::Usage;

use Data::Dumper;

srand(654321);

use version; our $VERSION = qv(1.3.1);

# Positionner la variable PERL5LIB si besoin

# PERL5LIB est une variable comme PATH, sauf qu'elle ne sert pas à
# trouver les programmes, mais plutôt les modules de perl
# En bash :
# export PERL5LIB=/chemin/vers/modules/perl:/chemin/vers/autres/modules
# En tcsh :
# setenv PERL5LIB /chemin/vers/modules/perl:/chemin/vers/autres/modules

# Si les modules sont installés dans les emplacements standard de
# perl, c'est inutile

use ALTree::Chi2 ();
use ALTree::Import;
use ALTree::Utils qw(erreur);
use ALTree::Input qw(PrepareTree);
#use Newchi2treeUtils;
use Math::TamuAnova;
#use NAnova;
use ALTree::Nanova;

###########################################
########  GLOBAL VARIABLES        #########
###########################################

# Variable $nodes
#   Ref on Hash of ('id' => Node)
# my $nodes;

# Variable $sites
#   Ref on Hash of ('site_nb' -> Site)
# my $sites;

###########################################
########  CONSTANTES              #########
###########################################

package SplitMode;
use constant NOSPLIT   => 0;
use constant CHI2SPLIT => 1;

package CoEvo;
use constant SIMPLE   => 0;
use constant DOUBLE => 1;

package RootMeth;
use constant OUTG  => 0;
use constant ANC   => 1;

package SignUtil;
use constant NO   => 0;
use constant YES  => 1;

package Seuil;
use constant SAMPLESIZE => 5;
use constant P_VAL_CHI2 => 0.01;
use constant P_VAL_TESTPROP => 0.01;
use constant ONLY_CASE  => 3;

package DataQual;
use constant QUALI   => 0;
use constant QUANTI => 1;

use ALTree::to_rewrite;

package main;

###########################################
#########  BUILDING OF THE TREE  ##########
###########################################


# Outgroup is not removed from the hash nodes.
# It is only removed from the list of children of it's father 
sub RemoveOutgroup
{
    my $tree=shift;
    my $outgroup=shift;

    my($father_outgr)=$outgroup->GetFather();
    $father_outgr->DeleteChild($outgroup);
}

# Outgroup is put again in the list of children of it's father
sub AddOutgroup 
{
    my($outgroup)=shift;
    my($father_outgr)=$outgroup->GetFather();
    $father_outgr->AddChild($outgroup);
}

sub NbFils
{
    my($node)=shift;
    return $node->NbChildren();
}

sub Name
{
    my($node)=shift;
    return $node->Name();
}

###########################################################
####### CHECK FUNCTIONS ###################################
###########################################################

# Do some check on the tree
# Return true (1) if the outgroup need to be removed
sub CheckCorrespondance 
{
    my($tree)=shift;
    my($correspondance)=shift;
    my($name_corres)=shift;
    my($outgroup)=shift;
    my($ret)=0;

#Check if all the leaf are defined in $correspondance
    foreach my $node ($tree->GetNodesList()) {
	my($nb_fils)=$node->NbChildren();
	if ($nb_fils == 0) { # We are on a leaf
	    if (not defined $correspondance->{$node->GetId()}) {
		if (defined($outgroup) && ($node->GetId() eq $outgroup)) {
		    $ret=1;
		} else {
		    #$node->SetCase(0);
		    #$node->SetControl(0);
		    if (not defined($outgroup)) {
			erreur("The leaf '". $node->Name().
			       "' is not in the input file ".
			       "'$name_corres'.\nPerhaps this is the".
			       " outgroup and you need to remove it (option".
			       " --remove-outgroup).\n".
			       "Please, check you data\n", 0);
		    } else {
			erreur("The leaf '". $node->Name().
			       "' is not in the input file ".
			       "'$name_corres'.\nPlease, check you data\n", 0);
		    }
		}
	    }
	}
	
    }

    # Check if all the entries for correspondance whose name begin by
    # H are leafs in the tree

    foreach my $clef (keys %{$correspondance}) {
	if (not $tree->HasNodeIndex($clef)) {
	    erreur("Node '$clef' found in '$name_corres' does".
		   " not exist in the tree. You have probably".
		   " forgot to remove the haplotype corresponding".
		   " to the ancestor.\nPlease, check the input file".
		   " '$name_corres'.\n", 0);
	}
	if ($tree->GetRoot()->Name() eq $clef) {
	    next;
	} else {
	    my($nb_fils)=$tree->GetNode($clef)->NbChildren();
	    
	    if ($nb_fils!=0) {
		erreur("'$clef' present in $name_corres is a".
		       " internal node (not a leaf) in the tree\n".
		       "Please, check your data");
	    }
	}
    }    
    return $ret;
}

#################################################
########  PARCOURS ET REMPLISSAGE ARBRE #########
#################################################

#REWRITE: sub CalculChi2

#REWRITE: sub parcours_nosplit_chi2split

sub FillCaseControl
{
    my($present_node)=shift;
    my($correspondance)=shift;
    if ($present_node->NbChildren()==0)  {
	my($id);
	$id=$present_node->{"id"};
	if (not defined $present_node->{"case"}) {# car sinon, pb pour H000
	    $present_node->{"case"} = $correspondance->{$id}->{"case"};
	}
	if (not defined $present_node->{"control"}) {
	   $present_node->{"control"} = $correspondance->{$id}->{"control"};
	}
	#print $present_node->{"id"}, " m:", $present_node->{"case"}, " c:", $present_node->{"control"}, " ";
    } else {
	my($child);
	# print $present_node->{"id"}, " "; #123456
	for $child ($present_node->GetChildrenList()) { 
	    FillCaseControl($child, $correspondance);
	    $present_node->{"case"}+=$child->{"case"};
	    $present_node->{"control"}+=$child->{"control"};
	}
	#print $present_node->{"id"}, " m:", $present_node->{"case"}," c:", $present_node->{"control"}," ";
    }
    
}

######## QUANTITATIF ########
sub FillQuanti
{
    my($present_node)=shift;
    my($correspondance)=shift;
    if ($present_node->NbChildren()==0)  {
	my($id);
	$id=$present_node->{"id"};
	#if (not defined $present_node->{"quanti"}) {
            # car sinon, pb pour H000
	    #$present_node->{"quanti"}=[];
	#print "cor: ", Dumper($correspondance->{$id}), "\n";
	foreach my $corresp (@{$correspondance->{$id}}) {
	    $present_node->AddQuanti($corresp->[0], $corresp->[1]);
	}
	#}
	
#DEBUG	print STDERR $present_node->{"id"}, " " ;
#DEBUG	for (my $i=0; $i< scalar (@{$present_node->{"quanti"}}); $i++) {
#DEBUG	    print STDERR $present_node->{"quanti"}->[$i]->[0], " (" ,$present_node->{"quanti"}->[$i]->[1], ") ";
#DEBUG	}
#DEBUG	print STDERR "\n";
    } else {
	my($child);
	for $child ($present_node->GetChildrenList()) { 
	    FillQuanti($child, $correspondance);
	    foreach my $quanti (@{$child->GetQuantiList()}) {
		$present_node->AddQuanti($quanti->[0], $quanti->[1]);
	    }
	}
#DEBUG	print STDERR $present_node->{"id"}, " " ;
#DEBUG	for (my $i=0; $i< scalar (@{$present_node->{"quanti"}}); $i++) {
#DEBUG	    print STDERR $present_node->{"quanti"}->[$i]->[0], " (" ,$present_node->{"quanti"}->[$i]->[1], ") ";
#DEBUG	}
#DEBUG	print STDERR "\n";
    }
    
}

#REWRITE: sub ParcoursQuanti

#REWRITE: sub CalculAnovaOneWay

##################################################
######## AFFICHAGE DE L' ARBRE  ##################
##################################################
sub LongueurTrait
{
    my($node)=shift;
    my($level)=shift;
    return("      "x$level."------");
}

sub AffichageParLevel # Ne prend pas une fonction pour l'affichage
{
    my ($racine)=shift;
    my ($len)=6;
    my $AffichageInterne;

    $AffichageInterne= sub
    {
	my ($node) = shift;
	my ($start) = shift; # début commun à tout ce noeud (et descendants)
	my ($up) = shift; # quand on est au dessus de ce noeud
	my ($here) =shift; # quand on affiche ce noeud
	my ($down) =shift; # quand on est au dessous de ce noeud
	
	my($nb_fils)=NbFils($node);
	my($i);
	
	if ($nb_fils >= 1) {
	    $AffichageInterne->($node->GetChild(0),
				$start.$up.(" "x$len), " ", "/", "|");
	}
	for ($i=1; $i<$nb_fils/2; $i++) {
	    $AffichageInterne->($node->GetChild($i),
				$start.$up.(" "x$len), "|", "|", "|");
	}
	print $start.$here.("-"x$len)."* ", Name($node), "\n";
	for ( ;$i < $nb_fils-1; $i++) {
	    $AffichageInterne->($node->GetChild($i),
				$start.$down.(" "x$len), "|", "|", " ");
	}
	if ($nb_fils > 1) {
	    $AffichageInterne->($node->GetChild($nb_fils-1),
				$start.$down.(" "x$len), "|", "\\", " ");
	}
    };
    $AffichageInterne->($racine, "", " ", "-", " ");
}

sub AffichageArbre # Prend une fonction pour l'affichage
{
    my ($racine)=shift;
    my ($function)=shift;
    my ($len)=4;
    my $AffichageInterne;
   
    
    $AffichageInterne= sub
    {
	my ($node) = shift;
	my ($start) = shift; # début commun à tout ce noeud (et descendants)
	my ($up) = shift; # quand on est au dessus de ce noeud
	my ($here) =shift; # quand on affiche ce noeud
	my ($down) =shift; # quand on est au dessous de ce noeud
	my ($at) =shift; # quand on les autres lignes de ce noeud
	
	my($nb_fils)=NbFils($node);
	my($i, $j, $sep);
	my (@tableau)=split (/\n/, $function->($node));
	if ($nb_fils >= 1) {
	    $AffichageInterne->($node->{"children"}->[0],
				$start.$up.(" "x$len), " ", "/", "|", "|");
	}
	for ($i=1; $i<$nb_fils/2; $i++) {
	    $AffichageInterne->($node->{"children"}->[$i],
				$start.$up.(" "x$len), "|", "|", "|", "|");
	}
	print $start.$here.("-"x$len)."* ", $tableau[0], "\n";
	if ($nb_fils > 1) {
	    $sep="|";
	} else {
	    $sep=" ";
	}
	for ($j=1; $j<=$#tableau; $j++) {
	    print $start.$at.(" "x$len).$sep." ", $tableau[$j], "\n";
	}
	#    print $start.$here.("-"x$len)."* ", $function->($node), "\n";
	for ( ;$i < $nb_fils-1; $i++) {
	    $AffichageInterne->($node->{"children"}->[$i],
				$start.$down.(" "x$len), "|", "|", "|", "|");
	}
	if ($nb_fils > 1) {
	    $AffichageInterne->($node->{"children"}->[$nb_fils-1],
				$start.$down.(" "x$len), "|", "\\", " ", " ");
	}
    };
    $AffichageInterne->($racine, "", " ", "-", " ", " ");
}

###############################################################
## FONCTION DEFINISSANT LES INFOS QUI VONT ETRE AFFICHEES #####
###############################################################

# Return results of the test: ddl, p_value, significatif or not, texte and warning
sub TestInfos
{
    my($node)=shift;
    return InfosAffichees($node, 2);

#    my($chaine)=Name($node)."\n";
#    my($lbl_test)=0;
#    my $test;
#    if (defined $node->{"res_test"}) {
#       	for $test (@{$node->{"res_test"}}) {
#	    $chaine.="[".$test->{"level"}."]"." ddl= ".$test->{"ddl"}.
#		" chi2= ".$test->{"chi2"}.
#		" p_value_chi2= ".$test->{"p_val"}.
#	}
#    }
#    return($chaine); 
    
}
sub AssociationInfos
{
    my($node)=shift;
    return InfosAffichees($node, 1);
}

sub TreeInfos
{
 my($node)=shift;
    return InfosAffichees($node, 0);
}

sub InfosQuanti
{
 my($node)=shift;
    return InfosAffichees($node, 3);
}

sub InfosQuantiNoperm
{
 my($node)=shift;
    return InfosAffichees($node, 4);
}
#Return ddl, level, pvalues and chi2

#REWRITE: sub InfosAffichees

##########################################################
######## MODIFICATIONS/CALCULS SUR L'ARBRE ###############
##########################################################

sub FusionBrNulles
{
    my($present_node)=shift;
    my($child);
    my($nb_fils)=NbFils($present_node);
    
    $present_node->{"label"}=$present_node->{"id"};
    $present_node->RecordOrigFather();
    if ($nb_fils != 0) { # on n'est pas dans une feuille
	$present_node->AddOldChild($present_node->GetChildrenList());
	$present_node->ForgetChildren();
	foreach $child ($present_node->GetOldChildrenList()) {
	    if (! FusionBrNulles($child)) {
		$present_node->AddChild($child);
	    }
	}
	if (not $present_node->HasFather()) {
	    return 0;
	}
	if (not defined $present_node->{"br_len"}) {
	    print STDERR "Branch lenght not defined for ", $present_node->{"id"}, "\n"; 
	    exit 1;
	} elsif ($present_node->GetBrLen() == 0) { # branche nulle
	    #print "brnulle ", $present_node->{"id"}, " ";
	    foreach $child  (@{$present_node->{"children"}}) {
		$child->{"father"}=$present_node->{"father"}; #remplace father
		#print "father's name ", $present_node->GetFather()->Name(),"\n"; 
		$present_node->GetFather()->AddChild($child);
	    }
	    $present_node->{"father"}->{"label"}.="+(".$present_node->{"label"}.")";
	    return 1;
	}
    }
    return 0;
}


##########################################################
################### CLEAN FUNCTION  ######################
##########################################################

sub CleanStats
{
    my($tree)=shift;
    
    foreach my $node ($tree->GetNodesList()) {
	$node->EraseCase();
	$node->EraseControl();
	$node->EraseQuanti();	
	delete $node->{"res_test"};
    }
}

##########################################################
########### FUNCTIONS FOR ASSOCIATION TEST ###############
##########################################################

############ QUALITATIF ###############

# From the hash correspondance, fill the variables necessary for Resampling
sub Correspond2Resampling
{
    my($correspondance)=shift;
    my($haploID, $ref_effectif, $total_mal, $total_tem);
    foreach $haploID (keys %{$correspondance}) {
	$ref_effectif->{$haploID}=$correspondance->{$haploID}->{"case"}+
	    $correspondance->{$haploID}->{"control"};
	$total_mal+=$correspondance->{$haploID}->{"case"};
	$total_tem+=$correspondance->{$haploID}->{"control"};
    }
   #  DEBUG print "total_mal=$total_mal total_tem=$total_tem\n";
    return ($total_mal, $total_tem, $ref_effectif);
}


sub Resampling # repompée intégralement de tree_resampling puis modifiée....
{
    my($total_mal) = shift;
    my($total_tem) = shift;
    my($ref_effectif) = shift; # ref on a hash: keys=H002 value= nbmal+nb_tem
    my($clefs, $alea, $i);
    my($new_correspondance);
    foreach $clefs (keys %{$ref_effectif}) { 
	$new_correspondance->{$clefs}->{"case"}=0;
	$new_correspondance->{$clefs}->{"control"}=0;
	for ($i=0; $i<$ref_effectif->{$clefs}; $i++) {
	    $alea=rand($total_mal+$total_tem);
	    # print "alea=$alea";
	    if ($alea < $total_mal) {
		#	print "inf\n";
		$total_mal--;
		$new_correspondance->{$clefs}->{"case"}++;
		$new_correspondance->{$clefs}->{"control"}+=0;
	    } else {
		#	    print "sup\n";
		$total_tem--;
		$new_correspondance->{$clefs}->{"control"}++;
		$new_correspondance->{$clefs}->{"case"}+=0;
	    }
	}
#	print "clefs:$clefs nb_mal=$nb_mal{$clefs}, nb_tem=$nb_tem{$clefs}\n";
    }
    return ($new_correspondance);
}


sub computeTreeStructure($) {
    my $racine=shift;
    my $nextIndex=shift;

    my (@leaf_refs, @leaf_depth, @leaf_parent, @nleaf_parent);

    my $nleaf2id={}; # Hash $node -> [index, $node->Father()]

    my $computeNLeaf;
    $computeNLeaf = sub {
	my $nodeList=shift;
	my $max_depth=0;
	my $nextIndex=0;
	my @childNodeList=();
	my $parent;

	foreach my $node (@{$nodeList}) {
	    push @childNodeList, $node->GetChildrenList();
	}
	if (scalar(@childNodeList) > 0) {
	    ($max_depth,$nextIndex)=$computeNLeaf->(\@childNodeList, $nextIndex);
	}
	foreach my $node (@{$nodeList}) {
	    # Si on est un noeud interne sans id, on définit un index
	    if (NbFils($node) != 0 && !exists($nleaf2id->{$node})) {
		$nleaf2id->{$node}=[$nextIndex++, $node->Father()];
	    }
	}
	if (scalar(@childNodeList) == 1) {
	    # Dans le cas où on a juste une seule branche, on "fusionne"
	    # avec le noeud père (on lui assigne à ce dernier le même index
	    # que nous)
	    my $parent=$childNodeList[0]->Father();
	    if (defined($parent)) {
		$nleaf2id->{$parent}=[$nleaf2id->{$childNodeList[0]}->[0], $parent->Father()] ;
	    }
	}
	if (scalar(@childNodeList) > 0) {
	    $max_depth++;
	}
	return ($max_depth, $nextIndex);
    };
    my ($max_depth,$nbIntNodes)=$computeNLeaf->([$racine], 0);
    foreach my $node (keys(%{$nleaf2id})) {
	my ($id_my, $parent)=@{$nleaf2id->{$node}};
	if (!defined($parent) || !exists($nleaf2id->{$parent})) {
	    $nleaf_parent[$id_my]=-1;
	} else {
	    my $id_parent=$nleaf2id->{$parent}->[0];
	    if ($id_parent != $id_my) {
		$nleaf_parent[$id_my]=$id_parent;
	    }
	}
    }

    my $computeLeaf;
    $computeLeaf = sub {
	my $nodeList=shift;
	my $depth=shift;

	my @childNodeList=();
	my @leafList=();

	foreach my $node (@{$nodeList}) {
	    if (NbFils($node) != 0) {
		push @childNodeList, $node->GetChildrenList();
	    } else {
		push @leafList, $node;
	    }
	}
	if (scalar(@childNodeList) > 0) {
	    my $rec_depth=$depth;
	    if (scalar(@childNodeList) > 1) {
		$rec_depth++;
	    }
	    $computeLeaf->(\@childNodeList, $rec_depth);
	}
	foreach my $node (@leafList) {
	    push @leaf_refs, $node;
	    push @leaf_depth, $depth;
	    push @leaf_parent, $nleaf2id->{$node->Father()}->[0];
	}
	return;
    };
    $computeLeaf->([$racine], 0);

    return (\@leaf_refs, \@leaf_depth, \@leaf_parent, \@nleaf_parent, $max_depth);
}

sub computeChi2 {
    my $leaf_refs=shift;
    my $leaf_depth=shift;
    my $leaf_parent=shift;
    my $nleaf_parent=shift;
    my $max_depth=shift;
    my $correspondance=shift;
    my $prolonge=shift;
    my $sign_util=shift;

    my @results;
    my $nb_leaves=scalar(@{$leaf_refs});
    my $first_leaf=0;
    my $nb_used_leaves=0;
    my @intNodes;
    my %next_intNodes_id;
    for(my $depth = $max_depth; $depth > 0; $depth--) {
	my @leaves;
	my @prev_intNodes_id=keys %next_intNodes_id;
	%next_intNodes_id=();
	my $next_first_leaf;
	{
	    my $i;
	    for ($i=$first_leaf; $i<$nb_leaves && $leaf_depth->[$i] == $depth; $i++) {
		my $id=$leaf_refs->[$i]->GetId();
		push @leaves, $correspondance->{$id};
		my $parent=$leaf_parent->[$i];
		if ($parent != -1) {
		    $next_intNodes_id{$parent} = 1;
		    $intNodes[$parent]->{"case"} += $correspondance->{$id}->{"case"};
		    $intNodes[$parent]->{"control"} += $correspondance->{$id}->{"control"};
		}
	    }
	    $next_first_leaf=$i;
	}
	if ($prolonge == 1) {
	    for (my $i=$next_first_leaf; $i<$nb_leaves; $i++) {
		my $id=$leaf_refs->[$i]->GetId();
		push @leaves, $correspondance->{$id};
	    }
	}
	foreach my $id (@prev_intNodes_id) {
	    push @leaves, $intNodes[$id];
	    my $parent=$nleaf_parent->[$id];
	    if ($parent != -1) {
		$next_intNodes_id{$parent} = 1;
		$intNodes[$parent]->{"case"} += $intNodes[$id]->{"case"};
		$intNodes[$parent]->{"control"} += $intNodes[$id]->{"control"};
	    }
	}

	my $nb_used_leaves=$next_first_leaf - $first_leaf;
	$first_leaf=$next_first_leaf;

	my($test_results)={};
	CalculChi2(\@leaves, scalar(@leaves)-1,
		   $test_results, $sign_util );
	unshift @results, $test_results->{"chi2"};
    }
    return \@results;
}

sub RepeatAssociation
{
    my($tree)=shift;
    my($dataqual)=shift;
    my($correspondance)=shift;
    my($prolonge)=shift;
    my($nb_permutation)= shift;
    my($sign_util)=SignUtil::NO;

    my($racine)=$tree->GetRoot();

    my($ligne_stats)=[];
    print "\n Number of permutation: $nb_permutation\n";
    my($value_per_line, $test_res);
    my($corrected_values);

    if ($dataqual == DataQual::QUALI) {
	my($leaf_refs,$leaf_depth,$leaf_parent,$nleaf_parent, $max_depth)
	    =computeTreeStructure($racine);
	$value_per_line=$max_depth;

      if(0) {
	my $res=computeChi2($leaf_refs,$leaf_depth,$leaf_parent,$nleaf_parent,
			    $max_depth, $correspondance,
			    $prolonge, $sign_util);
	
	push @{$ligne_stats}, @{$res};
	$value_per_line=scalar(@{$res});

	my($total_mal, $total_tem, $effectif);
	($total_mal, $total_tem, $effectif)=Correspond2Resampling($correspondance);

	for (my $i=0; $i<$nb_permutation; $i++) {
	    my $new_correspondance=Resampling($total_mal, $total_tem, $effectif);
	    my $res=computeChi2($leaf_refs,$leaf_depth,$leaf_parent,$nleaf_parent, $max_depth, $new_correspondance,
				$prolonge, $sign_util);
	    push @{$ligne_stats}, @{$res};
	}
     }
	my $resamp=ALTree::CUtils::ResamplingChi2(
	    $leaf_refs, $leaf_depth, $leaf_parent, $nleaf_parent,
	    $max_depth, $prolonge, $nb_permutation, 0
	    );
	$ligne_stats=$resamp->{"chi2s"};

	$corrected_values=ALTree::CUtils::DoublePermutation
	    ($nb_permutation+1, $value_per_line, $ligne_stats);

    } else {
	foreach $test_res (@{$racine->{"res_test"}}) {
	    # Si on n'a qu'une seule branche, la p-value n'est pas définie
	    if ($test_res->{"nb_facteurs"} > 1) {
		$value_per_line++;
	    }
	}
	StockeQuanti($ligne_stats, $racine); # F values corresponding to the real data are put to @{$ligne_F}
	#print STDERR Dumper($correspondance), "\n";
	for (my $i=0; $i<$nb_permutation; $i++) {
	    CleanStats($tree);
	    AssociationQuanti($racine, $correspondance, $prolonge, $sign_util);
	    StockeQuanti($ligne_stats, $racine);
	}
	
	$corrected_values=ALTree::CUtils::DoublePermutation
	    ($nb_permutation+1, $value_per_line, $ligne_stats);
    }



    return($value_per_line, $corrected_values);
}

############ QUANTITATIF ###############
# From the hash correspondance, fill the variables necessary for Resampling
sub Correspond2ResamplingQuanti
{
    my($correspondance)=shift;
    my($haploID, $nbval_per_haplo, @valeurs_tot);
    foreach $haploID (keys %{$correspondance}) {
	$nbval_per_haplo->{$haploID}=scalar(@{$correspondance->{$haploID}});
	foreach my $valeurs (@{$correspondance->{$haploID}}) {
	    push (@valeurs_tot, $valeurs->[0]);
	}
    }
#DEBUG  print STDERR" TOUTES VAL:\n";
#DEBUG  foreach my $val (@valeurs_tot) {
#DEBUG	print STDERR $val, " " ;
#DEBUG	}
#DEBUG  print STDERR"\n";
#DEBUG    print STDERR " NB_VAL_PER_HAPLO ";
#DEBUG    foreach $haploID (keys %{$nbval_per_haplo}) {
#DEBUG	print STDERR $haploID, " " ,$nbval_per_haplo->{$haploID}, "\n";
 #DEBUG   } 
    return (\@valeurs_tot, $nbval_per_haplo);
}

sub ResamplingQuanti
{
    my $valeurs_tot = shift;
    my $nbval_per_haplo = shift;
    my $new_correspondance;
    my $num_haplo=0;
    foreach my $haploID (keys %{$nbval_per_haplo}) {
	$num_haplo++;
	for (my $i=0; $i<$nbval_per_haplo->{$haploID}; $i++) {
	    my $nb= scalar (@{$valeurs_tot});
	    my $alea = int(rand($nb)); # Je récupère bien un nb entre 0 et le dernier élément du tab
	#DEBUG    print STDERR "ALEA$haploID ", $alea , " ";
	    $new_correspondance->{$haploID}->[$i]->[0]=$valeurs_tot->[$alea];
	    $new_correspondance->{$haploID}->[$i]->[1]=$num_haplo;
	    splice(@{$valeurs_tot}, $alea, 1);
	}
    }
#DEBUG    print "\n TEST\n";
    
#DEBUG    foreach my $haploID (keys(%{$new_correspondance})) {
#DEBUG	print STDERR $haploID , " " ;
#DEBUG	for (my $i=0; $i<scalar(@{$new_correspondance->{$haploID}}); $i++) {
#DEBUG	    print STDERR $new_correspondance->{$haploID}->[$i]->[0], " "; 
#DEBUG	}
#DEBUG	print STDERR "\n";
#DEBUG    }
    return ($new_correspondance);

}
sub StockeQuanti

{
    my($ligne_chi2)=shift;
    my $racine=shift;
    my $test_res;
    foreach $test_res (@{$racine->{"res_test"}}) {
	# Si on n'a qu'une seule branche, la p-value n'est pas définie
	if ($test_res->{"nb_facteurs"} > 1) {
	    push @{$ligne_chi2}, $test_res->{"F"};
	}
    }
    # Fill a 2*n table: each line containig chi2 and each columns
    # corresponding to one repetition
    #push (@{$table_of_line}, \@ligne_chi2);
}


sub AssociationQuanti 
{
    my($racine)=shift;
    my($correspondance)=shift;
    my($prolonge)=shift;
    my($sign_util)=shift;
    my($valeur_tot, $nbval_per_haplo, $new_correspondance);
    ($valeur_tot, $nbval_per_haplo)=Correspond2ResamplingQuanti($correspondance);
    ($new_correspondance)=ResamplingQuanti($valeur_tot, $nbval_per_haplo);

    # DEBUG  my($haploID);
    #foreach $haploID (keys %{$new_correspondance}) {
    #	print "Haplo: $haploID mal= ",$new_correspondance->{$haploID}->{"case"}, " tem=", $new_correspondance->{$haploID}->{"control"},"\n";
    #}  
    #print "\n";
    
    FillQuanti($racine,$new_correspondance);
    ParcoursQuanti($racine->{"children"}, $prolonge, 
			       SplitMode::NOSPLIT, $racine, $sign_util);
} 

##########################################################
################# LOCALISATION ###########################
##########################################################

sub CalculateRit 
{
    my($tree)=shift;
    my($s_site_nb)=shift;
    my($s_state)=shift;
    my($co_evo)=shift;
    my($clef);
    my($info_mutation); 
    my($s_t, $s_t_rev)=(0,0);
    
    my($s_sitesens_per_tree)=$tree->GetSite($s_site_nb)
	->GetSens($s_state);

    if (not defined($s_sitesens_per_tree)) {
	warn("No S site (number $s_site_nb)".
		" with sens '".$s_state->GetLabel()."' found\n");
    }
    # Calcul du Rit
    # Notation des commentaires: le s_site mute de T->M
    foreach my $node ($s_sitesens_per_tree->GetNodesList()) {
	foreach my $sitesens_per_tree ($node->GetApoList()) {
	    if ($sitesens_per_tree == $s_sitesens_per_tree) {
		# on profite du fait pour incrémenter le nombre de
		# mutation de s_site
		# (On aurait aussi pu globalement incrémenter s_t
		#  de la taille du tableau du premier foreach :
		#  en effet, on passe ici une fois pour chaque
		#  foreach du second niveau)
		$s_t++;
	    } else {
		# le sitesens co-mute par ex de 1->2
		$sitesens_per_tree->IncRit();
	    }
	}
    }

    if ($co_evo == CoEvo::DOUBLE){
	# On 'retourne' le sens s_state
	$s_state->Switch();
	my $s_sitesens_per_tree_rev=$tree->GetSite($s_site_nb)
	    ->GetSens($s_state);
	# On le remet dans le bon sens après l'avoir tout tourneboulé
	$s_state->Switch();
	
	if (defined($s_sitesens_per_tree_rev)) {
	    # Calcul du Rit
	    foreach my $node ($s_sitesens_per_tree_rev->GetNodesList()) {
		foreach my $sitesens_per_tree ($node->GetApoList()) {
		    if ($sitesens_per_tree == $s_sitesens_per_tree_rev) {
			# nombre de mutation inverse du s_site (M->T)
			# (même remarque que précédemment)
			$s_t_rev++;
		    } else {
			# le sitesens co-mute de 2->1 avec s_site M->T
			# donc on incrémente Rit pour sitesens 1->2
			$sitesens_per_tree->GetSensRev()->IncRit();
		    }
		}
	    }
	}
    } elsif ($co_evo == CoEvo::SIMPLE){
    } else {
	die "Invalid value for co_evo: $co_evo - should be 0 or 1\n";
    }
    return($s_t, $s_t_rev);
}

sub CalculateEit
{ 
    my($tree)=shift; # Comme d'hab
    my($s_site_nb)=shift; # Pas utile ici !
    my($s_t)=shift; # nombre de fois où S mute T->M
    my($s_t_rev)=shift; # nombre de fois où S mute M->T
    my($b_t)=shift; # nombre de branches au total dans l'arbre 
                    # (après fusion des branches nulles)

    # Calcul du Eit (et implicitement le Mit)
    foreach my $site_per_tree ($tree->GetSitesList()) {
	foreach my $sitesens_per_tree ($site_per_tree->GetSensList()) {
	    $sitesens_per_tree->
		SetEit(($sitesens_per_tree->GetMit()*$s_t+
			$sitesens_per_tree->GetSensRev()->GetMit()*$s_t_rev)
		       /$b_t);
	  #  print "  m_it= ", $sitesens_per_tree->GetMit(),"\n";
	  #  print $site_per_tree->GetSiteNb(),"  E_i ", $sitesens_per_tree->GetEit(), "\n";
	  # print  " R_it= ", $sitesens_per_tree->GetRit(), "\n"; 
	  #  print " V_it= ", $sitesens_per_tree->GetVit(),"\n";
	}
    }
}

sub PrintAllVit
{
    my($tree)=shift;
    my($s_site_nb)=shift;
    my($mutation , $sens); 
    foreach my $site_per_tree ($tree->GetSitesList()) {
	if ($site_per_tree->GetSiteNb() == $s_site_nb) {
	    next;
	}
	foreach my $sitesens_per_tree ($site_per_tree->GetSensList()) {
	    print "mutation= ",$site_per_tree->GetSiteNb(), "\t";
	    print "sens= ", $sitesens_per_tree->GetSensLabel(),"\n";
	    print "  m_it= ", $sitesens_per_tree->GetMit(),
	    " R_it= ", $sitesens_per_tree->GetRit(), 
	    " E_it= ", $sitesens_per_tree->GetEit(),
	    " V_it= ", $sitesens_per_tree->GetVit(),"\n";
	}
    }

}

sub PrintAllVi
{
    my($foret)=shift;
    my($s_site_nb)=shift;
    my($mutation , $sens); 
    foreach my $site_per_foret ($foret->GetSitesList()) {
	if ($site_per_foret->GetSiteNb() == $s_site_nb) {
	    next;
	}
	foreach my $sitesens_per_foret ($site_per_foret->GetSensList()) {
	    print "mutation= ",$site_per_foret->GetSiteNb(), "\t";
	    print "sens= ", $sitesens_per_foret->GetSensLabel(),"\t";
	    print "  V_i= ", $sitesens_per_foret->GetVi(),"\n";
	}
    }

}

# Pour chaque site, on a choisi un sens en fonction du Vi (max) et on
# affiche le tableau des Vi pour tous les sites et pour le sens choisi
sub PrintViMax
{
    my($foret)=shift;
    my($s_site_nb)=shift;
    foreach my $site ($foret->GetViMaxSiteList()) {
	if ($site->GetSiteNb() == $s_site_nb) {
	    next;
	}
	print "site number ", $site->GetSiteNb(), "\n";
	foreach my $sens ($site->GetViMaxSensList()) {
	    print "\tsens ", $sens->GetSensLabel(), "\t";
	    print "V_i = ", $sens->GetVi(), "\n";
	}
    }
}

# Affiche le tableau de tous les Vi, pour tes les sites et pour tous
# les sens (triés par ordre décroissant)
sub PrintViMaxSens
{
    my($foret)=shift;
    my($s_site_nb)=shift;
    foreach my $sens ($foret->GetViMaxSensList()) {
	if ($sens->GetSiteNb() == $s_site_nb) {
	    next;
	}
	print "site number ", $sens->GetSiteNb(), "\t";
	print "\tsens ", $sens->GetSensLabel(), "\t";
	print "V_i = ", $sens->GetVi(), "\n";
    }
}

###########################################################
# Fonctions du prog principal #############################
###########################################################
sub check_tree_numbers
{
    my $num_tree_in_file=shift;
    my $user_tree_numbers=shift;
    my @tree_numbers;
    my %selected_trees;

    for my $tree_num (@{$user_tree_numbers}) {
	if ($tree_num < 1) {
	    erreur("Invalid tree-to-analyse $tree_num\n", 0);
	  } elsif ($tree_num > $num_tree_in_file) {
	      erreur("Invalid tree-to-analyse $tree_num (only $num_tree_in_file in file)\n", 0);
	  }
	if (defined($selected_trees{$tree_num})) {
	    erreur("Invalid tree-to-analyse $tree_num (already selected)\n", 0);
	}
	$selected_trees{$tree_num}=1;
	# Correspondance entre les numéros de l'utilisateur et les
	# indices dans le tableau (ie commence à 1 ou à 0)
	push @tree_numbers, $tree_num-1;
    }
    return \@tree_numbers;
}

sub select_trees
{
    my $max=shift;
    my $nb=shift;

    my @tab;
    for (my $i=0; $i<$max; $i++) {
	$tab[$i]=$i;
    }
    my @selected;
    for (my $i=0; $i<$nb; $i++) {
	my $alea=int(rand($max--));
	push @selected, $tab[$alea];
	splice(@tab, $alea, 1);
    }
    return \@selected;
}

sub PrintTree { 
    my $tree=shift;

    my ($racine)=$tree->GetRoot();
    AffichageArbre($racine, \&TreeInfos);
    print "\n\n";
}

sub SwitchRoot {
    my $tree= shift;
    my $outgroup=shift;

    my $root=$outgroup->GetFather();
    $tree->ChangeRoot($root);    
    my $newracine=$tree->GetRoot();
}

################################################################
#################### TRAITE PLS FICHIERS #######################
################################################################

sub TraiteSeveralFiles
{
#DEBUG    print STDERR "Pool rés des fichiers\n";
    my $res_files = shift;
    my ($nb_files);
    
    my $p_finale=0;
    my $pmin_min=99;
    my @distrib=();
    for (my $i=0; $i<=$#$res_files; $i++) {
	shift(@{$res_files->[$i]->{"distrib_pmin"}});
	push (@distrib, @{$res_files->[$i]->{"distrib_pmin"}});
	if ( $res_files->[$i]->{"pmin"} < $pmin_min) {
	    $pmin_min = $res_files->[$i]->{"pmin"};
	}
    }
    print "La pmin minimale est $pmin_min, la taille de la distrib pmin est ", scalar(@distrib), "\n";
    
   # sort {$a <=> $b}  @distrib; # trie du plus petit au plus grand
    my $compte=0;
    foreach my $elem (@distrib) {
	$compte ++;
#DEBUG	print STDERR  "elem=$elem\n";
	if ($elem<=$pmin_min) {
	    $p_finale ++;
	}
    }
#	} else {
    return $p_finale/$compte;
#	}
}


###########################################################

sub manage_options
{
    my %options;
    my $result;
    
    my $choix={ "data-type" => 
		{ "snp" => DataType::SNP,
		  "dna" => DataType::DNA,
	        },
		 "data-qual" => 
		{ "quantitative" => DataQual::QUANTI,
		  "qualitative" => DataQual::QUALI,
	        },   
		"rootmeth" => 
		{ "outgroup" => RootMeth::OUTG,
		  "ancestor" => RootMeth::ANC,
	        },  
		"tree-building-program" =>
		{ "phylip" => PhylProg::PHYLIP,
		  "paup" => PhylProg::PAUP,
		  "paml" => PhylProg::PAML,
		  },
		"splitmode" =>
		{ "chi2split" => SplitMode::CHI2SPLIT,
		  "nosplit" => SplitMode::NOSPLIT,
		  },
		"co-evo" =>
		{ "simple" => CoEvo::SIMPLE,
		  "double" => CoEvo::DOUBLE,
		  },
		};

    my $handle_choix = sub {
	my $option=shift;
	my $value=shift;

	foreach my $key (keys %{$choix->{$option}}) {
	    if ($key=~/^$value/i) {
		$options{$option."-value"}=$choix->{$option}->{$key};
		return;
	    }
	}
	die "Option '$option': unauthorized value '$value'\n";
    };
    my $handle_args = sub {
	my $name=shift;
	die "What about '$name' ?\n";
    };
    my $handle_progname = sub {
	my $name=shift;
	die "What about '$name' ?\n";
    };
	

    %options=("<>" => $handle_args);

    foreach my $option (keys %{$choix}) {
	$options{$option}=$handle_choix;
    }

    GetOptions (\%options,
		"version",
		"short-help|h",
		"help",
		"man",
		"association|a!", # !!! demander pour le ! à Vince
		"s-localisation|l!",	
		"first-input-file|i=s",
                "second-input-file|j=s",
                "output-file|o=s",
                "data-type|t=s",
		"data-qual|q=s",
		"remove-outgroup!",
		"outgroup=s",
                "tree-building-program|p=s",
                "splitmode|s=s",
		"no-prolongation!",
		"chi2-threshold|n=f",
		"permutations|r=i",
		"number-of-trees-to-analyse=i",
		"tree-to-analyse=i@",
		"s-site-number=i",
		"s-site-characters=s",
		"co-evo|e=s",
		"print-tree!",
		"anc-seq=s",
		"nb-files=i",
		"<>",
		) or pod2usage(2);
    if (defined($options{"version"})) {
	print $0, " version ", $VERSION, "\n";
	print "(Perl version ", $], ")\n";
	exit 0;
    }
    if (defined($options{"short-help"})) {
	pod2usage(-exitstatus => 0, -verbose => 0);
    }
    if (defined($options{"help"})) {
	pod2usage(-exitstatus => 0, -verbose => 1);
    }
    if (defined($options{"man"})) {
	pod2usage(-exitstatus => 0, -verbose => 2);
    }

    delete($options{"<>"});
    foreach my $option (keys %{$choix}) {
	delete($options{$option});
	if (exists($options{$option."-value"})) {
	    $options{$option}=$options{$option."-value"};
	    delete($options{$option."-value"});
	}
    }
    
    return \%options;
}

sub main
{
    my($rec_program);
    my $result;
    my $options;
    
    my $option_require = sub {
	my $option=shift;
	my $texte=shift;
	if (not exists($options->{$option})) {
	    my $msg="Error: option '$option' needed";
	    if (defined($texte)) {
		$msg.="\n".$texte;
	    }		
	    pod2usage("Error: option '$option' required");
	}
    };
    my $option_value = sub {
	my $option=shift;
	my $default=shift;
	my $msg=shift;
	if (not exists($options->{$option})) {
	    if (defined $msg) {
		print STDERR $msg, "\n";
	    }
	    return $default;
	}
	return $options->{$option};
    };
    my $option_selected = sub {
	my $option=shift;
	my $texte=shift;
	if (not exists($options->{$option})) {
	    my $msg="option '$option' not selected";
	    if (defined($texte)) {
		$msg.="\n".$texte;
	    }		
	    return 0;
	} else {
	    return 1;
	}
    };

    $options=manage_options();
    
    my($prolonge)=$option_value->("prolongation", 1);
    my($print_tree)=$option_value->("print-tree", 0);

    $option_require->("tree-building-program");
    my($phylo_program)=$option_value->("tree-building-program");  
   
    my($association)=$option_value->("association", 0);

    my $nb_files = $option_value->("nb-files", 0);
    if ($nb_files == 0) {
	$nb_files=1
	}
    $option_require->("first-input-file");
    my $input_files=$option_value->("first-input-file");
    my @input_files=split(/:/, $input_files);
    my $nb_read_files = scalar(@input_files);
    if ($nb_read_files != $nb_files) {
	erreur("The number of files really given ($nb_read_files) is not equal to".
	       " the number of files specified in the option --nb-files ($nb_files)\n",0)
	}
    # name of the file containing haploID, nb case and nb control
    my($tot_name_corres)=$option_value->("second-input-file", "correspond.txt");
    my @tab_name_corres = split(/:/, $tot_name_corres);

    my($s_loc)=$option_value->("s-localisation", 0);
    if (!$association && !$s_loc) {
	erreur("Should I perform the association test or the".
	       " localisation test ?\n(use option '--association'".
	       " or '--s-localisation'.)\n", 1);
    }
    if ($option_selected->("output-file")) {
	my $out=$option_value->("output-file");
	open(STDOUT, '>', $out) 
	    or erreur("Unable to write in '$out': $!\n", 0);
    }

    my $datatype;
    my $ancetre="";
    if ($phylo_program == PhylProg::PHYLIP) {
	$option_require->("data-type");
	$datatype=$option_value->("data-type");
	if ($option_selected->("anc-seq")) {
	    $ancetre=$option_value->("anc-seq");
	}
    }
    
    my $dataqual="";
    if ($association == 1) {
	$option_require->("data-qual");
        $dataqual = $option_value->("data-qual");
    }
    ###########################################################
    # Récupération et précalcul des arbres qu'on va utiliser
    ###########################################################

    my @trees;
    my @res_files; #contient les ref sur les résultats de chaque fichier analysé
    my($all_files_foret)=ALTree::Foret->New();
    my $s_site_nb;
    {
	for (my $file=1; $file<=$nb_files; $file ++) {
	    my $input=$input_files[$file-1];
#DEBUG	    print STDERR $input, "\n";
	    print STDERR "Analyzing file number $file\n";
	   my ($input_file)=ALTree::Input::ReadInputFile1($input, 
							$phylo_program, 
							$datatype, $ancetre);
	    
	    my $num_trees_in_file=$input_file->{"nb_trees"};
	    
	    my(@no_tree);
	    my($user_tree_numbers)=$option_value->("tree-to-analyse", \@no_tree);
	    my $tree_numbers=check_tree_numbers($num_trees_in_file, $user_tree_numbers);
	    
	    my $nb_tree_selected=scalar(@{$tree_numbers});
	    if ($association) { 
		if ($nb_tree_selected == 0) {
		    $tree_numbers=select_trees($num_trees_in_file, 1);
		} elsif ($nb_tree_selected > 1) {
		    erreur("Only one tree can be selected for association\n");
		}
	    }
	    if ($s_loc) {
		if ($nb_tree_selected == 0) {
		    my($trees_to_analyse)=
			$option_value->("number-of-trees-to-analyse", undef);
		    if (not defined($trees_to_analyse)) {
			print STDERR "Warning: no option".
			    "number-of-trees-to-analyse or tree-to-analyse:".
			    " using all ($num_trees_in_file)".
			    " trees";
			$trees_to_analyse=$num_trees_in_file;
		    }
		    if ($trees_to_analyse>$num_trees_in_file) {
			erreur("Invalid number of trees to analyse :".
			       " your file contains only".
			       " $num_trees_in_file trees\n", 0);
		    }
		    if ($trees_to_analyse<1) {
			erreur("Not enought trees to analyse".
			       " ($trees_to_analyse)\n", 0);
		    }
		    $tree_numbers=select_trees($num_trees_in_file, $trees_to_analyse);
		} elsif ($nb_tree_selected > $num_trees_in_file) {
		    erreur("Invalid number of trees to analyse :".
			   " your file contains only".
			   " $num_trees_in_file trees\n", 0);
		}
	    }
	    
	    my $anctype=$input_file->{"anctype"};
	    if ($anctype == ANC::Rooted && 
		($option_selected->("outgroup") 
		 || $option_selected->("remove-outgroup"))) {
		erreur("You cannot use the options '--outgroup' or".
		       " '--remove-outgroup' because your input file '"
		       .$input_file->{"filename"}."' contains trees rooted".
		       " with an ancestral sequence.\n", 1);
	    }
	    if ($anctype == ANC::OutGroup) {
		my $outgroup=$input_file->{"outgroup"};
		if (defined($outgroup)) {
		    if ($option_value->("outgroup", $outgroup) ne $outgroup) {
			erreur("The option '--outgroup' tell me to use '".
			       $option_value->("outgroup")."' as an outgroup ".
			       "whereas the file '".$input_file->{"filename"}.
			       "' contains trees rooted with the outgroup '".
			       $outgroup."'.\n", 1);
		    }
		} else {
		    $input_file->{"outgroup"}=$option_value->("outgroup", undef);
		}
	    }
	    my ($switchroot)="";
	    if ($anctype == ANC::Unrooted) {
		$input_file->{"outgroup"}=$option_value->("outgroup", undef);
		$switchroot=$option_value->("outgroup", 0);
	    }
	    if ($option_selected->("remove-outgroup")
		&& not defined($input_file->{"outgroup"})) {
		if ($anctype == ANC::OutGroup) {
		    erreur("You tell me to remove the outgroup for the".
			   " analyses. ".
			   "The input file '".$input_file->{"filename"}.
			   "' contains trees rooted using ".
			   "an outgroup but I cannot automaticaly find which".
			   "one has been used.\nPlease, provide me the".
			   "outgroup (option '--outgroup')\n", 0);
		} else {
		    erreur("You tell me to remove the outgroup for the".
			   " analyses, ".
			   "however, I do not know what is the outgroup\n".
			   "Please, use the '--outgroup' option\n", 0);
		}
	    }
	    if ($association
		&& $anctype == ANC::Unrooted
		&& not defined($input_file->{"outgroup"})) {
		erreur("I need a rooted tree to perform the association".
		       " test. However, the input file '".
		       $input_file->{"filename"}."' contains unrooted trees\n".
		       "Please, provide me an outgroup (option '--outgroup')\n"
		       , 1);
		#erreur("I need a rooted tree to perform the association".
		#	   " test.\nThe input file '".$input_file->{"filename"}.
		#	   "' contains trees rooted using ".
		#	   "an outgroup\nbut I cannot automaticaly find which one".
		#	   " has been used.\nPlease, provide me the outgroup".
		#	   " (option '--outgroup')\n"
		#	   , 1);
	    }
	    
	    my $remove_outgroup=$option_value->("remove-outgroup", 0);
	    my @trees =();
	    for my $num_tree (@{$tree_numbers}) {
		my $file_tree=PrepareTree($phylo_program, $input_file, 
					  $datatype, $ancetre, $num_tree);
		push @trees, $file_tree;
		
		my ($tree)=$file_tree->{"tree"};
		my $outgroup;
		if (defined($input_file->{"outgroup"})) {
		    $outgroup=$tree->GetNode($input_file->{"outgroup"});
		    if (not defined($outgroup)) {
			erreur ("I cannot find the outgroup '".
				$input_file->{"outgroup"}."' in the".
				" tree number ".($file_tree->{"index"}+1).
				" in file '".$input_file->{"filename"}.
				"'. It does not correspond to any node!\n", 0);
		    }
		    if (!$switchroot) {
			if ($outgroup->GetFather() != $tree->GetRoot()) {
			    erreur("The outgroup '".$outgroup->Name().
				   "' is not just under the root ".
				   "for the tree number ".
				   ($file_tree->{"index"}+1).
				   " in file '".$input_file->{"filename"}.
				   "'.\nDo you choose the rigth outgroup ?\n", 0);
			}
		    }
		}
		if ($switchroot) {
		    SwitchRoot($tree, $outgroup);
		}
		if ($remove_outgroup) {
		    RemoveOutgroup($tree, $outgroup);
		}
		
		if ($print_tree) {
		    PrintTree($tree);
		}
		if ($s_loc && defined($file_tree->{"has_ambiguity"})) {
		    erreur("Some apomorphies are ambiguous in the".
			   " tree number ".($file_tree->{"index"}+1).
			   " in file '".$input_file->{"filename"}.
			   "' (I find the character state '?').".
			   " It cannot be used for the localisation test.\n", 0);
		}
	    }
	
	    
	    print STDERR "read done\n";
	    ###########################################################
	    # Let's go. D'abord pour l'association
	    ###########################################################
	    
	    if ($association == 1) { 
		my($splitmode)=SplitMode::NOSPLIT; # nosplit est imposé
		
		# name of the file containing haploID, nb case and nb control
		my($name_corres)= $tab_name_corres[$file-1];
		if ($name_corres eq "") {
		    die "check the number of file containing case/controle. There should ".
			"be as many such files as first input files\n";
		}
		my($correspondance); # ref on a hash containing haploID refferring to
		# a hash containing nb case and nb control 
		# or to a tab containing the different quantitative values
		
		if ($dataqual == DataQual::QUALI) {
		    $correspondance=ALTree::Input::ReadCorrespondQuali($name_corres);
		} elsif ($dataqual == DataQual::QUANTI) {
		    $correspondance=ALTree::Input::ReadCorrespondQuanti($name_corres);
		} else {	
		    erreur("invalid choice for option data-qual, $dataqual. It should be either qualitative or quantitative", 0);
		}
		
		my($file_tree)=$trees[0];
		my($tree)=$file_tree->{"tree"};

		my($outgroup)=$file_tree->{"file"}->{"outgroup"};
		if (CheckCorrespondance($tree, $correspondance, $name_corres,
					$outgroup)) {
		    if (!$option_selected->("remove-outgroup")) {
			print STDERR "Warning: assuming option '--remove-outgroup' as".
			    " the outgroup '".$outgroup.
			    "' is not in the file '".$name_corres."'\n";
			
			RemoveOutgroup($tree, $tree->GetNode($outgroup));
		    }
		} else {
		    if ($option_selected->("remove-outgroup")) {
			erreur("You tell me to remove the outgroup '$outgroup',".
			       " however it is present in the file '".
			       $name_corres.". Please, check your data.'\n", 0);
		    }
		}
		#print "\n";
		FusionBrNulles($tree->GetRoot());
		# Structure changée, on recalcul...
		$tree->FillHeight();
		$tree->FillLevel();
		if ($dataqual == DataQual::QUALI) {
		    FillCaseControl($tree->GetRoot(),$correspondance); 
		} else {
		    FillQuanti($tree->GetRoot(),$correspondance); 
		}
		print "\n";
		
		my $racine=$tree->GetRoot();
		my @children=$racine->GetChildrenList();
		
# EN DEVELOPPEMENT: NANOVA
#		{
#		    my $mat;
#		    print "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT\n";
#		    $mat=ALTree::Nanova::Tree2mat($tree->GetRoot());
#		    ALTree::Nanova::WriteMat($mat);
#		    my $values = [];
#		    my $groups = [];
#		    my $nb_term=0;
#		    ALTree::Nanova::FillTableaux($tree->GetRoot(), $values, $groups, $nb_term);
#		    print "Values:\n";
#		    foreach my $val (@{$values}) {
#			print $val, "\t";
#		    }
#		    print "\n";
#		    print "Groups:\n";
#		    foreach my $group (@{$groups}) {
#			print $group, "\t";
#		    }
#		    print "\n";
#		    print "*********\n";
#		    my $una_results; 
#		    $una_results = NAnova::UnbalancedNestedANOVA($mat, $groups, $values);
#		    my $i=1;
#		    foreach my $res (@{$una_results}) {
#			print "Level $i\tF=",$res->{"F"}, "\tp-value=",$res->{"p_value"}, "\n";
#			$i++;
	#	    }
#		}
#		print "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT\n";

		$option_require->("permutations", 
				  "The number of permutations used to calculate".
				  " exact p-values must be specified or set to 0\n");
		
		my $permutation=$option_value->("permutations");
		my($sign_util);
		if ($permutation==0) {
		    #TODO, passer quanti
		    my($seuil_chi2)=$option_value->("chi2-threshold", 0.01, 
						    "Using default chi2 threshold 0.01");
		    ALTree::Chi2::definition_p_chi2($seuil_chi2, 0.01); # mettre une option 
		    # pour seuil test_prop
		    $sign_util = SignUtil::YES # on a besoin de la significativité
		} elsif ($permutation>0) {
		    $sign_util = SignUtil::NO; # on n'a pas besoin de la sign
	        } else {
		    die "invalid value for the number of permutation: $permutation\n";
		}
		print STDERR "Starting tree analysis\n";
	    
		if ($dataqual == DataQual::QUALI) {
		    parcours_nosplit_chi2split(\@children, $prolonge, $splitmode, 
					       $racine, $sign_util);
		} else {
		    ParcoursQuanti(\@children, $prolonge, $splitmode, 
				   $racine, $sign_util); 
		}
		
		{
		    print STDERR "Starting permutations\n";
		    if ($permutation==0) {
			if ($dataqual == DataQual::QUALI) {
			    AffichageArbre($racine, \&TestInfos);
			} else {
			    AffichageArbre($racine, \&InfosQuantiNoperm);
			}
		    } elsif ($permutation>0) {
			my($value_per_line, $corrected_values);
			if ($dataqual == DataQual::QUALI) {
			    AffichageArbre($racine, \&AssociationInfos);
			} else {
			    AffichageArbre($racine, \&InfosQuanti);
			}
			($value_per_line, $corrected_values)=RepeatAssociation
				($tree, $dataqual, $correspondance, $prolonge,$permutation);
			print "\n";
			print "p_val for each level:\n";
			my($i);
			for ($i=0; $i<$value_per_line; $i++) {
			    print "level ", $i+1, " p-value (non corrected) ",
			    $corrected_values->{"chi2"}->[$i], "\n";
			}
			print "corrected minimal p_value in the tree: ", 
			$corrected_values->{"pmin"}, "\n";# at level TODO\n";
			if ($nb_files>1) {
			    $res_files[$file-1]=$corrected_values;
			    #print STDERR "RRR", $res_files[$file-1]->{"distrib_pmin"}->[O], "\n";
			}
		    } else {
			die "invalid value for the number of permutation: $permutation\n";
		    }
		}
	    } 
	    

	    ###########################################################
	    # Let's go. Et pour la localisation
	    ###########################################################
	    
	    if ($s_loc==1) {
		#$option_require->("splitmode");
		#my($splitmode)=$option_value->("splitmode");
		print "Localisation method using S-character\n";
		$option_require->("s-site-number");
		$s_site_nb=$option_value->("s-site-number");
		$option_require->("s-site-characters");
		my($s_char_state)=$option_value->("s-site-characters");
		$option_require->("co-evo");
		my($co_evo_type)=$option_value->("co-evo");

#	print "co_evo_type=$co_evo_type\n";
		
#DEBUG	print "s_char_state=$s_char_state\n";
		if (not ($s_char_state =~ 
			 m/([A-Za-z0-9]+)\s*[-=_]*[>]\s*([A-Za-z0-9]+)/)) {
		    erreur("Invalid character change for character".
			   " S. It should be Anc -> Der\n", 0);
		}
		my($s_anc)=$1;
		my($s_der)=$2;
		my($s_state)= ALTree::Sens->New($s_anc." --> ".$s_der);
		my($foret)=ALTree::Foret->New();
		
		#$Data::Dumper::Indent = 0;
		for my $file_tree (@trees) {
		    my $tree=$file_tree->{"tree"};
		    
		    my ($b_t)=$tree->GetNbBrNonNulle();
		  #  print "NB branches: $b_t\n";
		    my $site=$tree->GetSite($s_site_nb);
		    if (not defined($site)) {
			erreur("Invalid value ($s_site_nb) for".
			       " --s-site-number\n", 0);
		    }
		    
		    my($s_t, $s_t_rev)=CalculateRit($tree, $s_site_nb, $s_state, 
						    $co_evo_type); 
		    CalculateEit($tree, $s_site_nb, $s_t, $s_t_rev, $b_t);
		    $foret->AddTree($tree);
		    $all_files_foret->AddTree($tree);
		}
		$foret->CalculVi();
		#PrintAllVit($foret, $s_site_nb);# Non trié
		#PrintViMax($foret, $s_site_nb); # Affiche la liste en choisissant
		#pour chaque sit, juste le meilleur sens
		print "\n";
		print "Results:\n";
		PrintViMaxSens($foret, $s_site_nb);

	    }
	}
	
	if ($association == 1) { 
	    if ($nb_files>1) {
		my $p_finale = TraiteSeveralFiles(\@res_files, $nb_files);
		print "\n############################################\n";
		print "p_value for the $nb_files files = $p_finale\n";
	    }
	}
	if ($s_loc == 1) {
	    if ($nb_files>1) {
		$all_files_foret->CalculVi();
		print "\n############################################\n";
		print "Results for all the files\n";
		PrintViMaxSens($all_files_foret, $s_site_nb);
	    }
	}
    }
}


sub PleinInfos {
    my $node=shift;
    
    my $suite=""; #"\nligne suivante\net encore après\n ";
    if (defined ($node->{"father"})) {
	return Name($node). "\nFrom: ". Name($node->{"father"}).$suite;
    } else {
	return Name($node).$suite;
    }
}

use FileHandle;
use IPC::Open2;
sub my_test {
    my $pid;
    
    $pid = open2(*Reader, *Writer, "phylip mix" );
    print Writer "r\no\n6\nw\na\n5\ny\nweight\nr\n";
    while (<Reader>) {
	print "PHYL: ", $_;
    }
    
}

#my_test;
#exit 0;

main;

#man perlpod

__END__
    
=head1 NAME

altree - Association and Localisation tests using  phylogenetic Trees

=head1 SYNOPSIS

altree [options]

 Options:
    --version        program version
    --short-help|h   brief help message
    --help           help message with options descriptions
    --man            full documentation
    --association|a  perform the association test
    --s-localisation|l   perform the localisation using the S character
    --first-input-file|i result_file from phylogeny reconstruction programs
    --second-input-file|j file containing the nb of cases/controls carrying an haplotype
    --output-file|o output_file
    --data-type|t DNA|SNP
    --data-qual|q qualitative|quantitative
    --outgroup outgroup_name
    --remove-outgroup 
    --tree-building-program|p phylip|paup|paml
    --splitmode|s nosplit|chi2split
    --no-prolongation
    --chi2-threshold|n value
    --permutations|r number
    --number-of-trees-to-analyse number
    --tree-to-analyse number
    --s-site-number number
    --s-site-characters ancestral state -> derived state
    --co-evo|e simple|double
    --print-tree 
    --anc-seq ancestral sequence (only with phylip)
    --nb-files number of input files to analyse (only for association test)

=head1 OPTIONS

=over 8

=item B<--version>

Print the program version and exits.

=item B<--short-help>

Print a brief help message and exits.

=item B<--help>

Print a help message with options descriptions and exits.

=item B<--man>

Prints the manual page and exits.

=item B<--association|a>

Perform the association test

=item B<--s-localisation|l>

Localise the susceptibility locus using the "S-character method"

=item B<--first-input-file|i> F<result_file>

Input file 1 (paup, phylip or paml results file).
If several input files are analysed, their names must be separated by colons. 
Example: input1:input2 etc

=item B<--second-input-file|j> F<correspond_file>

Input file 2, default F<correspond.txt>.
The number of input file 2 must be the same as the number of input file 1.
The name of the different input file 2 must be separated by colons

=item B<--output-file|o> F<outfile>

Output file

=item B<--data-type|t> C<DNA>|C<SNP>

Type of data: DNA (ATGCU) or SNP (0-1)

=item B<--data-qual|q> C<qualitative>|C<quantitative>

Analyse qualitative (case/control) or quantitative data

=item B<--outgroup>  outgroup

Root the tree with this outgroup

=item B<--remove-outgroup>

Remove the outgroup of the tree before performing the tests

=item B<--tree-building-program|p> C<phylip>|C<paup>|C<paml>

Phylogeny reconstruction program

=item B<--splitmode|s> C<nosplit>|C<chi2split>
    
how tests are performed from a level to another

=item B<--no-prolongation>

No prolongation of branches in the tree

=item B<--chi2-threshold|n> value

Significance threshold for chi2 (default value 0.01)

=item B<--permutations|r> number

Number of permutations used to calculate exact p_values
(Only for association test)

=item B<--number-of-trees-to-analyse> number

Number of trees to analyse in the localisation analysis 
(only for localisation method using S-character)

=item B<--tree-to-analyse> number

With this option, you can specify the tree to use (instead of
random). Can be used several times to specify multiple trees.

=item B<--s-site-number> number

Number of the S character site in the sequence
(only for localisation method using S-character)

=item B<--s-site-characters> transition

Character states for the S character: ancestral state -> derived state
ex: G->C or 0->1  (only for localisation method using S-character)

=item B<--co-evo|e> C<simple>|C<double>

Type of co-evolution indice 
  simple: only the anc -> der transition of S is used 
  double: the two possible transitions are used

=item B<--print-tree>

If this option is selected, the tree will be printed to the output

=item B<--anc-seq> anc_seq

With this option, you can specify the ancestral sequence.
This option is only useful when the tree is reconstructed using the mix program 
of phylip with the ancestral states specified in the file "ancestors"

=item B<--nb-files> number

With this option, you specify the number of input files (1 and 2)  to analyse
This option only works for the association test. Be careful if the number of trees is not
the same for the different input files: if the chosen tree doesn't exist in one file, 
the program wil not work correctly


=back

=head1 DESCRIPTION

B<This program> performs 

(a) an association test between a candidate gene and disease or a quantitative trait

(b) a localsation tests: it allows to detect which SNP is involved in the determinism of the disease or the quantitative trait

These two tests are based on the analysis of haplotype phylogenetic trees.

=cut