enzyme.c

/*********************************************************************************
 * MIT License                                                                   *
 *                                                                               *
 * Copyright (c) 2021 Chenxi Zhou <chnx.zhou@gmail.com>                          *
 *                                                                               *
 * Permission is hereby granted, free of charge, to any person obtaining a copy  *
 * of this software and associated documentation files (the "Software"), to deal *
 * in the Software without restriction, including without limitation the rights  *
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell     *
 * copies of the Software, and to permit persons to whom the Software is         *
 * furnished to do so, subject to the following conditions:                      *
 *                                                                               *
 * The above copyright notice and this permission notice shall be included in    *
 * all copies or substantial portions of the Software.                           *
 *                                                                               *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR    *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,      *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE   *
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER        *
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, *
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE *
 * SOFTWARE.                                                                     *
 *********************************************************************************/

/********************************** Revision History *****************************
 *                                                                               *
 * 14/12/21 - Chenxi Zhou: Created                                               *
 *                                                                               *
 *********************************************************************************/
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <ctype.h>
#include <float.h>

#include "kvec.h"
#include "enzyme.h"
#include "sdict.h"
#include "asset.h"

#undef DEBUG_ENZ

static double MIN_RE_DENS = .1;
static double MAX_RE_DENS = DBL_MAX;

re_cuts_t *re_cuts_init(uint32_t n)
{
    re_cuts_t *re_cuts = (re_cuts_t *) malloc(sizeof(re_cuts_t));
    re_cuts->n = n;
    re_cuts->density = 0;
    re_cuts->re = (re_t *) malloc(n * sizeof(re_t));
    uint32_t i;
    for (i = 0; i < n; ++i)
        re_cuts->re[i].sites = 0;
    return re_cuts;
}

void re_cuts_destroy(re_cuts_t *re_cuts)
{
    uint32_t i;
    for (i = 0; i < re_cuts->n; ++i)
        if (re_cuts->re[i].sites)
            free(re_cuts->re[i].sites);
    free(re_cuts->re);
    free(re_cuts);
}

typedef struct {size_t n, m; uint32_t *a;} u32_v;
int u32_cmp(const void *p, const void *q)
{
    uint32_t a, b;
    a = *(uint32_t *) p;
    b = *(uint32_t *) q;
    return a == b? 0 : (a > b? 1 : -1);
}

re_cuts_t *find_re_from_seqs(const char *f, uint32_t ml, char **enz_cs, int enz_n)
{
    // now find all RE cutting sites
    int i, j, c, c0, c1;
    uint32_t len, n, p;
    int64_t n_re, genome_size;
    char *pch, *seq, *re;
    sdict_t *sdict;
    re_cuts_t *re_cuts;

    sdict = make_sdict_from_fa(f, ml);
    n = sdict->n;
    n_re = genome_size = 0;
    re_cuts = re_cuts_init(n);
    for (i = 0; i < n; ++i) {
        u32_v enz_cs_pos = {0, 0, 0};

        seq = sdict->s[i].seq;
        len = sdict->s[i].len;
        for (p = 0; p < len; ++p) {
            c = seq[p];
            if (!isalpha(c)) {
                fprintf(stderr, "[E::%s] non-alphabetic chacrater in FASTA file: %c\n", __func__, c);
                re_cuts_destroy(re_cuts);
                sd_destroy(sdict);
                return 0;
            }
            seq[p] = nucl_toupper[c];
        }
        
        for (j = 0; j < enz_n; ++j) {
            re = enz_cs[j];
            pch = strstr(seq, re);
            while (pch != NULL) {
                kv_push(uint32_t, enz_cs_pos, pch - seq);
                pch = strstr(pch + 1, re);
            }
        }

        // reverse complement sequence
        for (p = 0; p < len>>1; ++p) {
            c0 = comp_table[(int) seq[p]];
            c1 = comp_table[(int) seq[len - 1 - p]];
            seq[p] = c1;
            seq[len - 1 - p] = c0;
        }
        if (len & 1) // complement the remaining base
            seq[len>>1] = comp_table[(int) seq[len>>1]];
        
        for (j = 0; j < enz_n; ++j) {
            re = enz_cs[j];
            pch = strstr(seq, re);
            while (pch != NULL) {
                kv_push(uint32_t, enz_cs_pos, len - 1 - (pch - seq));
                pch = strstr(pch + 1, re);
            }
        }

        qsort(enz_cs_pos.a, enz_cs_pos.n, sizeof(uint32_t), u32_cmp);
        re_cuts->re[i].sites = enz_cs_pos.a;
        re_cuts->re[i].n = enz_cs_pos.n;
        re_cuts->re[i].l = len;
        
        n_re += enz_cs_pos.n;
        genome_size += len;
    }

    re_cuts->density = (double) n_re / genome_size;

    fprintf(stderr, "[I::%s] NO. restriction enzyme cutting sites found in sequences: %ld\n", __func__, n_re);
    fprintf(stderr, "[I::%s] restriction enzyme cutting sites density: %.6f\n", __func__, re_cuts->density);
#ifdef DEBUG_ENZ
    printf("[I::%s] restriction enzyme cutting sites for individual sequences (n = %d)\n", __func__, n);
    for (i = 0; i < n; ++i)
        printf("[I::%s] %s %u %u %.6f\n", __func__, sdict->s[i].name, sdict->s[i].len, re_cuts->re[i].n, (double) re_cuts->re[i].n / sdict->s[i].len);
#endif

    sd_destroy(sdict);

    return re_cuts;
}

double **calc_re_cuts_density(re_cuts_t *re_cuts, uint32_t resolution)
{
    if (!re_cuts)
        return 0;

    uint32_t i, j, b;
    double **dens, *ds;
    re_t re;
    dens = (double **) malloc(re_cuts->n * sizeof(double *));
    for (i = 0; i < re_cuts->n; ++i) {
        re = re_cuts->re[i];
        b = div_ceil(re.l, resolution);
        ds = (double *) calloc(b, sizeof(double));
        for (j = 0; j < re.n; ++j)
            ds[(MAX(re.sites[j], 1) - 1) / resolution] += 1.;
        for (j = 0; j < b - 1; ++j)
            ds[j] /= (double) resolution * re_cuts->density;
        ds[b - 1] /= ((double) re.l - (double) (b - 1) * resolution) * re_cuts->density;
        
        for (j = 0; j < b; ++j)
            if (ds[j] < MIN_RE_DENS || ds[j] > MAX_RE_DENS)
                ds[j] = .0;
        
        dens[i] = ds;
#ifdef DEBUG_ENZ
        printf("DENS [%u/%u] (%u):", i, re_cuts->n, b);
        for (j = 0; j < b; ++j)
            printf(" %.6f", ds[j]);
        printf("\n");
#endif
    }
    return dens;
}

static uint32_t bin_search(uint32_t *a, uint32_t n, uint32_t s)
{
    uint32_t low, high, mid;
    low = 0;
    high = n;
    while (low != high) {
        mid = (low >> 1) + (high >> 1);
        if (a[mid] < s)
            low = mid + 1;
        else
            high = mid;
    }
    return low;
} 

double **calc_re_cuts_density1(re_cuts_t *re_cuts, uint32_t resolution, asm_dict_t *dict)
{
    if (!re_cuts)
        return 0;

    uint32_t i, j, b, n, e, a;
    double **dens, *ds;
    re_t re;
    sd_aseq_t seq;
    sd_seg_t seg;
    
    n = dict->n;
    dens = (double **) malloc(n * sizeof(double *));
    for (i = 0; i < n; ++i) {
        seq = dict->s[i];
        b = div_ceil(seq.len, resolution);
        ds = (double *) calloc(b, sizeof(double));
        for (j = 0; j < seq.n; ++j) {
            seg = dict->seg[seq.s + j];
            re = re_cuts->re[seg.c >> 1]; // get seq re cuts
            e = seg.x + seg.y; // seq end position, exclusive
            // binary search to get starting re cuts
            // first cutting site no smaller than seg.x
            // not many seq breaks - linear search might be faster on average?
            a = bin_search(re.sites, re.n, seg.x);
            
            if (seg.c & 1) {
                // reverse complement
                while (a < re.n && re.sites[a] < e) {
                    ds[(MAX(seg.a + seg.y - (re.sites[a] - seg.x), 1) - 1) / resolution] += 1.;
                    ++a;
                }
            } else {
                while (a < re.n && re.sites[a] < e) {
                    ds[(MAX(seg.a + re.sites[a] - seg.x, 1) - 1) / resolution] += 1.;
                    ++a;
                }
            }
        }
        for (j = 0; j < b - 1; ++j)
            ds[j] /= (double) resolution * re_cuts->density;
        ds[b - 1] /= ((double) seq.len - (double) (b - 1) * resolution) * re_cuts->density;
        
        for (j = 0; j < b; ++j)
            if (ds[j] < MIN_RE_DENS || ds[j] > MAX_RE_DENS)
                ds[j] = .0;
        
        dens[i] = ds;
#ifdef DEBUG_ENZ
        printf("DENS1 [%u/%u] (%u):", i, n, b);
        for (j = 0; j < b; ++j)
            printf(" %.6f", ds[j]);
        printf("\n");
#endif
    }

    return dens;
}

double **calc_re_cuts_density2(re_cuts_t *re_cuts, uint32_t resolution, asm_dict_t *dict)
{
    if (!re_cuts)
        return 0;

    uint32_t i, j, b, n, e, a;
    uint64_t l, p;
    double **dens, *ds;
    re_t re;
    sd_aseq_t seq;
    sd_seg_t seg;
    n = dict->n;
    dens = (double **) malloc(n * sizeof(double *));
    for (i = 0; i < n; ++i) {
        seq = dict->s[i];
        l = seq.len >> 1; // split sequence into two parts
        b = div_ceil(l, resolution);
        ds = (double *) calloc(b << 1, sizeof(double));
        for (j = 0; j < seq.n; ++j) {
            seg = dict->seg[seq.s + j];
            re = re_cuts->re[seg.c >> 1]; // get seq re cuts
            e = seg.x + seg.y; // seq end position, exclusive
            // binary search to get starting re cuts
            // first cutting site no smaller than seg.x
            // not many seq breaks - linear search might be faster on average?
            a = bin_search(re.sites, re.n, seg.x);
            if (seg.c & 1) {
                // reverse complement
                while (a < re.n && re.sites[a] < e) {
                    p = seg.a + seg.y - (re.sites[a] - seg.x); // position on seq
                    if (p < l)
                        ds[(MAX(p, 1) - 1) / resolution << 1] += 1.;
                    else
                        ds[(MAX(seq.len - p, 1) - 1) / resolution << 1 | 1] += 1.;
                    ++a;
                }
            } else {
                while (a < re.n && re.sites[a] < e) {
                    p = seg.a + re.sites[a] - seg.x; // position on seq
                    if (p < l)
                        ds[(MAX(p, 1) - 1) / resolution << 1] += 1.;
                    else
                        ds[(MAX(seq.len - p, 1) - 1) / resolution << 1 | 1] += 1.;
                    ++a;
                }
            }
        }
        for (j = 0; j < (b - 1) << 1; ++j)
            ds[j] /= (double) resolution * re_cuts->density;
        ds[(b - 1) << 1] /= ((double) l - (double) (b - 1) * resolution) * re_cuts->density;
        ds[(b - 1) << 1 | 1] /= ((double) l - (double) (b - 1) * resolution) * re_cuts->density;
        
        for (j = 0; j < b << 1; ++j)
            if (ds[j] < MIN_RE_DENS || ds[j] > MAX_RE_DENS)
                 ds[j] = .0;

        dens[i] = ds;
#ifdef DEBUG_ENZ
        printf("DENS2 [%u/%u] (%u):", i, n, b);
        for (j = 0; j < b << 1; ++j)
            printf(" %.6f", ds[j]);
        printf("\n");
#endif

    }

    return dens;
}