% read consensus and build the model
consV = fastaread ('consV.fa');
hmmLen = sum(consV.Sequence ~= '-');
hmmInit   = hmmprofstruct (147, 'ALPHABET', 'AA');
hmmGlobin  = hmmprofestimate (hmmInit, msaGlobins);

% Read query seqeunces

query = fastaread('query_hmm.fa');
pfamGlobin = gethmmprof(42);

fprintf(1, "\nResults \n");
fprintf(1, "No. \t \t Simple \t \t \t PFAM \t  \n");
for i = 1:length(query)
   score_simple_model(i) =  hmmprofalign (hmmGlobin, query(i).Sequence);
   score_pfam_model(i) =  hmmprofalign (pfamGlobin, query(i).Sequence);
   if score_simple_model(i) > 0 
       class_simple_model{i} = 'globin';
   else
       class_simple_model{i} = 'non-globin';
   end
 
   if score_pfam_model(i) > 0 
       class_pfam_model{i} = 'globin';
   else
       class_pfam_model{i} = 'non-globin';
   end
   
   if strcmp(class_simple_model{i}, class_pfam_model{i})
       match{i} = 'MATCH';
   else
       match{i} = 'No';
   end
   fprintf(1, "%d\t  %6.3f \t %s \t %6.3f \t %s \t %s\n", i, ...
        score_simple_model(i), class_simple_model{i}, ...
        score_pfam_model(i), class_pfam_model{i}, match{i});
       
end

result = table(score_simple_model', class_simple_model', ...
                    score_pfam_model', class_pfam_model', match');
result.Properties.VariableNames = {'score_simple', 'class_simple', ...
               'score_pfam', 'class_pfam', 'Results_Matched'};
result
