diff --git a/scripts/set_processPooledVcf.pl b/scripts/set_processPooledVcf.pl index e01c1fb..131a6b1 100755 --- a/scripts/set_processPooledVcf.pl +++ b/scripts/set_processPooledVcf.pl @@ -56,28 +56,36 @@ sub main{ system("pairwiseDistances.pl --numcpus $$settings{numcpus} < $filteredAlignment | sort -k3,3n | tee $pairwise | pairwiseTo2d.pl > $pairwiseMatrix"); die if $?; - # Process a resulting VCF: make an MSA; remove uninformative sites; find pairwise distances; - # find Fst; make a tree; calculate the eigenvector + my $numSamples=`grep -c ">" $filteredAlignment` + 0; + die "ERROR: could not count the number of samples in $filteredAlignment" if $?; + + # Process a resulting VCF: make an MSA; remove uninformative sites; find pairwise distances; + # find Fst; make a tree; calculate the eigenvector my $indexCase="$$settings{prefix}.eigen.tsv"; system("set_indexCase.pl $pairwise | sort -k2,2n > $indexCase"); die "ERROR with set_indexCase.pl" if $?; # TODO Fst - # RAxML puts everything into the CWD and so we have to work around that. - system("cp $filteredAlignment $$settings{tempdir}/; cd $$settings{tempdir}; launch_raxml.sh -n $$settings{numcpus} $filteredAlignment suffix"); - die "ERROR with launch_raxml.sh" if $?; - for (qw(RAxML_bestTree RAxML_bipartitionsBranchLabels RAxML_bipartitions RAxML_bootstrap RAxML_info)){ - system("mv -v $$settings{tempdir}/$_.suffix $$settings{prefix}.$_"); - die "ERROR: could not move $$settings{tempdir}/$_.suffix: $!" if $?; - } + # Can work with trees if there are enough samples + if($numSamples >= 4){ + # RAxML puts everything into the CWD and so we have to work around that. + system("cp $filteredAlignment $$settings{tempdir}/; cd $$settings{tempdir}; launch_raxml.sh -n $$settings{numcpus} $filteredAlignment suffix"); + die "ERROR with launch_raxml.sh" if $?; + for (qw(RAxML_bestTree RAxML_bipartitionsBranchLabels RAxML_bipartitions RAxML_bootstrap RAxML_info)){ + system("mv -v $$settings{tempdir}/$_.suffix $$settings{prefix}.$_"); + die "ERROR: could not move $$settings{tempdir}/$_.suffix: $!" if $?; + } - # TODO: reroot the tree - #rerootLongestBranch("$$settings{prefix}.RAxML_bipartitions",$settings); + # TODO: reroot the tree + #rerootLongestBranch("$$settings{prefix}.RAxML_bipartitions",$settings); - # Get combined distance statistics on the tree - system("cladeDistancesFromTree.pl -t $$settings{prefix}.RAxML_bipartitions -p $$settings{prefix}.pairwise.tsv --outprefix $$settings{prefix}"); - die "ERROR with cladeDistancesFromTree.pl" if $?; + # Get combined distance statistics on the tree + system("cladeDistancesFromTree.pl -t $$settings{prefix}.RAxML_bipartitions -p $$settings{prefix}.pairwise.tsv --outprefix $$settings{prefix}"); + die "ERROR with cladeDistancesFromTree.pl" if $?; + } else { + logmsg "WARNING: only $numSamples samples are in the alignment; skipping tree-building."; + } return 0; }