#!/grid/gillis/home/nfox/software/miniconda3/bin/python #!/sonas-hs/gillis/hpc/home/nfox/software/miniconda3/bin/python import argparse import os import shutil import numpy as np import pandas as pd import numpy as np import scipy.spatial.distance as ssd import gc import h5py import itertools import scipy.stats as stats from scipy.io import mmread from scipy import sparse, io from scipy.sparse import csr_matrix, triu, coo_matrix import numpy as np #from hic_normalization import * import sys class HiCNorm(object): def __init__(self, matrix): self.bias = None self.matrix = matrix self.norm_matrix = None self._bias_var = None def iterative_correction(self, max_iter=400): #mat = self.matrix * 1.0 ## use float to avoid int overflow mat_lil = sparse.lil_matrix(self.matrix) # remove low count rows for normalization row_sum = np.asarray(np.sum(mat_lil, axis=1)).ravel() low_count = np.quantile(row_sum, 0.15) # 0.15 gives best correlation between KR/SP/VC mask_row = row_sum < low_count #mask_row = np.asarray(mask_row.ravel()) #mask_row = [count for count, value in enumerate(mask_row) if value] #print (mask_row) #mat_lil = sparse.lil_matrix(mat) mat_lil[mask_row, :] = 0 mat_lil[:, mask_row] = 0 self.bias = np.ones(mat_lil.shape[0]) print (self.bias.dtype) self._bias_var = [] x, y = np.nonzero(mat_lil) mat_sum = np.sum(mat_lil) # force the sum of matrix to stay the same, otherwise ICE may not converge as fast for i in range(max_iter): #bias =sparse.csr_matrix.sum(mat, axis=1) bias = np.asarray(np.sum(mat_lil, axis=1)).ravel() bias_mean = np.mean(bias[bias > 0]) bias_var = np.var(bias[bias > 0]) self._bias_var.append(bias_var) bias = bias / bias_mean bias[bias == 0] = 1 #print (np.asarray(bias)[:,:].ravel()) #print (bias[x]) #print (bias[x]*bias[y]) mat = sparse.csr_matrix(mat_lil) mat[x, y] = mat[x, y] / (bias[x]*bias[y]) #print (mat) new_sum = np.sum(mat) mat = mat * (mat_sum / new_sum) mat_lil = sparse.lil_matrix(mat) # update total bias self.bias = self.bias * bias * np.sqrt(new_sum / mat_sum) if bias_var < .001: break else: pass self.norm_matrix = csr_matrix(self.matrix, dtype=np.float64) self.norm_matrix[x, y] = self.norm_matrix[x, y] / (self.bias[x] * self.bias[y]) def vanilla_coverage(self): self.iterative_correction(max_iter=1) def make_tss_tss(gene_data, chrom, dist_folder, RAWobserved, resoulution_in_kb, ind_network): RAWobserved_csr = sparse.csr_matrix(RAWobserved) #converted coordiante format to csr format RAWobserved_csr.data = np.nan_to_num(RAWobserved_csr.data) raw_matrix_ice_sparse = sparse.csr_matrix.transpose(RAWobserved_csr) + sparse.triu(RAWobserved_csr, k=1, format='csr') row_sum = np.asarray(np.sum(raw_matrix_ice_sparse, axis=1)).ravel() map_resolution = np.quantile(row_sum, 0.20) maximum_read = (np.amax(raw_matrix_ice_sparse)) mat_ic = HiCNorm(raw_matrix_ice_sparse) mat_ic.iterative_correction() bias_var = mat_ic._bias_var maximum_bias = bias_var[-1] return maximum_read, map_resolution, maximum_bias def make_distance_matrix(gene_folder, dist_folder, resoulution_in_kb=100, ind_network=""): print (resoulution_in_kb) SRA_list = [] chr_list = [] resolution_list = [] maximum_read_list = [] map_resolution_list = [] maximum_bias_list = [] if os.path.isdir(dist_folder): # shutil.rmtree(dist_folder) pass else: os.mkdir(dist_folder) #for subdir in ['tss_tss']: # os.mkdir(os.path.join(dist_folder, subdir)) #chrom_list = list(range(1,6)) chrom_list = list(range(1,23)) +['X'] chrom_list.reverse() #chrom_list = list(range(3,23)) + ['X'] #chrom_list = list(range(1,3)) #chrom_list = ['X'] #chrom_list = ['X'] chrom_list = ['chr%s' %x for x in chrom_list ] gene_file_list = ["gencodev32_hg38_known_gene_%s_order.tsv" %x for x in chrom_list] for gene_file in gene_file_list: print (gene_file) gene_data = pd.read_csv("%s/%s" %(gene_folder, gene_file), sep='\t', header=0, index_col=None) gene_data.set_index('Gene stable ID', inplace=True) #exp_file = parse_expression_hd5_lee("/data/johlee/CoCoCoNet/networks/human_prioAggNet.hdf5") #subset_cols = list (set(gene_data.index) & set(exp_file.index)) #gene_data = gene_data.loc[subset_cols, :] if not gene_data.empty: chrom = gene_data.iloc[0]['chrom'] if chrom in chrom_list: SRA_list.append(ind_network) chr_list.append(chrom) resolution_list.append(resoulution_in_kb) prefix = '/'.join(dist_folder.split('/')[:-1]) RAWobserved = mmread("/grid/gillis/data/lohia/hi_c_data_processing/data/%s_%sKB/%s.mtx" %(ind_network, resoulution_in_kb, chrom)) #print ("/sonas-hs/gillis/hpc/data/lohia/hi_c_data_processing/data/%s/%s_%sKB/%s_binned.mtx" %(ind_network, ind_network, resoulution_in_kb, chrom)) maximum_read, map_resolution, maximum_bias = make_tss_tss(gene_data, chrom, dist_folder, RAWobserved, resoulution_in_kb, ind_network) maximum_read_list.append(maximum_read) map_resolution_list.append(map_resolution) maximum_bias_list.append(maximum_bias) SRA_info = pd.DataFrame(list(zip(SRA_list , chr_list , resolution_list , maximum_read_list , map_resolution_list , maximum_bias_list)), columns = ['SRA', 'chr', 'resolution', 'maximum_read', 'map_resolution', 'maximum_bias']) SRA_info.to_csv(os.path.join(dist_folder, f'SRR_info_{ind_network}_{resoulution_in_kb}.csv'), index=False, sep='\t') else: pass else: pass if __name__ == '__main__': print ("temp") #for resolution in [500, 50, 100, 10, 5]: resolution = sys.argv[1] ind_network = sys.argv[2] outdir = "/grid/gillis/data/lohia/hi_c_data_processing/software/SRR_info" #for resolution in [500, 50, 10]: make_distance_matrix('/grid/gillis/data/lohia/hi_c_data_processing/software/gene_tx_map', outdir, resoulution_in_kb=int(resolution), ind_network=ind_network)