22 files changed, 181 insertions, 71 deletions

diff --git a/src/assign/assign_protein_type.c b/src/assign/assign_protein_type.c
new file mode 100644
index 0000000..0a11c23
--- a/dev/null
+++ b/src/assign/assign_protein_type.c
@@ -0,0 +1,253 @@
+#include "assign_protein_type.h"
+#include "error/check_ncbi_error.h"
+#include "error/check_h5_error.h"
+#include "model/gi_type_data.h"
+#include "model/gi_type_data_init.h"
+#include "model/sequence_data.h"
+#include "model/sequence_data_init.h"
+#include <ncbi.h>
+#include <readdb.h>
+#include <blast.h>
+#include <salpacc.h>
+#include <stdbool.h>
+#include <hdf5_hl.h>
+#include <error.h>
+
+/*
+ * todo: Write results to an independent table by GI so that they can
+ * be reused when the input tables influenza.faa and influenza_aa.dat
+ * are replaced with new data.  Since GI values are unique and
+ * persistent the assignment of a protein type should stay valid
+ * unless the list of reference sequences is updated.
+ *
+ * todo: Save the species type as a new field { A, B, C }.
+ */
+
+/*
+ * BLAST database containing all of the influenza protein sequences.
+ */
+#define SEQDB "/home/don/exp004/influenzadb/influenzadb"
+
+/*
+ * BLAST reference database of prototypical protein types.
+ */
+#define REFDB "/home/don/exp004/influenzadb/proteinnames"
+
+#define BUFFER_LEN 50
+
+void
+assign_protein_type (hid_t file_id)
+{
+  /*
+   * Iterate through the records for which no protein type has been
+   * assigned.  Create a BioSeq Pointer to the data and then use this
+   * as input to the BLAST run against the reference database of
+   * prototypical sequence records.
+   */
+
+  /* 
+   * Make BLAST messages reportable.
+   */
+  ErrSetMessageLevel (SEV_WARNING);
+
+  /*
+   * Open the BLAST sequence database.
+   */
+  ReadDBFILEPtr seqdb = readdb_new (SEQDB, true);
+
+  /*
+   * Get default BLAST options.
+   */
+  BLAST_OptionsBlkPtr options = BLASTOptionNew ("blastp", true);
+
+  /*
+   * Set to single hit mode.
+   */
+  options->window_size = 0;
+  options->multiple_hits_only = false;
+
+  ValNodePtr error_returns = NULL;
+
+  /*
+   * Read the data from HDF5 influenza.faa.
+   */
+  hsize_t faa_nfields;
+  hsize_t faa_nrecords;
+  herr_t status = H5TBget_table_info (file_id, "influenza.faa", &faa_nfields,
+                                      &faa_nrecords);
+  if (status < 0)
+    check_h5_error (status, __FILE__, __LINE__);
+
+  sequence_data* faa_buf = malloc (sizeof(sequence_data) * faa_nrecords);
+
+  size_t faa_size;
+  size_t faa_offset[SEQUENCE_DATA_FIELD_NUM];
+  size_t faa_sizes[SEQUENCE_DATA_FIELD_NUM];
+  hid_t faa_field_type[SEQUENCE_DATA_FIELD_NUM];
+  sequence_data_init (&faa_size, faa_offset, faa_sizes, faa_field_type);
+
+  status = H5TBread_table (file_id, "influenza.faa", faa_size, faa_offset,
+                           faa_sizes, faa_buf);
+  if (status < 0)
+    check_h5_error (status, __FILE__, __LINE__);
+
+  /*
+   * Read the data from HDF5 gi_type_data.
+   */
+  hsize_t gi_nfields;
+  hsize_t gi_nrecords;
+  status = H5TBget_table_info (file_id, "gi_type_data", &gi_nfields,
+                               &gi_nrecords);
+  if (status < 0)
+    check_h5_error (status, __FILE__, __LINE__);
+  
+  gi_type_data* gi_buf = malloc (sizeof(gi_type_data) * gi_nrecords);
+
+  size_t gi_size;
+  size_t gi_offset[GI_TYPE_DATA_FIELD_NUM];
+  size_t gi_sizes[GI_TYPE_DATA_FIELD_NUM];
+  hid_t gi_field_type[GI_TYPE_DATA_FIELD_NUM];
+  gi_type_data_init (&gi_size, gi_offset, gi_sizes, gi_field_type);
+
+  status = H5TBread_table (file_id, "gi_type_data", gi_size, gi_offset,
+                           gi_sizes, gi_buf);
+  if (status < 0)
+    check_h5_error (status, __FILE__, __LINE__);
+
+  /*
+   * Assign protein types to records for which the field is empty.
+   */
+  printf ("Records to process: %i\n", (int)faa_nrecords);
+  bool updates_pending = false;
+  for (int i = 0; i < faa_nrecords; i++)
+    {
+
+      if (dst_buf[i].protein_type[0] == '\0') {
+        /*
+         * Read the sequence from the database by GI.
+         */
+        Int4 sequence_number = readdb_gi2seq (seqdb, dst_buf[i].gi, NULL);
+        BioseqPtr bsp = readdb_get_bioseq (seqdb, sequence_number);
+        if (bsp == NULL) 
+          {
+            error_at_line (EXIT_FAILURE, 0, __FILE__, __LINE__,
+                           "Unable to find BLAST record for gi|%i.  Ensure the BLAST "
+                           "database is up-to-date with the HDF5 record set.  See the "
+                           "BLAST formatdb.log file for details.\n",
+                           dst_buf[i].gi);
+          }
+        
+        SeqAlignPtr seqalign = BioseqBlastEngine (bsp,
+                                                  "blastp",
+                                                  REFDB,
+                                                  options,
+                                                  NULL,
+                                                  &error_returns,
+                                                  NULL);
+
+        /*
+         * BLAST reported an error.  Write it out and continue processing.
+         */
+        if (error_returns != NULL)
+          {
+            printf ("Warning: An error has been reported by the NCBI Toolkit API "
+                    "for sequence gi|%i: %s",
+                    dst_buf[i].gi, BlastErrorToString (error_returns));
+
+          }
+
+        /*
+         * A hit was found.  Record the first hit as the protein type.
+         * Skip the first 6 characters and eat the "lcl|x_".
+         */
+        else if (seqalign != NULL)
+          {
+            Char target_id_buf[BUFFER_LEN + 1];
+            SeqIdPtr target_id = SeqAlignId (seqalign, 1);
+            SeqIdWrite (target_id, target_id_buf, PRINTID_FASTA_SHORT, BUFFER_LEN);
+
+            // Species Type
+            dst_buf[i].type[0] = &target_id_buf[4];
+            dst_buf[i].type[1] = '\0';
+
+            // Protein Type
+            strncpy (dst_buf[i].protein, &target_id_buf[6], 
+                     sizeof (dst_buf[i].protein));
+
+            updates_pending = true;
+          }
+
+        /*
+         * BLAST did not find any hits.
+         */
+        else
+          {
+            printf ("Warning: Unable to identify protein type for sequence gi|%i\n",
+                    dst_buf[i].gi);
+          }
+        
+        /*
+         * Clean up memory for the next ieration.
+         */
+        seqalign = SeqAlignSetFree (seqalign);
+        bsp = BioseqFree (bsp);
+
+      }
+
+      /*
+       * Write the data out to the file.
+       */
+      if ( (i % 1000 == 0) && (i > 0) && updates_pending)
+        {
+          status = H5TBwrite_records (file_id, "influenza.faa", i - 1000, 1000,
+                                      dst_size, dst_offset, dst_sizes, 
+                                      &dst_buf[i-1000]);
+          if (status < 0)
+            check_h5_error (status, __FILE__, __LINE__);
+          
+          status = H5Fflush (file_id, H5F_SCOPE_GLOBAL);
+          if (status < 0)
+            check_h5_error (status, __FILE__, __LINE__);
+          
+          updates_pending = false;
+
+          printf ("Processed %i of %i records.\n", i, (int)nrecords);
+        }
+      
+    }
+
+  /*
+   * Write out records from the last bin if it was less than 1000
+   * records in size.
+   */
+  if (updates_pending)
+    {
+      if ((int)nrecords < 1000) 
+        {
+          status = H5TBwrite_records (file_id, "influenza.faa", 0, nrecords,
+                                      dst_size, dst_offset, dst_sizes,
+                                      dst_buf);
+        }
+      else
+        {
+          status = H5TBwrite_records (file_id, "influenza.faa", nrecords - 1000, 1000,
+                                      dst_size, dst_offset, dst_sizes, 
+                                      &dst_buf[nrecords-1000]);
+        }
+      if (status < 0)
+        check_h5_error (status, __FILE__, __LINE__);
+      
+      status = H5Fflush (file_id, H5F_SCOPE_GLOBAL);
+      if (status < 0)
+        check_h5_error (status, __FILE__, __LINE__);
+      
+      updates_pending = false;
+    }
+  
+  free (faa_buf);
+  free (gi_buf);
+  
+  options = BLASTOptionDelete (options);
+  
+  return;
+}