--- /dev/null
+/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
+
+/*********************************************************************
+ * Clustal Omega - Multiple sequence alignment
+ *
+ * Copyright (C) 2010 University College Dublin
+ *
+ * Clustal-Omega is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) any later version.
+ *
+ * This file is part of Clustal-Omega.
+ *
+ ********************************************************************/
+
+/*
+ * RCS $Id: hhhalfalignment-C.h 227 2011-03-28 17:03:09Z fabian $
+ */
+
+// hhfullalignment.C
+
+#ifndef MAIN
+#define MAIN
+#include <iostream> // cin, cout, cerr
+#include <fstream> // ofstream, ifstream
+#include <stdio.h> // printf
+#include <stdlib.h> // exit
+#include <string> // strcmp, strstr
+#include <math.h> // sqrt, pow
+#include <limits.h> // INT_MIN
+#include <float.h> // FLT_MIN
+#include <time.h> // clock
+#include <ctype.h> // islower, isdigit etc
+using std::ios;
+using std::ifstream;
+using std::ofstream;
+using std::cout;
+using std::cerr;
+using std::endl;
+#include "util-C.h" // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc.
+#include "list.h" // list data structure
+#include "hash.h" // hash data structure
+#include "hhdecl-C.h" // constants, class
+#include "hhutil-C.h" // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc.
+#include "hhhmm.h" // class HMM
+#include "hhalignment.h" // class Alignment
+#include "hhhit.h"
+#endif
+
+/////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////
+// Methods of class HalfAlignment
+/////////////////////////////////////////////////////////////////////////////
+/////////////////////////////////////////////////////////////////////////////
+
+
+
+/////////////////////////////////////////////////////////////////////////////
+// Constructor
+HalfAlignment::HalfAlignment(int maxseqdis)
+{
+ n=0;
+ sname=seq=NULL;
+ nss_dssp = nss_pred = nss_conf = nsa_dssp = ncons= -1;
+ h = new(int[maxseqdis]); //h[k] = next position of sequence k to be written
+ s = new(char*[maxseqdis]); //s[k][h] = character in column h, sequence k of output alignment
+ l = new(int*[maxseqdis]); //counts non-gap residues: l[k][i] = index of last residue AT OR BEFORE match state i in seq k
+ m = new(int*[maxseqdis]); //counts positions: m[k][i] = position of match state i in string seq[k]
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+// Destructor
+HalfAlignment::~HalfAlignment()
+{
+ Unset();
+ delete[] h; h = NULL;
+ delete[] s; s = NULL;
+ delete[] l; l = NULL;
+ delete[] m; m = NULL;
+}
+
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Free memory in HalfAlignment arrays s[][], l[][], and m[][]
+ */
+void
+HalfAlignment::Unset()
+{
+ // Free memory for alignment characters and residue counts
+ for (int k=0; k<n; k++)
+ {
+ delete[] s[k]; s[k] = NULL;
+ delete[] l[k]; l[k] = NULL;
+ delete[] m[k]; m[k] = NULL;
+ }
+ n=0;
+ sname=seq=NULL;
+ nss_dssp = nss_pred = nss_conf = nsa_dssp = ncons= -1;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Prepare a2m/a3m alignment:
+ * Calculate l[k][i] (residue indices) and m[k][i] (position in seq[k])
+*/
+void
+HalfAlignment::Set(char* name, char** seq_in, char** sname_in, int n_in, int L_in, int n1, int n2, int n3, int n4, int nc, int L_in2/*<--FS*/)
+{
+ int i; /* counts match states in seq[k] */
+ int ll; /* counts residues LEFT from or at current position in seq[k] */
+ int mm; /* counts postions in string seq[k] */
+ int k; /* counts sequences */
+ char c;
+ char warned=0;
+
+ nss_dssp=n1; nss_pred=n2; nss_conf=n3; nsa_dssp=n4; ncons=nc;
+ seq=seq_in; /* flat copy of sequences */
+ sname=sname_in; /* flat copy of sequence names */
+ n=n_in;
+ L=L_in;
+ pos=0;
+
+ /* Allocate memory for alignment characters and residue counts */
+ for (k=0; k<n; k++) {
+ s[k]=new char[LINELEN];
+ l[k]=new int[L+10+L_in2/*<--FS*/];
+ m[k]=new int[L+10+L_in2/*<--FS*/];
+ if (!s[k] || !l[k] || !m[k]) MemoryError("space for formatting HMM-HMM alignment");
+ h[k]=0; //starting positions in alignment = 0
+ } /* k <= 0 < n (= n_in) */
+
+ for (k=0; k<n; k++) {
+ m[k][0]=0; // 0'th match state (virtual) is begin state at 0
+ //if k is consensus sequence
+ if (k==nc) {
+ for (i=1; i<=L; i++) m[k][i]=l[k][i]=i;
+ m[k][L+1]=l[k][L+1]=L;
+ continue;
+ }
+ i=1; mm=1; ll=1;
+ while ((c=seq[k][mm]))
+ {
+ if (MatchChr(c)==c) //count match/delete states
+ {
+ l[k][i]=ll;
+ m[k][i]=mm;
+ i++;
+ }
+ if (WordChr(c)) ll++; //index of next residue
+ mm++;
+ }
+ l[k][i]=ll-1; //set l[k][L+1] eq number of residues in seq k (-1 since there is no residue at L+1st match state)
+ m[k][i]=mm; //set m[k][L+1]
+
+ if ((i-1)!=L && !warned)
+ {
+ cerr<<"Warning: sequence "<<sname[k]<<" in HMM "<<name<<" has "<<i<<" match states but should have "<<L<<"\n";
+ warned=1;
+ }
+ } /* k <= 0 < n (= n_in) */
+ //DEBUG
+ if (v>=5)
+ {
+ printf(" i chr m l\n");
+ for(i=0;i<=L+1;i++) printf("%3i %1c %3i %3i\n",i,seq[0][m[0][i]],m[0][i],l[0][i]);
+ printf("\n");
+ }
+
+} /*** end HalfAlignment::Set() ***/
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Fill in insert states following match state i (without inserting '.' to fill up)
+ */
+void
+HalfAlignment::AddInserts(int i)
+{
+ for (int k=0; k<n; k++) // for all sequences...
+ for (int mm=m[k][i]+1; mm<m[k][i+1]; mm++) // for all inserts between match state i and i+1...
+ s[k][h[k]++]=seq[k][mm]; // fill inserts into output alignment s[k]
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Fill up alignment with gaps '.' to generate flush end (all h[k] equal)
+ */
+void
+HalfAlignment::FillUpGaps()
+{
+ int k; //counts sequences
+ pos=0;
+
+ // Determine max position h[k]
+ for (k=0; k<n; k++) pos = imax(h[k],pos);
+
+ // Fill in gaps up to pos
+ for (k=0; k<n; k++)
+ {
+ for (int hh=h[k]; hh<pos; hh++) s[k][hh]='.';
+ h[k]=pos;
+ }
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Fill in insert states following match state i and fill up gaps with '.'
+ */
+void
+HalfAlignment::AddInsertsAndFillUpGaps(int i)
+{
+ AddInserts(i);
+ FillUpGaps();
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Add gap column '.'
+ */
+void
+HalfAlignment::AddChar(char c)
+{
+ for (int k=0; k<n; k++) s[k][h[k]++]=c;
+ pos++;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Add match state column i as is
+ */
+void
+HalfAlignment::AddColumn(int i)
+{
+ for (int k=0; k<n; k++) s[k][h[k]++]=seq[k][m[k][i]];
+ pos++;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Add match state column i as insert state
+ */
+void
+HalfAlignment::AddColumnAsInsert(int i)
+{
+ char c;
+ for (int k=0; k<n; k++)
+ if ((c=seq[k][m[k][i]])!='-' && (c<'0' || c>'9'))
+ s[k][h[k]++]=InsertChr(c);
+ pos++;
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Remove all characters c from template sequences
+ */
+void
+HalfAlignment::RemoveChars(char c)
+{
+ int k,h,hh;
+ for (k=0; k<n; k++)
+ {
+ for (h=hh=0; h<pos; h++)
+ if (s[k][h]!=c) s[k][hh++]=s[k][h];
+ s[k][++hh]='\0';
+ }
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Transform alignment sequences from A3M to A2M (insert ".")
+ */
+void
+HalfAlignment::BuildFASTA()
+{
+ AddInserts(0);
+ FillUpGaps();
+ for (int i=1; i<=L; i++)
+ {
+ AddColumn(i);
+ AddInserts(i);
+ FillUpGaps();
+ }
+ ToFASTA();
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Transform alignment sequences from A3M to A2M (insert ".")
+ */
+void
+HalfAlignment::BuildA2M()
+{
+ AddInserts(0);
+ FillUpGaps();
+ for (int i=1; i<=L; i++)
+ {
+ AddColumn(i);
+ AddInserts(i);
+ FillUpGaps();
+ }
+ AddChar('\0');
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Transform alignment sequences from A3M to A2M (insert ".")
+ */
+void
+HalfAlignment::BuildA3M()
+{
+ AddInserts(0);
+ for (int i=1; i<=L; i++)
+ {
+ AddColumn(i);
+ AddInserts(i);
+ }
+ AddChar('\0');
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Transform alignment sequences from A2M to FASTA ( lowercase to uppercase and '.' to '-')
+ */
+void
+HalfAlignment::ToFASTA()
+{
+ for (int k=0; k<n; k++)
+ {
+ uprstr(s[k]);
+ strtr(s[k],".","-");
+ }
+}
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Align query (HalfAlignment) to template (i.e. hit) match state structure
+ */
+void
+HalfAlignment::AlignToTemplate(Hit& hit)
+{
+ int i,j;
+ int step; // column of the HMM-HMM alignment (first:nstep, last:1)
+ char state;
+
+ if(0) { //par.loc==0) { //////////////////////////////////////////// STILL NEEDED??
+ // If in global mode: Add part of alignment before first MM state
+ AddInserts(0); // Fill in insert states before first match state
+ for (i=1; i<hit.i[hit.nsteps]; i++)
+ {
+ AddColumnAsInsert(i);
+ AddInserts(i);
+ if (par.outformat<=2) FillUpGaps();
+ }
+ }
+
+ // Add endgaps (First state must be an MM state!!)
+ for (j=1; j<hit.j[hit.nsteps]; j++)
+ {
+ AddChar('-');
+ }
+
+ // Add alignment between first and last MM state
+ for (step=hit.nsteps; step>=1; step--)
+ {
+ state = hit.states[step];
+ i = hit.i[step];
+
+ switch(state)
+ {
+ case MM: //MM pair state (both query and template in Match state)
+ AddColumn(i);
+ AddInserts(i);
+ break;
+ case DG: //D- state
+ case MI: //MI state
+ AddColumnAsInsert(i);
+ AddInserts(i);
+ break;
+ case GD: //-D state
+ case IM: //IM state
+ AddChar('-');
+ break;
+ }
+ if (par.outformat<=2) FillUpGaps();
+
+ }
+
+ if(0) { //par.loc==0) { //////////////////////////////////////////// STILL NEEDED??
+
+ // If in global mode: Add part of alignment after last MM state
+ for (i=hit.i[1]+1; i<=L; i++)
+ {
+ AddColumnAsInsert(i);
+ AddInserts(i);
+ if (par.outformat==2) FillUpGaps();
+ }
+ }
+
+ // Add endgaps
+ for (j=hit.j[1]+1; j<=hit.L; j++)
+ {
+ AddChar('-');
+ }
+
+ // Add end-of-string character
+ AddChar('\0');
+}
+
+
+/////////////////////////////////////////////////////////////////////////////////////
+/**
+ * @brief Write the a2m/a3m alignment into alnfile
+ */
+void
+HalfAlignment::Print(char* alnfile)
+{
+ int k; //counts sequences
+ int omitted=0; // counts number of sequences with no residues in match states
+ FILE *outf;
+ if (strcmp(alnfile,"stdout"))
+ {
+ if (par.append) outf=fopen(alnfile,"a"); else outf=fopen(alnfile,"w");
+ if (!outf) OpenFileError(alnfile);
+ }
+ else
+ outf = stdout;
+ if (v>=3) cout<<"Writing alignment to "<<alnfile<<"\n";
+
+ for (k=0; k<n; k++)
+ {
+ // Print sequence only if it contains at least one residue in a match state
+ if (1) //strpbrk(s[k],"ABCDEFGHIKLMNPQRSTUVWXYZ1234567890"))
+ {
+ fprintf(outf,">%s\n",sname[k]);
+ fprintf(outf,"%s\n",s[k]);
+ } else {
+ omitted++;
+ if (v>=3) printf("%-14.14s contains no residue in match state. Omitting sequence\n",sname[k]);
+ }
+ }
+ if (v>=2 && omitted) printf("Omitted %i sequences in %s which contained no residue in match state\n",omitted,alnfile);
+ fclose(outf);
+}
+
+
+/** EOF hhhalfalignment-C.h **/
git://source.jalview.org / jabaws.git / blobdiff