Merge pull request #68 from iqtree/huaiyan

Update sequence type detection function

Author: bqminh

alignment/alignment.cpp

@@ -1583,53 +1583,70 @@ void Alignment::regroupSitePattern(int groups, IntVector& site_group)
 	@return the data type of the input sequences
 */
 SeqType Alignment::detectSequenceType(StrVector &sequences) {
+    size_t num_proper_nuc = 0;
     size_t num_nuc   = 0;
-    size_t num_ungap = 0;
+    size_t num_aa    = 0;
     size_t num_bin   = 0;
-    size_t num_alpha = 0;
     size_t num_digit = 0;
+    size_t num_alpha = 0;
     double detectStart = getRealTime();
     size_t sequenceCount = sequences.size();
+    std::unordered_set<char> proper_nucleotides = {'A', 'C', 'G', 'T', 'U'};
+    std::unordered_set<char> nucleotides = {'A', 'C', 'G', 'T', 'U', 'R', 'Y', 'W', 'S', 'M', 'K', 'B', 'H', 'D', 'V', 'N', 'X'};
+    std::unordered_set<char> proper_amino_acids = {'A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I', 'L', 'K', 'M', 'F', 'P', 'S',  'T', 'W', 'Y', 'V'};
+    std::unordered_set<char> binaries = {'0', '1'};
+//    std::unordered_set<char> gap_miss = {'?', '-', '.', '~'};
+    
 #ifdef _OPENMP
-#pragma omp parallel for reduction(+:num_nuc,num_ungap,num_bin,num_alpha,num_digit)
+#pragma omp parallel for reduction(+:num_proper_nuc, num_nuc, num_aa, num_bin, num_digit, num_alpha)
 #endif
     for (size_t seqNum = 0; seqNum < sequenceCount; ++seqNum) {
         auto start = sequences.at(seqNum).data();
         auto stop  = start + sequences.at(seqNum).size();
         for (auto i = start; i!=stop; ++i) {
-            if ((*i) == 'A' || (*i) == 'C' || (*i) == 'G' || (*i) == 'T' || (*i) == 'U') {
-                ++num_nuc;
-                ++num_ungap;
-                continue;
-            }
-            if ((*i)=='?' || (*i)=='-' || (*i) == '.' ) {
-                continue;
-            }
-            if (*i != 'N' && *i != 'X' &&  (*i) != '~') {
-                num_ungap++;
-                if (isdigit(*i)) {
-                    num_digit++;
-                    if ((*i) == '0' || (*i) == '1') {
-                        num_bin++;
-                    }
-                }
-            }
-            if (isalpha(*i)) {
+//            if (gap_miss.find(*i) != gap_miss.end()) {
+//                continue;
+//            }
+            if (proper_nucleotides.find(*i) != proper_nucleotides.end())
+                num_proper_nuc++;
+            if (nucleotides.find(*i) != nucleotides.end())
+                num_nuc++;
+            if (proper_amino_acids.find(*i) != proper_amino_acids.end())
+                num_aa++;
+            if (binaries.find(*i) != binaries.end())
+                num_bin++;
+            if (isdigit(*i))
+                num_digit++;
+            if (isalpha(*i))
                 num_alpha++;
-            }
         }
     }
     if (verbose_mode >= VB_MED) {
         cout << "Sequence Type detection took " << (getRealTime()-detectStart) << " seconds." << endl;
     }
-    if (((double)num_nuc) / num_ungap > 0.9)
-        return SEQ_DNA;
-    if (num_bin == num_ungap) // For binary data, only 0, 1, ?, -, . can occur
-        return SEQ_BINARY;
-    if (((double)num_alpha + num_nuc) / num_ungap > 0.9)
-        return SEQ_PROTEIN;
-    if (((double)(num_alpha + num_digit + num_nuc)) / num_ungap > 0.9)
+    if (num_digit == 0) {
+        if (num_alpha < 10) // two few occurences to decide
+            return SEQ_UNKNOWN;
+        if (num_nuc == num_alpha) {
+            if ((double)num_proper_nuc / num_alpha > 0.9) {
+                return SEQ_DNA;      
+            } else {
+                return SEQ_UNKNOWN;
+            }
+        }       
+        // likely protein if there are only letters
+        // TODO: check if this condition is OK
+        if ((double)num_aa / num_alpha > 0.9)
+            return SEQ_PROTEIN;
+        else
+            return SEQ_UNKNOWN;
+    } else if (num_digit >= 10) {
+        // there are some digit(s) in the data
+        if (num_bin == (num_digit+num_alpha)) // For binary data, only 0, 1, ?, -, . can occur
+            return SEQ_BINARY;
         return SEQ_MORPH;
+    }
+    // can't decide
     return SEQ_UNKNOWN;
 }
 

main/phyloanalysis.cpp

@@ -1364,7 +1364,6 @@ void reportSubstitutionProcess(ostream &out, Params &params, IQTree &tree)
 
         PhyloSuperTree *stree = (PhyloSuperTree*) &tree;
         PhyloSuperTree::iterator it;
-        it = stree->begin();
 
         int part;
 
@@ -1479,6 +1478,11 @@ void reportPhyloAnalysis(Params &params, IQTree &tree, ModelCheckpoint &model_in
 
         out << "SEQUENCE ALIGNMENT" << endl << "------------------" << endl
                 << endl;
+        
+        if (tree.aln->sequence_type.empty()) {
+            out << "NOTE: Alignment sequence type is auto-detected. If in doubt, specify it via -st option." << endl;
+        }
+        
         if (tree.isSuperTree()) {
       // TODO DS: Changes may be needed here for PoMo.
             out << "Input data: " << tree.aln->getNSeq()+tree.removed_seqs.size() << " taxa with "

model/partitionmodel.cpp

@@ -466,7 +466,7 @@ double PartitionModel::computeMarginalLh() {
                     sub_tree2->copyTree(tree2);
                 }
 
-                if (inter_seqs_id.size() == 2) {
+                if (inter_seqs_id.size() == 2 && tree2->getModel()->isReversible()) {
                     // if only two seqs in the subset
                     // add a gappy seq two the sub_aln
                     string gappy_seq = "gappy_seq";
@@ -565,113 +565,8 @@ double PartitionModel::computeMarginalLh() {
                 sub_tree2->aln = NULL;
                 delete sub_tree2;
                 delete[] ptn_lh_array;
-                /*
-            } else if (inter_seqs_id.size() == 2) {
-                //AlignmentPairwise aln_pair(tree2, inter_seqs_id[0], inter_seqs_id[1]);
-                //double site_lh = aln_pair.computeFunction(1);
-
-                //compute the distance between the tree nodes of the two "overlapping" taxa
-                auto it = inter_seqs.begin();
-                string inter_seq1 = *it;
-                Node *node1 = tree2->findLeafName(inter_seq1);
-                ++it;
-                string inter_seq2 = *it;
-                Node *node2 = tree2->findLeafName(inter_seq2);
-                double branch_len = tree2->pairDist(node1, node2, node1);
-
-                //compute the transition probobility matrix based on tree2
-                RateHeterogeneity *site_rate = tree2->getRate();
-                int ncat = site_rate->getNDiscreteRate();
-                int trans_size = tree2->getModel()->getTransMatrixSize();
-                double *trans_matrix = new double[trans_size];
-                double *sum_trans_matrix = new double[trans_size];
-
-                if (tree2->getModelFactory()->site_rate->getGammaShape() == 0.0)
-                    tree2->getModelFactory()->computeTransMatrix(branch_len, sum_trans_matrix);
-                else {
-                    tree2->getModelFactory()->computeTransMatrix(branch_len * site_rate->getRate(0), sum_trans_matrix);
-                    for (int cat = 1; cat < ncat; cat++) {
-                        tree2->getModelFactory()->computeTransMatrix(branch_len * site_rate->getRate(cat), trans_matrix);
-                        for (int i = 0; i < trans_size; i++)
-                            sum_trans_matrix[i] += trans_matrix[i];
-                    }
-                }
-
-                // take the log for the transition probobility matrix
-                vector<double> log_sum_trans_matrix(trans_size);
-                for (int n = 0; n < trans_size; n++) {
-                    log_sum_trans_matrix[n] = log(sum_trans_matrix[n]);
-                }
-
-                for (int l = 0; l < tree1_nsite; l++) {
-                    double site_lh = 0.0;
-                    Pattern p = tree1_aln->at(tree1_aln->getPatternID(l));
-                    int state1 = p[inter_seqs_id[0]];
-                    int state2 = p[inter_seqs_id[1]];
-
-                    // compute the site log-likelihood of the two "overlaping" sequences.
-                    if (state1 < n_states && state2 < n_states) {
-                        site_lh = log_sum_trans_matrix[state1 * n_states + state2] + log_state_freq[state1];
-                    } else if (state1 >= n_states && state2 < n_states) {
-                        // compute ambiguous frequencies
-                        int cstate1 = state1 - n_states + 1;
-                        double amb_lh = 0.0;
-                        for (int m = 0; m < n_states; m++) {
-                            if ((cstate1) & (1 << m)) {
-                                amb_lh = sum_trans_matrix[m * n_states + state2] * state_freq[m];
-                                site_lh *= amb_lh;
-                            }
-                        }
-                        site_lh = log(site_lh);
-                    } else if (state2 >= n_states && state1 < n_states) {
-                        // compute ambiguous frequencies
-                        int cstate2 = state2-n_states+1;
-                        double amb_lh = 0.0;
-                        for (int m = 0; m < n_states; m++) {
-                            if ((cstate2) & (1 << m)) {
-                                amb_lh = sum_trans_matrix[state1 * n_states + m] * state_freq[state1];
-                                site_lh *= amb_lh;
-                            }
-                        }
-                        site_lh = log(site_lh);
-                    } else {
-                        // compute ambiguous frequencies
-                        int cstate1 = state1-n_states+1;
-                        int cstate2 = state2-n_states+1;
-                        double amb_lh = 0.0;
-                        for (int m = 0; m < n_states; m++) {
-                            for (int n = 0; n < n_states; n++) {
-                                if (((cstate1) & (1 << m)) && ((cstate2) & (1 << n))) {
-                                    amb_lh = sum_trans_matrix[m * n_states + n] * state_freq[m];
-                                    site_lh *= amb_lh;
-                                }
-                            }
-                        }
-                        site_lh = log(site_lh);
-                    }
-
-                    for (int missing_id: missing_seqs_id) {
-                        int char_id = p[missing_id];
-                        if (char_id  < n_states) {
-                            site_lh += log_state_freq[char_id];
-                        } else {
-                            // compute ambiguous frequencies
-                            int cstate = char_id-n_states+1;
-                            double amb_freq = 0;
-                            for (int m = 0; m < n_states; m++) {
-                                if ((cstate) & (1 << m)) {
-                                    amb_freq += state_freq[m];
-                                }
-                            }
-                            site_lh += log(amb_freq);
-                        }
-                    }
-                    lh_array[tree1_nsite * k + l] = site_lh;
-                }
-                delete[] trans_matrix;
-                delete[] sum_trans_matrix;
-                 */
-            } else {
+            } else if (tree2->getModel()->isReversible()) {
+                // case when the intersection of taxon sets is 1
                 for (int l = 0; l < tree1_nsite; l++) {
                     double site_lh = 0.0;
                     Pattern p = tree1_aln->at(tree1_aln->getPatternID(l));
@@ -705,6 +600,9 @@ double PartitionModel::computeMarginalLh() {
                     }
                     lh_array[tree1_nsite * k + l] = site_lh;
                 }
+            } else {
+                // intersection has only 1 taxon and non-reversible model
+                outError("mAIC calculation doesn't work yet for intersection of 1 taxon and non-reversible model");
             }
             delete[] state_freq;
         }