Commit 0f96be77 authored by Andreas Tille's avatar Andreas Tille
Browse files

New upstream version 1.6.7.1+dfsg

parent a1731c9c

CMakeLists.txt

+12 −3
Original line number Diff line number Diff line
@@ -51,7 +51,7 @@ add_definitions(-DIQ_TREE) 
# The version number.

set (iqtree_VERSION_MAJOR 1)

set (iqtree_VERSION_MINOR 6)

set (iqtree_VERSION_PATCH "5")

set (iqtree_VERSION_PATCH "7.1")



set(BUILD_SHARED_LIBS OFF)
@@ -513,12 +513,21 @@ if (IQTREE_FLAGS MATCHES "mpi") 
	add_library(mympi utils/TreeCollection.cpp utils/ObjectStream.cpp)

endif()



if (NOT IQTREE_FLAGS MATCHES "single")

    if (BINARY32)

        link_directories(${PROJECT_SOURCE_DIR}/lib32)

    else()

        link_directories(${PROJECT_SOURCE_DIR}/lib)

    endif()

endif()



add_executable(iqtree

main/main.cpp

main/phyloanalysis.cpp

main/phyloanalysis.h

main/phylotesting.cpp

main/phylotesting.h

obsolete/parsmultistate.cpp

)



if(Backtrace_FOUND)
@@ -612,7 +621,7 @@ endif() 
set_target_properties(iqtree PROPERTIES OUTPUT_NAME "iqtree${EXE_SUFFIX}")



# strip the release build

if (CMAKE_BUILD_TYPE STREQUAL "Release" AND (GCC OR CLANG) AND NOT APPLE) # strip is not necessary for MSVC

if (NOT IQTREE_FLAGS MATCHES "nostrip" AND CMAKE_BUILD_TYPE STREQUAL "Release" AND (GCC OR CLANG) AND NOT APPLE) # strip is not necessary for MSVC

	if (WIN32)

		ADD_CUSTOM_COMMAND(TARGET iqtree POST_BUILD COMMAND strip $<TARGET_FILE:iqtree>)

	elseif (NOT APPLE)

alignment/alignment.cpp

+80 −74
Original line number Diff line number Diff line
@@ -532,18 +532,26 @@ int Alignment::readNexus(char *filename) { 
		return 0;

	}



    if (char_block->GetNTax() == 0) { char_block = data_block; }



    if (char_block->GetNTax() == 0) {

        outError("No data is given in the input file");

    if (data_block->GetNTax() == 0 && char_block->GetNTax() == 0) {

        outError("No DATA or CHARACTERS blocks found");

        return 0;

    }

    if (verbose_mode >= VB_DEBUG)

        char_block->Report(cout);





    if (char_block->GetNTax() > 0) {

        extractDataBlock(char_block);

        if (verbose_mode >= VB_DEBUG)

            char_block->Report(cout);

    } else {

        extractDataBlock(data_block);

        if (verbose_mode >= VB_DEBUG)

            data_block->Report(cout);

    }



    delete trees_block;

    delete char_block;

    delete data_block;

    delete assumptions_block;

    delete taxa_block;

    return 1;

}
@@ -567,8 +575,8 @@ int getDataBlockMorphStates(NxsCharactersBlock *data_block) { 
    char ch;

    int nstates = 0;



    for (site = 0; site < nsite; site++)

        for (seq = 0; seq < nseq; seq++) {

    for (seq = 0; seq < nseq; seq++)

        for (site = 0; site < nsite; site++) {

            int nstate = data_block->GetNumStates(seq, site);

            if (nstate == 0)

                continue;
@@ -580,23 +588,11 @@ int getDataBlockMorphStates(NxsCharactersBlock *data_block) { 
                else if (ch >= 'A' && ch <= 'Z')

                    ch = ch - 'A' + 11;

                else

                    outError(data_block->GetTaxonLabel(seq) + " has invalid state at site " + convertIntToString(site));

                    outError(data_block->GetTaxonLabel(seq) + " has invalid single state " + ch + " at site " + convertIntToString(site+1));

                if (ch > nstates) nstates = ch;

                continue;

            }

            for (int state = 0; state < nstate; state++) {

                ch = data_block->GetState(seq, site, state);

                if (!isalnum(ch)) continue;

                if (ch >= '0' && ch <= '9') ch = ch - '0' + 1;

                if (ch >= 'A' && ch <= 'Z') ch = ch - 'A' + 11;

                if (ch >= '0' && ch <= '9')

                    ch = ch - '0' + 1;

                else if (ch >= 'A' && ch <= 'Z')

                    ch = ch - 'A' + 11;

                else

                    outError(data_block->GetTaxonLabel(seq) + " has invalid state at site " + convertIntToString(site));

                if (ch > nstates) nstates = ch;

            }

            //cout << "NOTE: " << data_block->GetTaxonLabel(seq) << " has ambiguous state at site " << site+1 << " which is treated as unknown" << endl;

        }

    return nstates;

}
@@ -609,6 +605,9 @@ void Alignment::extractDataBlock(NxsCharactersBlock *data_block) { 
    char char_to_state[NUM_CHAR];

    char state_to_char[NUM_CHAR];

    

    if (!data_block->GetMatrix())

        outError("MATRIX command undeclared or invalid");

    

    NxsCharactersBlock::DataTypesEnum data_type = (NxsCharactersBlock::DataTypesEnum)data_block->GetDataType();

    if (data_type == NxsCharactersBlock::continuous) {

        outError("Continuous characters not supported");
@@ -674,14 +673,24 @@ void Alignment::extractDataBlock(NxsCharactersBlock *data_block) { 
                pat.push_back(STATE_UNKNOWN);

            else if (nstate == 1) {

                pat.push_back(char_to_state[(int)data_block->GetState(seq, site, 0)]);

            } else {

                ASSERT(data_type != NxsCharactersBlock::dna || data_type != NxsCharactersBlock::rna || data_type != NxsCharactersBlock::nucleotide);

            } else if (data_type == NxsCharactersBlock::dna || data_type == NxsCharactersBlock::rna || data_type == NxsCharactersBlock::nucleotide) {

                // 2018-06-07: correctly interpret ambiguous nucleotide

                char pat_ch = 0;

                for (int state = 0; state < nstate; state++) {

                    pat_ch |= (1 << char_to_state[(int)data_block->GetState(seq, site, state)]);

                }

                pat_ch += 3;

                pat.push_back(pat_ch);

            } else {

                // other ambiguous characters are treated as unknown

                stringstream str;

                str << "Sequence " << seq_names[seq] << " site " << site+1 << ": {";

                for (int state = 0; state < nstate; state++) {

                    str << data_block->GetState(seq, site, state);

                }

                str << "} treated as unknown character";

                outWarning(str.str());

                pat.push_back(STATE_UNKNOWN);

            }

        }

        num_gaps_only += addPattern(pat, site);
@@ -1103,6 +1112,7 @@ void Alignment::buildStateMap(char *map, SeqType seq_type) { 
        map[(unsigned char)'J'] = 22; // I or L

        map[(unsigned char)'*'] = STATE_UNKNOWN; // stop codon

        map[(unsigned char)'U'] = STATE_UNKNOWN; // 21st amino acid

        map[(unsigned char)'O'] = STATE_UNKNOWN; // 22nd amino acid



        return;

    case SEQ_MULTISTATE:
@@ -1191,6 +1201,7 @@ char Alignment::convertState(char state, SeqType seq_type) { 
		if (state == 'J') return 22;

        if (state == '*') return STATE_UNKNOWN; // stop codon

        if (state == 'U') return STATE_UNKNOWN; // 21st amino-acid

        if (state == 'O') return STATE_UNKNOWN; // 22nd amino-acid

        loc = strchr(symbols_protein, state);



        if (!loc) return STATE_INVALID; // unrecognize character
@@ -1287,13 +1298,13 @@ char Alignment::convertStateBack(char state) { 
    }

}



string Alignment::convertStateBackStr(char state) {

string Alignment::convertStateBackStr(StateType state) {

	string str;

    if (seq_type == SEQ_POMO)

        return string("POMO")+convertIntToString(state);

	if (seq_type != SEQ_CODON) {

		str = convertStateBack(state);

	} else {

    if (seq_type == SEQ_MULTISTATE)

        return " " + convertIntToString(state);

	if (seq_type == SEQ_CODON) {

        // codon data

        if (state >= num_states) return "???";

        assert(codon_table);
@@ -1301,7 +1312,10 @@ string Alignment::convertStateBackStr(char state) { 
        str = symbols_dna[state/16];

        str += symbols_dna[(state%16)/4];

        str += symbols_dna[state%4];

        return str;

	}

    // all other data types

    str = convertStateBack(state);

	return str;

}
@@ -1395,9 +1409,9 @@ SeqType Alignment::getSeqType(const char *sequence_type) { 
        user_seq_type = SEQ_PROTEIN;

    } else if (strncmp(sequence_type, "NT2AA", 5) == 0) {

        user_seq_type = SEQ_PROTEIN;

    } else if (strcmp(sequence_type, "NUM") == 0 || strcmp(sequence_type, "MORPH") == 0 || strcmp(sequence_type, "MULTI") == 0) {

    } else if (strcmp(sequence_type, "NUM") == 0 || strcmp(sequence_type, "MORPH") == 0) {

        user_seq_type = SEQ_MORPH;

    } else if (strcmp(sequence_type, "TINA") == 0) {

    } else if (strcmp(sequence_type, "TINA") == 0 || strcmp(sequence_type, "MULTI") == 0) {

        user_seq_type = SEQ_MULTISTATE;

    } else if (strncmp(sequence_type, "CODON", 5) == 0) {

        user_seq_type = SEQ_CODON;
@@ -1492,11 +1506,11 @@ int Alignment::buildPattern(StrVector &sequences, char *sequence_type, int nseq, 
            num_states = 20;

            nt2aa = true;

            cout << "Translating to amino-acid sequences with genetic code " << &sequence_type[5] << " ..." << endl;

        } else if (strcmp(sequence_type, "NUM") == 0 || strcmp(sequence_type, "MORPH") == 0 || strcmp(sequence_type, "MULTI") == 0) {

        } else if (strcmp(sequence_type, "NUM") == 0 || strcmp(sequence_type, "MORPH") == 0) {

            num_states = getMorphStates(sequences);

            if (num_states < 2 || num_states > 32) throw "Invalid number of states";

            user_seq_type = SEQ_MORPH;

        } else if (strcmp(sequence_type, "TINA") == 0) {

        } else if (strcmp(sequence_type, "TINA") == 0 || strcmp(sequence_type, "MULTI") == 0) {

            cout << "Multi-state data with " << num_states << " alphabets" << endl;

            user_seq_type = SEQ_MULTISTATE;

        } else if (strncmp(sequence_type, "CODON", 5) == 0) {
@@ -1623,7 +1637,7 @@ int Alignment::readPhylip(char *filename, char *sequence_type) { 
    string line;

    // remove the failbit

    in.exceptions(ios::badbit);

    bool tina_state = (sequence_type && strcmp(sequence_type,"TINA") == 0);

    bool tina_state = (sequence_type && (strcmp(sequence_type,"TINA") == 0 || strcmp(sequence_type,"MULTI") == 0));

    num_states = 0;



    for (; !in.eof(); line_num++) {
@@ -2924,6 +2938,8 @@ void convert_range(const char *str, int &lower, int &upper, int &step_size, char 
    //int d_save = d;

    upper = d;

    step_size = 1;

    // skip blank chars

    for (; *endptr == ' '; endptr++) {}

    if (*endptr != '-') return;



    // parse the upper bound of the range
@@ -2938,6 +2954,9 @@ void convert_range(const char *str, int &lower, int &upper, int &step_size, char 


    //lower = d_save;

    upper = d;

    // skip blank chars

    for (; *endptr == ' '; endptr++) {}



    if (*endptr != '\\') return;



    // parse the step size of the range
@@ -2977,7 +2996,7 @@ void extractSiteID(Alignment *aln, const char* spec, IntVector &site_id) { 
            for (i = lower; i <= upper; i+=step)

                site_id.push_back(i);

            if (*str == ',' || *str == ' ') str++;

            else break;

            //else break;

        }

        if (aln->seq_type == SEQ_CODON && nchars % 3 != 0)

            throw (string)"Range " + spec + " length is not multiple of 3 (necessary for codon data)";
@@ -3044,18 +3063,11 @@ void Alignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq 
    if (!spec) {

		// standard bootstrap

        int added_sites = 0;

		for (site = 0; site < nsite; site++) {

            int site_id;

            if (Params::getInstance().jackknife_prop == 0.0) {

                // bootstrap sampling with replacement

                site_id = random_int(nsite);

            } else {

                // jacknife without replacement

                if (random_double() < Params::getInstance().jackknife_prop)

                    continue;

                site_id = site;

            }

			int ptn_id = aln->getPatternID(site_id);

        IntVector sample;

        random_resampling(nsite, sample);

		for (site = 0; site < nsite; site++)

        for (int rep = 0; rep < sample[site]; rep++) {

			int ptn_id = aln->getPatternID(site);

			Pattern pat = aln->at(ptn_id);

            int nptn = getNPattern();

			addPattern(pat, added_sites);
@@ -3165,26 +3177,16 @@ void Alignment::createBootstrapAlignment(int *pattern_freq, const char *spec, in 
    if (Params::getInstance().jackknife_prop > 0.0 && spec)

        outError((string)"Unsupported jackknife with " + spec);



    if (!spec ||  strncmp(spec, "SCALE=", 6) == 0) {

    if (!spec) {



        if (spec) {

            double scale = convert_double(spec+6);

            nsite = (int)round(scale * nsite);

        }

        int nptn = getNPattern();



        if (nsite/8 < nptn || Params::getInstance().jackknife_prop > 0.0) {

            int orig_nsite = getNSite();

            for (site = 0; site < nsite; site++) {

                int site_id;

                if (Params::getInstance().jackknife_prop == 0.0)

                    site_id = random_int(orig_nsite, rstream);

                else {

                    if (random_double() < Params::getInstance().jackknife_prop)

                        continue;

                    site_id = site;

                }

                int ptn_id = getPatternID(site_id);

            IntVector sample;

            random_resampling(nsite, sample, rstream);

            for (site = 0; site < nsite; site++)

            for (int rep = 0; rep < sample[site]; rep++) {

                int ptn_id = getPatternID(site);

                pattern_freq[ptn_id]++;

            }

        } else {
@@ -3243,7 +3245,11 @@ void Alignment::createBootstrapAlignment(int *pattern_freq, const char *spec, in 
		}

	} else {

		// resampling sites within genes

        try {

            convert_int_vec(spec, site_vec);

        } catch (...) {

            outError("-bsam not allowed for non-partition model");

        }

		if (site_vec.size() % 2 != 0)

			outError("Bootstrap specification length is not divisible by 2");

		int part, begin_site = 0, out_site = 0;

alignment/alignment.h

+1 −1
Original line number Diff line number Diff line
@@ -258,7 +258,7 @@ public: 
	 * @param state internal state code

	 * @return user-readable state string

	 */

	string convertStateBackStr(char state);

	string convertStateBackStr(StateType state);



	/**

            get alignment site range from the residue range relative to a sequence

alignment/superalignment.cpp

+90 −19
Original line number Diff line number Diff line
@@ -799,6 +799,7 @@ void SuperAlignment::printSiteInfo(const char* filename) { 
    }

}



/*

void SuperAlignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq, const char *spec) {

	ASSERT(aln->isSuperAlignment());

	Alignment::copyAlignment(aln);
@@ -881,6 +882,71 @@ void SuperAlignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern 
	taxa_index = super_aln->taxa_index;

    countConstSite();

}

*/



void SuperAlignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq, const char *spec) {

    ASSERT(aln->isSuperAlignment());

    SuperAlignment *super_aln = (SuperAlignment*) aln;

    ASSERT(partitions.empty());

    name = aln->name;

    model_name = aln->model_name;

    sequence_type = aln->sequence_type;

    position_spec = aln->position_spec;

    aln_file = aln->aln_file;

    

    if (!spec) {

        // resampling sites within genes

        Alignment::copyAlignment(aln);

        partitions.reserve(super_aln->partitions.size());

        for (vector<Alignment*>::iterator it = super_aln->partitions.begin(); it != super_aln->partitions.end(); it++) {

            Alignment *boot_aln = new Alignment;

            if (pattern_freq) {

                IntVector part_pattern_freq;

                boot_aln->createBootstrapAlignment(*it, &part_pattern_freq);

                pattern_freq->insert(pattern_freq->end(), part_pattern_freq.begin(), part_pattern_freq.end());

            } else {

                boot_aln->createBootstrapAlignment(*it);

            }

            partitions.push_back(boot_aln);

        }

        taxa_index = super_aln->taxa_index;

        countConstSite();

    } else if (strcmp(spec, "GENE") == 0) {

        ASSERT(!pattern_freq);

        // resampling whole genes

        IntVector gene_freq;

        random_resampling(super_aln->partitions.size(), gene_freq);

        for (int i = 0; i < gene_freq.size(); i++)

            if (gene_freq[i] > 0) {

                Alignment *boot_aln = new Alignment;

                boot_aln->copyAlignment(super_aln->partitions[i]);

                if (gene_freq[i] > 1) {

                    for (auto it = boot_aln->begin(); it != boot_aln->end(); it++)

                        it->frequency *= gene_freq[i];

                    auto site_pattern = boot_aln->site_pattern;

                    for (int j = 1; j < gene_freq[i]; j++)

                        boot_aln->site_pattern.insert(boot_aln->site_pattern.end(), site_pattern.begin(), site_pattern.end());

                }

                partitions.push_back(boot_aln);

            }

        init();

    } else if (strcmp(spec, "GENESITE") == 0) {

        ASSERT(!pattern_freq);

        // resampling whole genes then sites within resampled genes

        IntVector gene_freq;

        random_resampling(super_aln->partitions.size(), gene_freq);

        for (int i = 0; i < gene_freq.size(); i++)

            for (int rep = 0; rep < gene_freq[i]; rep++) {

            Alignment *boot_aln = new Alignment;

            boot_aln->createBootstrapAlignment(super_aln->partitions[i]);

            boot_aln->name = boot_aln->name + "." + convertIntToString(rep);

            partitions.push_back(boot_aln);

        }

        init();

    } else {

        outError("Wrong -bsam, either -bsam GENE or -bsam GENESITE");

    }

}



void SuperAlignment::createBootstrapAlignment(IntVector &pattern_freq, const char *spec) {

	ASSERT(isSuperAlignment());
@@ -910,28 +976,29 @@ void SuperAlignment::createBootstrapAlignment(int *pattern_freq, const char *spe 
			nptn += (*it)->getNPattern();

		}

		memset(pattern_freq, 0, nptn * sizeof(int));

		for (int i = 0; i < partitions.size(); i++) {

			int part = random_int(partitions.size(), rstream);

        IntVector gene_freq;

        random_resampling(partitions.size(), gene_freq, rstream);

		for (int part = 0; part < partitions.size(); part++)

        for (int rep = 0; rep < gene_freq[part]; rep++){

			Alignment *aln = partitions[part];

			if (strncmp(spec,"GENESITE",8) == 0) {

				// then resampling sites in resampled gene

				for (int j = 0; j < aln->getNSite(); j++) {

					int ptn_id = aln->getPatternID(random_int(aln->getNPattern(), rstream));

                IntVector sample;

                random_resampling(aln->getNSite(), sample, rstream);

				for (int site = 0; site < sample.size(); site++)

                for (int rep2 = 0; rep2 < sample[site]; rep2++) {

					int ptn_id = aln->getPatternID(site);

					pattern_freq[ptn_id + part_pos[part]]++;

				}



			} else {

				for (int j = 0; j < aln->getNPattern(); j++)

					pattern_freq[j + part_pos[part]] += aln->at(j).frequency;

				for (int ptn = 0; ptn < aln->getNPattern(); ptn++)

					pattern_freq[ptn + part_pos[part]] += aln->at(ptn).frequency;

			}

		}

	} else {

		// resampling sites within genes

		int offset = 0;

		for (vector<Alignment*>::iterator it = partitions.begin(); it != partitions.end(); it++) {

            if (spec && strncmp(spec, "SCALE=", 6) == 0)

                (*it)->createBootstrapAlignment(pattern_freq + offset, spec, rstream);

            else

            (*it)->createBootstrapAlignment(pattern_freq + offset, NULL, rstream);

			offset += (*it)->getNPattern();

		}
@@ -1148,21 +1215,25 @@ Alignment *SuperAlignment::concatenateAlignments(set<int> &ids) { 
    int site = 0;

    for (it = ids.begin(); it != ids.end(); it++) {

    	int id = *it;

		string taxa_set;

        Pattern taxa_pat = getPattern(id);

        taxa_set.insert(taxa_set.begin(), taxa_pat.begin(), taxa_pat.end());

        // 2018-08-23: important bugfix in v1.6: taxa_set has wrong correspondance

		//string taxa_set;

        //Pattern taxa_pat = getPattern(id);

        //taxa_set.insert(taxa_set.begin(), taxa_pat.begin(), taxa_pat.end());

    	for (Alignment::iterator it = partitions[id]->begin(); it != partitions[id]->end(); it++) {

    		Pattern pat;

    		int part_seq = 0;

    		//int part_seq = 0;

    		for (int seq = 0; seq < union_taxa.size(); seq++)

    			if (union_taxa[seq] == 1) {

    				char ch = aln->STATE_UNKNOWN;

    				if (taxa_set[seq] == 1) {

    					ch = (*it)[part_seq++];

    				}

                    int seq_part = taxa_index[seq][id];

                    if (seq_part >= 0)

                        ch = (*it)[seq_part];

                    //if (taxa_set[seq] == 1) {

                    //    ch = (*it)[part_seq++];

                    //}

    				pat.push_back(ch);

    			}

    		ASSERT(part_seq == partitions[id]->getNSeq());

    		//ASSERT(part_seq == partitions[id]->getNSeq());

    		aln->addPattern(pat, site, (*it).frequency);

    		// IMPORTANT BUG FIX FOLLOW

    		int ptnindex = aln->pattern_index[pat];

main/main.cpp

+3 −3
Original line number Diff line number Diff line
@@ -54,7 +54,7 @@ 
#include "phyloanalysis.h"

#include "tree/matree.h"

//#include "ngs.h"

//#include "parsmultistate.h"

#include "obsolete/parsmultistate.h"

//#include "gss.h"

#include "alignment/maalignment.h" //added by MA

#include "tree/ncbitree.h"
@@ -2540,8 +2540,8 @@ int main(int argc, char *argv[]) { 
	// call the main function

	if (Params::getInstance().tree_gen != NONE) {

		generateRandomTree(Params::getInstance());

//	} else if (Params::getInstance().do_pars_multistate) {

//		doParsMultiState(Params::getInstance());

    } else if (Params::getInstance().do_pars_multistate) {

        doParsMultiState(Params::getInstance());

	} else if (Params::getInstance().rf_dist_mode != 0) {

		computeRFDist(Params::getInstance());

	} else if (Params::getInstance().test_input != TEST_NONE) {