deps/kb_meme/meme-5.0.1/src/mhmm.c
void load_filter_process_motifs_for_hmm
(
  int request_n, // IN select the top n positive strand motifs
  RBTREE_T *request_nums, // IN select the requested motif numbers
  RBTREE_T *request_ids, // IN select the requested motif ids
  double request_ev, // IN select motifs with a <= evalue
  double request_complexity, // IN select motifs with a >= complexity
  const char *request_order, // IN specifiy motifs and their ordering and spacing
  double hit_pv, // IN specify the p value for keeping motif occurences
  const char *motif_filename, // IN file containing motifs to load
  const char *bg_source, // IN source of background
  double motif_pseudo, // IN pseudocount applied to motifs
  double trans_pseudo,  // IN pseudocount applied to transitions
  double spacer_pseudo, // IN pseudocount applied to spacers
  bool keep_unused, // IN keep the unused motifs
  ALPH_T **alphabet, // OUT alphabet of the motifs
  RBTREE_T **motifs, // OUT motifs with file position as key
  ARRAY_T **background, // OUT background from motif file
  ORDER_T **order_spacing, // OUT order and spacing (uses motif key as identifiers)
  MATRIX_T **transp_freq, // OUT transition frequencies
  MATRIX_T **spacer_ave // OUT spacer averages
) {
  int (*filter)(void *context, int motif_num, MOTIF_T *motif);
  void *filter_context;
  char *filter_name;
  REQMOT_T requested;
  RBTREE_T *ids2nums; // mapping of motif IDs to numbers
  RBNODE_T *node;
  MREAD_T *mread;
  MOTIF_T *motif;
  int motif_i, rc_motif_i, strands, i;
  bool created, has_reverse_strand;
  ARRAYLST_T *motif_occurrences;

  /**********************************************
   * SELECT MOTIF FILTER
   **********************************************/
  filter = NULL;
  filter_context = NULL;
  filter_name = "<no filter>";
  // filter top n
  if (request_n != 0) {
    filter = select_top_n;
    filter_context = &request_n;
    filter_name = "-nmotifs";
  }
  // filter requested
  if (rbtree_size(request_nums) != 0 || rbtree_size(request_ids) != 0) {
    if (filter != NULL) {
      die("Can't use -motif or -motif-id with %s.\n", filter_name);
    }
    requested.nums = request_nums;
    requested.ids = request_ids;
    filter = select_requested;
    filter_context = &requested;
    filter_name = "-motif or -motif-id";
  }
  // filter e threshold
  if (request_ev != 0.0) {
    if (filter != NULL) {
      die("Can't use -ethresh with %s.\n", filter_name);
    }
    filter = select_better_e_value;
    filter_context = &request_ev;
  }
  // filter complexity threshold
  if (request_complexity != 0.0) {
    if (filter != NULL) {
      die("Can't use -lowcomp with %s.\n", filter_name);
    }
    filter = select_more_complex;
    filter_context = &request_complexity;
  }
  // filter order (but not when reading motifs)
  if (request_order != NULL && filter != NULL) {
    die("Can't use -order with %s.\n", filter_name);
  }
  if (filter == NULL) filter = select_all;

  /**********************************************
   * READING THE MOTIFS
   **********************************************/

  mread = mread_create(motif_filename, OPEN_MFILE | CALC_AMBIGS | SCANNED_SITES, true);
  mread_set_bg_source(mread, bg_source, NULL);
  mread_set_pseudocount(mread, motif_pseudo);

  has_reverse_strand = (mread_get_strands(mread) == 2);

  *motifs = rbtree_create(motif_num_cmp, rbtree_intcpy, free, NULL, destroy_motif);
  ids2nums = rbtree_create(rbtree_strcmp, rbtree_strcpy, free, rbtree_intcpy, free);
  for (motif_i = 1, rc_motif_i = -1; mread_has_motif(mread); ++motif_i, --rc_motif_i) {
    motif = mread_next_motif(mread);
    strands = filter(filter_context, motif_i, motif);
    if (strands & MOTIF_PLUS) rbtree_make(*motifs, &motif_i, motif);
    if (strands & MOTIF_MINUS) rbtree_make(*motifs, &rc_motif_i, dup_rc_motif(motif));
    if (strands) rbtree_make(ids2nums, get_motif_id(motif), &motif_i);
    else destroy_motif(motif);
  }

  *background = mread_get_background(mread);
  *alphabet = alph_hold(mread_get_alphabet(mread));
  motif_occurrences = mread_get_motif_occurrences(mread);

  mread_destroy(mread);

  // Parse the order string. 
  *order_spacing = create_order_from_ids(request_order, ids2nums);
  rbtree_destroy(ids2nums);

  if (*order_spacing) {
    // Remove unused motifs
    for (node = rbtree_first(*motifs); node != NULL; node = rbtree_next(node)) {
      motif_i = *((int*)rbtree_key(node));
      if (!order_contains(motif_i, *order_spacing)) {
        rbtree_delete(*motifs, node, NULL, NULL);
      }
    }
  }

  // Turn the raw motifs and motif occurrences into the elements of the model
  if (motif_occurrences != NULL && arraylst_size(motif_occurrences) > 0) {
    compute_transitions_and_spacers(
        *motifs,
        motif_occurrences,
        has_reverse_strand,
        hit_pv,
        order_spacing,
        transp_freq,
        spacer_ave);
  } else {
    // If no occurrences are found, initialize matrices uniformly.
    compute_naive_transitions_and_spacers(
        rbtree_size(*motifs), 
        transp_freq, 
        spacer_ave);
  }

  // Convert spacer info to probabilities.
  normalize_spacer_counts(
      trans_pseudo, 
      spacer_pseudo,
      keep_unused,
      *transp_freq, 
      *spacer_ave);

  // Throw out unused motifs.
  throw_out_unused_motifs(*motifs, *transp_freq, *spacer_ave);

  destroy_motif_occurrences(motif_occurrences);
}
/***********************************************************************
 * Read the motifs, filter them, create an order and spacing if
 * possible and create the transition frequency matrix and the
 * spacer average matrix.
 ***********************************************************************/
Read the motifs, filter them, create an order and spacing if
possible and create the transition frequency matrix and the
spacer average matrix.
Read the motifs filter them create an order and spacing if possible and create the transition frequency matrix and the spacer average matrix .
void load_filter_process_motifs_for_hmm
(
  int request_n, 
  RBTREE_T *request_nums, 
  RBTREE_T *request_ids, 
  double request_ev, 
  double request_complexity, 
  const char *request_order, 
  double hit_pv, 
  const char *motif_filename, 
  const char *bg_source, 
  double motif_pseudo, 
  double trans_pseudo,  
  double spacer_pseudo, 
  bool keep_unused, 
  ALPH_T **alphabet, 
  RBTREE_T **motifs, 
  ARRAY_T **background, 
  ORDER_T **order_spacing, 
  MATRIX_T **transp_freq, 
  MATRIX_T **spacer_ave 
) {
  int (*filter)(void *context, int motif_num, MOTIF_T *motif);
  void *filter_context;
  char *filter_name;
  REQMOT_T requested;
  RBTREE_T *ids2nums; 
  RBNODE_T *node;
  MREAD_T *mread;
  MOTIF_T *motif;
  int motif_i, rc_motif_i, strands, i;
  bool created, has_reverse_strand;
  ARRAYLST_T *motif_occurrences;
  filter = NULL;
  filter_context = NULL;
  filter_name = "<no filter>";
  if (request_n != 0) {
    filter = select_top_n;
    filter_context = &request_n;
    filter_name = "-nmotifs";
  }
  if (rbtree_size(request_nums) != 0 || rbtree_size(request_ids) != 0) {
    if (filter != NULL) {
      die("Can't use -motif or -motif-id with %s.\n", filter_name);
    }
    requested.nums = request_nums;
    requested.ids = request_ids;
    filter = select_requested;
    filter_context = &requested;
    filter_name = "-motif or -motif-id";
  }
  if (request_ev != 0.0) {
    if (filter != NULL) {
      die("Can't use -ethresh with %s.\n", filter_name);
    }
    filter = select_better_e_value;
    filter_context = &request_ev;
  }
  if (request_complexity != 0.0) {
    if (filter != NULL) {
      die("Can't use -lowcomp with %s.\n", filter_name);
    }
    filter = select_more_complex;
    filter_context = &request_complexity;
  }
  if (request_order != NULL && filter != NULL) {
    die("Can't use -order with %s.\n", filter_name);
  }
  if (filter == NULL) filter = select_all;
  mread = mread_create(motif_filename, OPEN_MFILE | CALC_AMBIGS | SCANNED_SITES, true);
  mread_set_bg_source(mread, bg_source, NULL);
  mread_set_pseudocount(mread, motif_pseudo);
  has_reverse_strand = (mread_get_strands(mread) == 2);
  *motifs = rbtree_create(motif_num_cmp, rbtree_intcpy, free, NULL, destroy_motif);
  ids2nums = rbtree_create(rbtree_strcmp, rbtree_strcpy, free, rbtree_intcpy, free);
  for (motif_i = 1, rc_motif_i = -1; mread_has_motif(mread); ++motif_i, --rc_motif_i) {
    motif = mread_next_motif(mread);
    strands = filter(filter_context, motif_i, motif);
    if (strands & MOTIF_PLUS) rbtree_make(*motifs, &motif_i, motif);
    if (strands & MOTIF_MINUS) rbtree_make(*motifs, &rc_motif_i, dup_rc_motif(motif));
    if (strands) rbtree_make(ids2nums, get_motif_id(motif), &motif_i);
    else destroy_motif(motif);
  }
  *background = mread_get_background(mread);
  *alphabet = alph_hold(mread_get_alphabet(mread));
  motif_occurrences = mread_get_motif_occurrences(mread);
  mread_destroy(mread);
  *order_spacing = create_order_from_ids(request_order, ids2nums);
  rbtree_destroy(ids2nums);
  if (*order_spacing) {
    for (node = rbtree_first(*motifs); node != NULL; node = rbtree_next(node)) {
      motif_i = *((int*)rbtree_key(node));
      if (!order_contains(motif_i, *order_spacing)) {
        rbtree_delete(*motifs, node, NULL, NULL);
      }
    }
  }
  if (motif_occurrences != NULL && arraylst_size(motif_occurrences) > 0) {
    compute_transitions_and_spacers(
        *motifs,
        motif_occurrences,
        has_reverse_strand,
        hit_pv,
        order_spacing,
        transp_freq,
        spacer_ave);
  } else {
    compute_naive_transitions_and_spacers(
        rbtree_size(*motifs), 
        transp_freq, 
        spacer_ave);
  }
  normalize_spacer_counts(
      trans_pseudo, 
      spacer_pseudo,
      keep_unused,
      *transp_freq, 
      *spacer_ave);
  throw_out_unused_motifs(*motifs, *transp_freq, *spacer_ave);
  destroy_motif_occurrences(motif_occurrences);
}
void load_filter_process_motifs_for_hmm ( int request_n , RBTREE_T * request_nums , RBTREE_T * request_ids , double request_ev , double request_complexity , const char * request_order , double hit_pv , const char * motif_filename , const char * bg_source , double motif_pseudo , double trans_pseudo , double spacer_pseudo , bool keep_unused , ALPH_T * * alphabet , RBTREE_T * * motifs , ARRAY_T * * background , ORDER_T * * order_spacing , MATRIX_T * * transp_freq , MATRIX_T * * spacer_ave ) { int ( * filter ) ( void * context , int motif_num , MOTIF_T * motif ) ; void * filter_context ; char * filter_name ; REQMOT_T requested ; RBTREE_T * ids2nums ; RBNODE_T * node ; MREAD_T * mread ; MOTIF_T * motif ; int motif_i , rc_motif_i , strands , i ; bool created , has_reverse_strand ; ARRAYLST_T * motif_occurrences ; filter = NULL ; filter_context = NULL ; filter_name = " " ; if ( request_n != 0 ) { filter = select_top_n ; filter_context = & request_n ; filter_name = " " ; } if ( rbtree_size ( request_nums ) != 0 || rbtree_size ( request_ids ) != 0 ) { if ( filter != NULL ) { die ( " \n " , filter_name ) ; } requested . nums = request_nums ; requested . ids = request_ids ; filter = select_requested ; filter_context = & requested ; filter_name = " " ; } if ( request_ev != 0.0 ) { if ( filter != NULL ) { die ( " \n " , filter_name ) ; } filter = select_better_e_value ; filter_context = & request_ev ; } if ( request_complexity != 0.0 ) { if ( filter != NULL ) { die ( " \n " , filter_name ) ; } filter = select_more_complex ; filter_context = & request_complexity ; } if ( request_order != NULL && filter != NULL ) { die ( " \n " , filter_name ) ; } if ( filter == NULL ) filter = select_all ; mread = mread_create ( motif_filename , OPEN_MFILE | CALC_AMBIGS | SCANNED_SITES , true ) ; mread_set_bg_source ( mread , bg_source , NULL ) ; mread_set_pseudocount ( mread , motif_pseudo ) ; has_reverse_strand = ( mread_get_strands ( mread ) == 2 ) ; * motifs = rbtree_create ( motif_num_cmp , rbtree_intcpy , free , NULL , destroy_motif ) ; ids2nums = rbtree_create ( rbtree_strcmp , rbtree_strcpy , free , rbtree_intcpy , free ) ; for ( motif_i = 1 , rc_motif_i = -1 ; mread_has_motif ( mread ) ; ++ motif_i , -- rc_motif_i ) { motif = mread_next_motif ( mread ) ; strands = filter ( filter_context , motif_i , motif ) ; if ( strands & MOTIF_PLUS ) rbtree_make ( * motifs , & motif_i , motif ) ; if ( strands & MOTIF_MINUS ) rbtree_make ( * motifs , & rc_motif_i , dup_rc_motif ( motif ) ) ; if ( strands ) rbtree_make ( ids2nums , get_motif_id ( motif ) , & motif_i ) ; else destroy_motif ( motif ) ; } * background = mread_get_background ( mread ) ; * alphabet = alph_hold ( mread_get_alphabet ( mread ) ) ; motif_occurrences = mread_get_motif_occurrences ( mread ) ; mread_destroy ( mread ) ; * order_spacing = create_order_from_ids ( request_order , ids2nums ) ; rbtree_destroy ( ids2nums ) ; if ( * order_spacing ) { for ( node = rbtree_first ( * motifs ) ; node != NULL ; node = rbtree_next ( node ) ) { motif_i = * ( ( int * ) rbtree_key ( node ) ) ; if ( ! order_contains ( motif_i , * order_spacing ) ) { rbtree_delete ( * motifs , node , NULL , NULL ) ; } } } if ( motif_occurrences != NULL && arraylst_size ( motif_occurrences ) > 0 ) { compute_transitions_and_spacers ( * motifs , motif_occurrences , has_reverse_strand , hit_pv , order_spacing , transp_freq , spacer_ave ) ; } else { compute_naive_transitions_and_spacers ( rbtree_size ( * motifs ) , transp_freq , spacer_ave ) ; } normalize_spacer_counts ( trans_pseudo , spacer_pseudo , keep_unused , * transp_freq , * spacer_ave ) ; throw_out_unused_motifs ( * motifs , * transp_freq , * spacer_ave ) ; destroy_motif_occurrences ( motif_occurrences ) ; }
Read the motifs, filter them, create an order and spacing if
 possible and create the transition frequency matrix and the
 spacer average matrix.
Read the motifs filter them create an order and spacing if possible and create the transition frequency matrix and the spacer average matrix .
// IN select the top n positive strand motifs // IN select the requested motif numbers // IN select the requested motif ids // IN select motifs with a <= evalue // IN select motifs with a >= complexity // IN specifiy motifs and their ordering and spacing // IN specify the p value for keeping motif occurences // IN file containing motifs to load // IN source of background // IN pseudocount applied to motifs // IN pseudocount applied to transitions // IN pseudocount applied to spacers // IN keep the unused motifs // OUT alphabet of the motifs // OUT motifs with file position as key // OUT background from motif file // OUT order and spacing (uses motif key as identifiers) // OUT transition frequencies // OUT spacer averages // mapping of motif IDs to numbers /**********************************************
   * SELECT MOTIF FILTER
   **********************************************/ // filter top n // filter requested // filter e threshold // filter complexity threshold // filter order (but not when reading motifs) /**********************************************
   * READING THE MOTIFS
   **********************************************/ // Parse the order string.  // Remove unused motifs // Turn the raw motifs and motif occurrences into the elements of the model // If no occurrences are found, initialize matrices uniformly. // Convert spacer info to probabilities. // Throw out unused motifs.