<filename>tad.py<gh_stars>0
# calling TAD using deTOKI
# migrate to notebook compatible 
import numpy as np
import pandas as pd
from sklearn import decomposition
import math
import multiprocessing

## Run NMF in several times with random initialisation and output consensus matrix
def corate(A,n,time):
    # Input:
    #   A : matrix, window-sized submatrix
    #   n : number of components, optimized by 
    #       {average_length_of_TAD, length_of_window}
    #   time : times of NMF running
    # Output:
    #   consensus NMF matrix
    S=np.zeros([np.shape(A)[0],np.shape(A)[1]])
    for i in range(time):
        K=np.zeros([np.shape(A)[0],np.shape(A)[1]])
        estimator=decomposition.NMF(n_components=n, init='random',random_state=i, max_iter=10_000)
        estimator.fit(A)
        B=estimator.transform(A)
        index=B.argmax(axis=1)
        for j in range(n):
            for s in np.where(index==j)[0]:
                for t in np.where(index==j)[0]:
                    K[s,t]=1
        S=S+K
    return S.astype(np.float64)/time
## Define silhouette coefficient of consensus matrix and detected TAD boundaries
def silhou(R, pos):
    n=np.shape(R)[0]
    silhou=0
    for i in range(len(pos)-1):
        for j in range(pos[i],pos[i+1]):
            a=np.sum((1-R)[j,pos[i]:pos[i+1]])
            b=np.sum((1-R)[i,:])-a
            silhou+=(-a/(pos[i+1]-pos[i])+b/(n+pos[i]-pos[i+1]))/max((a/(pos[i+1]-pos[i]),b/(n+pos[i]-pos[i+1])))
    return silhou/n

## Calculate clustering rate of each bin
def IS(R, length):
    # Calculate clustering rate(number of contacts 
    # locatingin  inter_domain triangle field of the bin)

    # Input:
    #   R: contact matrix, typically consensus_NMF_matrix
    #   length : length of bin's flanking regions considered
    bias=np.zeros([np.shape(R)[0]])
    for i in range(1,np.shape(R)[0]-1):
        bias[i]=np.mean(
            R[max(0,i-length):(i+1), i:min(i+length+1,np.shape(R)[0])]
            )
    return bias

## Find the best n_components by comparing silhouette coefficient
def bestco(F, reso, length, delta, min_val, max_val):
    # argument picking
    # find best n(number of components)
    # judged by silhouette coefficient

    # n ~ number of TADs in 1 window

    # rely on:
    #   corate, silhou, zero
    # Input:
    #   F, reso, min_val, max_val: for n determining
    #   length, delta: for boundary calling
    # Output:
    #   final consensus NMF matrix, and its n
    x=-1
    R0=0
    n1=0
    for n in range(max(int(np.shape(F)[0]*reso/max_val),1),int(np.shape(F)[0]*reso/min_val)+1):
        # get NMF_consensus_matrix
        R1=corate(F,n,10)
        # calc tad under above n/NMF_consensus_matrix
        x1=silhou(R1, zero(R1, n-1, length, delta))
        if x1>=x:
            R0=R1
            x=x1
            n1=n
    # final R0, best consensus NMF matrix
    # fianl n1, best n_component used in NMF
    return R0,n1

## Find bins with local minimal clustering rate and global comparative low clustering rate （these bins are detected TAD boundaries)
def zero(R, t, length, delta):
    # pick at most t TAD boundary from contact matrix
    # Input:
    #   R : small piece of matrix
    #   t : number of output boundaries

    # find local valley of IS result

    #bias: number of contacts in "anti-diag" field
    bias=IS(R, length)
    #delta_list: right_bias - left_bias
    delta_list=[]
    for i in range(delta,len(bias)-delta):
        delta_list.append(-sum(bias[(i-delta):i])+sum(bias[i:(i+delta)]))
    # zero: position of local valleys
    zero=[0]
    for i in range(len(delta_list)-1):
        if delta_list[i]<0 and delta_list[i+1]>=0:
            #zero.append(i+5-np.argmin([bias[i+x] for x in range(5,0,-1)]))
            zero.append(i+delta) # why "+delta"?
    zero.append(len(bias))
    zero=sorted(np.unique(zero))

    # pick boundaries near strong TAD

    # enrich: subset of zero, picked by index
    enrich=[0]
    # strength: max summit value upstream and downstream the valley
    #   aka inner tad strength of tad's alongside of the boundary(zero pos)
    strength=[]
    # dicard left and right end
    for j in range(1,len(zero)-1):
        strength.append(
            max(
                # upstream summit
                max( bias[zero[j-1]:(zero[j]+1)] ) - bias[zero[j]],
                # downstream summit
                max( bias[(zero[j]-1):zero[j+1]] ) - bias[zero[j]]
                )
            )
    strength=np.array(strength)
    # select globally from zeros, using strength
    # top t, and value>0.3
    index = sorted(list(set(np.argsort(-strength)[:t])|set(np.where(strength>0.3)[0])))
    # store result in enrich, adding left and right end
    for k in index:
        enrich.append(zero[k+1])
    enrich.append(len(bias))
    return np.array(enrich)
## Split huge contact matrix to windows and detected TAD boundaries in each window
def part_zero(F:np.ndarray, core:int, window:float, reso:int, delta, length, min_val, max_val):
    # depend on : bestco, zero
    
    #Input:
    #   F: full contact_matrix
    #   core: number of threads
    #   window, reso, length, delta -> bestco
    #   length, delta -> zero
    
    pos=[]
    # length of chromosome/input_matrix, in num of reso_bins
    n=np.shape(F)[0]

    # deTOKI's core function

    global task
    def task(i):
        # i : work id
        # (implicit)F :  full mat
        # (implicit)window : window size (uint: reso_bins)
        # (implicit)n : length of full chromosome 

        # P : a window-size-tile of F, different near chromosome ends
        P=F[max(0,window//2*i-window//4):min(n,window//2*i+window-window//4),
            max(0,window//2*i-window//4):min(n,window//2*i+window-window//4)]
        if np.sum(P)<100:
            # poor signal
            # windows with less than 100 contacts yield 0 tad boundaries
            return []
        # find best NMF_consensus_matrix
        R, t = bestco(P, reso, length, delta, min_val, max_val)
        if t==0:
            # all n candidates are not good
            return []
        # call TAD on consensus matrix
        p=zero(R, t-1, length, delta)
        pos=[]
        # transform TAD boundary from local(in window) index to global index
        for j in p:
            if j in range(window//2*i-max(0,window//2*i-window//4),window//2*(i+1)-max(0,window//2*i-window//4)):
                pos.append(j+max(0,window//2*i-window//4))
        return pos

    # multi-thread calc
    
    pool = multiprocessing.Pool(processes=core)
    res_list=[]
    for i in range(math.ceil(2*n/window)):
        # iterate by window's mid-points
        # window size is 8mb,
        # each window overlaps neighbor windows by 4mb,
        # yield 2*n/window mid-points
        result = pool.apply_async(task, args=(i,))
        res_list.append(result)
    pool.close()
    pool.join()

    # async get and join every window's tad boundary

    for i in range(math.ceil(2*n/window)):
        result=res_list[i]
        pos=np.append(pos,result.get())
    return pos.astype('int32')
## entry function
def tad(F, core:int=1, reso:int=40, min_val:int=600, max_val:int=1000, window_size_raw:int=8000, delta_raw:int=100):
#Input:
#   F: input contact matrix : np.ndarray
#   reso: resolution of the matrix (kbp)
#   min, max: tad mean size (kbp)
#   split: window size (kbp) 
#   core: treads used

    # length of flanking regions considered
    # by default 400kbp
    length=400//reso
    # small intervals for calculating curve's summints and valleys
    delta=int(math.ceil(delta_raw/reso))
    # tiling windows' size(number of res_bins), for parallel calculation 
    window_size = window_size_raw//reso

    return part_zero(F, core, window, reso, delta, length, min_val, max_val)
<filename>tad.py<gh_stars>0
# calling TAD using deTOKI
# migrate to notebook compatible 
import numpy as np
import pandas as pd
from sklearn import decomposition
import math
import multiprocessing

## Run NMF in several times with random initialisation and output consensus matrix
def corate(A,n,time):
    # Input:
    #   A : matrix, window-sized submatrix
    #   n : number of components, optimized by 
    #       {average_length_of_TAD, length_of_window}
    #   time : times of NMF running
    # Output:
    #   consensus NMF matrix
    S=np.zeros([np.shape(A)[0],np.shape(A)[1]])
    for i in range(time):
        K=np.zeros([np.shape(A)[0],np.shape(A)[1]])
        estimator=decomposition.NMF(n_components=n, init='random',random_state=i, max_iter=10_000)
        estimator.fit(A)
        B=estimator.transform(A)
        index=B.argmax(axis=1)
        for j in range(n):
            for s in np.where(index==j)[0]:
                for t in np.where(index==j)[0]:
                    K[s,t]=1
        S=S+K
    return S.astype(np.float64)/time
## Define silhouette coefficient of consensus matrix and detected TAD boundaries
def silhou(R, pos):
    n=np.shape(R)[0]
    silhou=0
    for i in range(len(pos)-1):
        for j in range(pos[i],pos[i+1]):
            a=np.sum((1-R)[j,pos[i]:pos[i+1]])
            b=np.sum((1-R)[i,:])-a
            silhou+=(-a/(pos[i+1]-pos[i])+b/(n+pos[i]-pos[i+1]))/max((a/(pos[i+1]-pos[i]),b/(n+pos[i]-pos[i+1])))
    return silhou/n

## Calculate clustering rate of each bin
def IS(R, length):
    # Calculate clustering rate(number of contacts 
    # locatingin  inter_domain triangle field of the bin)

    # Input:
    #   R: contact matrix, typically consensus_NMF_matrix
    #   length : length of bin's flanking regions considered
    bias=np.zeros([np.shape(R)[0]])
    for i in range(1,np.shape(R)[0]-1):
        bias[i]=np.mean(
            R[max(0,i-length):(i+1), i:min(i+length+1,np.shape(R)[0])]
            )
    return bias

## Find the best n_components by comparing silhouette coefficient
def bestco(F, reso, length, delta, min_val, max_val):
    # argument picking
    # find best n(number of components)
    # judged by silhouette coefficient

    # n ~ number of TADs in 1 window

    # rely on:
    #   corate, silhou, zero
    # Input:
    #   F, reso, min_val, max_val: for n determining
    #   length, delta: for boundary calling
    # Output:
    #   final consensus NMF matrix, and its n
    x=-1
    R0=0
    n1=0
    for n in range(max(int(np.shape(F)[0]*reso/max_val),1),int(np.shape(F)[0]*reso/min_val)+1):
        # get NMF_consensus_matrix
        R1=corate(F,n,10)
        # calc tad under above n/NMF_consensus_matrix
        x1=silhou(R1, zero(R1, n-1, length, delta))
        if x1>=x:
            R0=R1
            x=x1
            n1=n
    # final R0, best consensus NMF matrix
    # fianl n1, best n_component used in NMF
    return R0,n1

## Find bins with local minimal clustering rate and global comparative low clustering rate （these bins are detected TAD boundaries)
def zero(R, t, length, delta):
    # pick at most t TAD boundary from contact matrix
    # Input:
    #   R : small piece of matrix
    #   t : number of output boundaries

    # find local valley of IS result

    #bias: number of contacts in "anti-diag" field
    bias=IS(R, length)
    #delta_list: right_bias - left_bias
    delta_list=[]
    for i in range(delta,len(bias)-delta):
        delta_list.append(-sum(bias[(i-delta):i])+sum(bias[i:(i+delta)]))
    # zero: position of local valleys
    zero=[0]
    for i in range(len(delta_list)-1):
        if delta_list[i]<0 and delta_list[i+1]>=0:
            #zero.append(i+5-np.argmin([bias[i+x] for x in range(5,0,-1)]))
            zero.append(i+delta) # why "+delta"?
    zero.append(len(bias))
    zero=sorted(np.unique(zero))

    # pick boundaries near strong TAD

    # enrich: subset of zero, picked by index
    enrich=[0]
    # strength: max summit value upstream and downstream the valley
    #   aka inner tad strength of tad's alongside of the boundary(zero pos)
    strength=[]
    # dicard left and right end
    for j in range(1,len(zero)-1):
        strength.append(
            max(
                # upstream summit
                max( bias[zero[j-1]:(zero[j]+1)] ) - bias[zero[j]],
                # downstream summit
                max( bias[(zero[j]-1):zero[j+1]] ) - bias[zero[j]]
                )
            )
    strength=np.array(strength)
    # select globally from zeros, using strength
    # top t, and value>0.3
    index = sorted(list(set(np.argsort(-strength)[:t])|set(np.where(strength>0.3)[0])))
    # store result in enrich, adding left and right end
    for k in index:
        enrich.append(zero[k+1])
    enrich.append(len(bias))
    return np.array(enrich)
## Split huge contact matrix to windows and detected TAD boundaries in each window
def part_zero(F:np.ndarray, core:int, window:float, reso:int, delta, length, min_val, max_val):
    # depend on : bestco, zero
    
    #Input:
    #   F: full contact_matrix
    #   core: number of threads
    #   window, reso, length, delta -> bestco
    #   length, delta -> zero
    
    pos=[]
    # length of chromosome/input_matrix, in num of reso_bins
    n=np.shape(F)[0]

    # deTOKI's core function

    global task
    def task(i):
        # i : work id
        # (implicit)F :  full mat
        # (implicit)window : window size (uint: reso_bins)
        # (implicit)n : length of full chromosome 

        # P : a window-size-tile of F, different near chromosome ends
        P=F[max(0,window//2*i-window//4):min(n,window//2*i+window-window//4),
            max(0,window//2*i-window//4):min(n,window//2*i+window-window//4)]
        if np.sum(P)<100:
            # poor signal
            # windows with less than 100 contacts yield 0 tad boundaries
            return []
        # find best NMF_consensus_matrix
        R, t = bestco(P, reso, length, delta, min_val, max_val)
        if t==0:
            # all n candidates are not good
            return []
        # call TAD on consensus matrix
        p=zero(R, t-1, length, delta)
        pos=[]
        # transform TAD boundary from local(in window) index to global index
        for j in p:
            if j in range(window//2*i-max(0,window//2*i-window//4),window//2*(i+1)-max(0,window//2*i-window//4)):
                pos.append(j+max(0,window//2*i-window//4))
        return pos

    # multi-thread calc
    
    pool = multiprocessing.Pool(processes=core)
    res_list=[]
    for i in range(math.ceil(2*n/window)):
        # iterate by window's mid-points
        # window size is 8mb,
        # each window overlaps neighbor windows by 4mb,
        # yield 2*n/window mid-points
        result = pool.apply_async(task, args=(i,))
        res_list.append(result)
    pool.close()
    pool.join()

    # async get and join every window's tad boundary

    for i in range(math.ceil(2*n/window)):
        result=res_list[i]
        pos=np.append(pos,result.get())
    return pos.astype('int32')
## entry function
def tad(F, core:int=1, reso:int=40, min_val:int=600, max_val:int=1000, window_size_raw:int=8000, delta_raw:int=100):
#Input:
#   F: input contact matrix : np.ndarray
#   reso: resolution of the matrix (kbp)
#   min, max: tad mean size (kbp)
#   split: window size (kbp) 
#   core: treads used

    # length of flanking regions considered
    # by default 400kbp
    length=400//reso
    # small intervals for calculating curve's summints and valleys
    delta=int(math.ceil(delta_raw/reso))
    # tiling windows' size(number of res_bins), for parallel calculation 
    window_size = window_size_raw//reso

    return part_zero(F, core, window, reso, delta, length, min_val, max_val)
# calling TAD using deTOKI # migrate to notebook compatible ## Run NMF in several times with random initialisation and output consensus matrix # Input: #   A : matrix, window-sized submatrix #   n : number of components, optimized by #       {average_length_of_TAD, length_of_window} #   time : times of NMF running # Output: #   consensus NMF matrix ## Define silhouette coefficient of consensus matrix and detected TAD boundaries ## Calculate clustering rate of each bin # Calculate clustering rate(number of contacts # locatingin  inter_domain triangle field of the bin) # Input: #   R: contact matrix, typically consensus_NMF_matrix #   length : length of bin's flanking regions considered ## Find the best n_components by comparing silhouette coefficient # argument picking # find best n(number of components) # judged by silhouette coefficient # n ~ number of TADs in 1 window # rely on: #   corate, silhou, zero # Input: #   F, reso, min_val, max_val: for n determining #   length, delta: for boundary calling # Output: #   final consensus NMF matrix, and its n # get NMF_consensus_matrix # calc tad under above n/NMF_consensus_matrix # final R0, best consensus NMF matrix # fianl n1, best n_component used in NMF ## Find bins with local minimal clustering rate and global comparative low clustering rate （these bins are detected TAD boundaries) # pick at most t TAD boundary from contact matrix # Input: #   R : small piece of matrix #   t : number of output boundaries # find local valley of IS result #bias: number of contacts in "anti-diag" field #delta_list: right_bias - left_bias # zero: position of local valleys #zero.append(i+5-np.argmin([bias[i+x] for x in range(5,0,-1)])) # why "+delta"? # pick boundaries near strong TAD # enrich: subset of zero, picked by index # strength: max summit value upstream and downstream the valley #   aka inner tad strength of tad's alongside of the boundary(zero pos) # dicard left and right end # upstream summit # downstream summit # select globally from zeros, using strength # top t, and value>0.3 # store result in enrich, adding left and right end ## Split huge contact matrix to windows and detected TAD boundaries in each window # depend on : bestco, zero #Input: #   F: full contact_matrix #   core: number of threads #   window, reso, length, delta -> bestco #   length, delta -> zero # length of chromosome/input_matrix, in num of reso_bins # deTOKI's core function # i : work id # (implicit)F :  full mat # (implicit)window : window size (uint: reso_bins) # (implicit)n : length of full chromosome # P : a window-size-tile of F, different near chromosome ends # poor signal # windows with less than 100 contacts yield 0 tad boundaries # find best NMF_consensus_matrix # all n candidates are not good # call TAD on consensus matrix # transform TAD boundary from local(in window) index to global index # multi-thread calc # iterate by window's mid-points # window size is 8mb, # each window overlaps neighbor windows by 4mb, # yield 2*n/window mid-points # async get and join every window's tad boundary ## entry function #Input: #   F: input contact matrix : np.ndarray #   reso: resolution of the matrix (kbp) #   min, max: tad mean size (kbp) #   split: window size (kbp) #   core: treads used # length of flanking regions considered # by default 400kbp # small intervals for calculating curve's summints and valleys # tiling windows' size(number of res_bins), for parallel calculation