diff --git a/README.md b/README.md index c4b128a..f835302 100644 --- a/README.md +++ b/README.md @@ -1,25 +1,20 @@ # sleps Supervised Learning of Enhancer Promoter Specificity. +## Overview +Sleps software is a pipeline that utilizes machine learning and graph relationships to predict the effect of an enhancer on a particular promoter. -## Data +Currently, the software is in alpha and needs further development to generalize enough for training on new datasets. However, our experiment and trained sleps model can be reproduced using this guide. There are three main steps needed to complete the training of the model. -Data download: [link](https://drive.google.com/drive/folders/1270MmEk8oF3VpJ5llkaKSCqPy0i8-qLW?usp=drive_link) +1. [Download Data](https://drive.google.com/drive/folders/1270MmEk8oF3VpJ5llkaKSCqPy0i8-qLW?usp=drive_link) + - First, download requisite data. Training data includes Hi-C, DHS, H3K27ac, and CRISPRi datasets as well as a suite of ChIPSeq data from ENCODE. +2. [Generate Enhancer Networks](https://github.com/HanLabUNLV/abic/blob/master/network_generation_process.md) + - Next generate networks using Hi-C data and tools. This step requires a singularity environment that we have provided as well as HPC capabilities. +3. [Train Model & Apply](https://github.com/HanLabUNLV/abic/blob/master/learning.md) + - Lastly train the model again using singularity and an HPC. With the trained model in hand, you can apply the model as well to new data. -## Network Generation - -container used during network generation can be downloaded here: -[network singularity container](https://drive.google.com/drive/folders/13WP9gLttNaa3HQLAs5Of-PB4ZVqbmwUJ?usp=sharing) - -documentation: -[How to generate enhancer networks](https://github.com/HanLabUNLV/abic/blob/master/network_generation_process.md) - - -## Machine Learning -container used during learning can be downloaded here: -[learning singularity container](https://drive.google.com/drive/folders/1QTNEvYx6T5kXfspyx4w_OEKo8dJ8cJTQ?usp=sharing) - -documentation: -[How to train and apply the xgboost model to predict positive EP pairs](https://github.com/HanLabUNLV/abic/blob/master/learning.md) +## Singularity Containers +- [Network singularity container](https://drive.google.com/drive/folders/13WP9gLttNaa3HQLAs5Of-PB4ZVqbmwUJ?usp=sharing) +- [Machine Learning singularity container](https://drive.google.com/drive/folders/1QTNEvYx6T5kXfspyx4w_OEKo8dJ8cJTQ?usp=sharing) diff --git a/missingfiles.txt b/missingfiles.txt new file mode 100644 index 0000000..fb70a55 --- /dev/null +++ b/missingfiles.txt @@ -0,0 +1,3 @@ +data/enhancers.gas.class.tsv +data/gene_tss.uniq.tsv +raw_data/hic/5kb_resolution_intrachromosomal diff --git a/src/network/hic.py b/src/network/hic.py index 89cac41..5dac88a 100644 --- a/src/network/hic.py +++ b/src/network/hic.py @@ -136,10 +136,19 @@ def apply_kr_threshold(hic_mat, hic_norm_file, kr_cutoff): #Convert all entries in the hic matrix corresponding to low kr norm entries to NaN #Note that in scipy sparse matrix multiplication 0*nan = 0 #So this doesn't convert 0's to nan only nonzero to nan - + norms = np.loadtxt(hic_norm_file) - norms[norms < kr_cutoff] = np.nan - norms[norms >= kr_cutoff] = 1 + + #NANs are fine, but break the code, fixed it + for i in range(0,len(norms)): + if norms[i]!=norms[i]: + pass + elif norms[i] < kr_cutoff: + norms[i]=np.nan + else: + norms[i] = 1 + #norms[norms < kr_cutoff] = np.nan + #norms[norms >= kr_cutoff] = 1 norm_mat = ssp.dia_matrix(( 1.0/norms, [0]), (len(norms), len(norms))) return norm_mat * hic_mat * norm_mat diff --git a/src/network/predictor.py b/src/network/predictor.py index f935b06..45b6261 100644 --- a/src/network/predictor.py +++ b/src/network/predictor.py @@ -471,7 +471,7 @@ def validated_network_around_gene(enhancers, gene, valid_loops): print([v for v in network.vs]) exit('network has less than 2 vertices') add_attr_network(network, enhancers, gene) - hic_dir = 'raw_data/hic/5kb_resolution_intrachromosomal/' + hic_dir = 'data/raw_data/hic/5kb_resolution_intrachromosomal/' hic_resolution = 5000 hic_type = 'juicebox' tss_hic_contribution = 100 @@ -480,8 +480,8 @@ def validated_network_around_gene(enhancers, gene, valid_loops): args = {'hic_dir':hic_dir,'hic_resolution':hic_resolution,'hic_type':hic_type,'tss_hic_contribution':tss_hic_contribution,'window':window,'hic_gamma':hic_gamma} chromosome = enhancers.chr.tolist()[0] #hic_file, hic_norm_file, hic_is_vc = get_hic_file(chromosome, args['hic_dir'], hic_type = args['hic_type']) - hic_file = 'raw_data/hic/5kb_resolution_intrachromosomal/'+chromosome+'/'+chromosome+'_5kb.RAWobserved' - hic_norm_file = 'raw_data/hic/5kb_resolution_intrachromosomal/'+chromosome+'/'+chromosome+'_5kb.KRnorm' + hic_file = 'data/raw_data/hic/5kb_resolution_intrachromosomal/'+chromosome+'/'+chromosome+'_5kb.RAWobserved' + hic_norm_file = 'data/raw_data/hic/5kb_resolution_intrachromosomal/'+chromosome+'/'+chromosome+'_5kb.KRnorm' hic_is_vc = True hic_threshold = 0.#2.788652e-03 * 0.01 #75% of contacts #load HIC matrix @@ -572,7 +572,7 @@ def network_around_gene(enhancers, gene): print([v for v in network.vs]) exit('network has less than 2 vertices') add_attr_network(network, enhancers, gene) - hic_dir = 'raw_data/hic/5kb_resolution_intrachromosomal/' + hic_dir = 'data/raw_data/hic/5kb_resolution_intrachromosomal/' hic_resolution = 5000 hic_type = 'juicebox' tss_hic_contribution = 100 @@ -581,10 +581,10 @@ def network_around_gene(enhancers, gene): args = {'hic_dir':hic_dir,'hic_resolution':hic_resolution,'hic_type':hic_type,'tss_hic_contribution':tss_hic_contribution,'window':window,'hic_gamma':hic_gamma} chromosome = enhancers.chr.tolist()[0] #hic_file, hic_norm_file, hic_is_vc = get_hic_file(chromosome, args['hic_dir'], hic_type = args['hic_type']) - hic_file = 'raw_data/hic/5kb_resolution_intrachromosomal/'+chromosome+'/'+chromosome+'_5kb.RAWobserved' - hic_norm_file = 'raw_data/hic/5kb_resolution_intrachromosomal/'+chromosome+'/'+chromosome+'_5kb.KRnorm' + hic_file = 'data/raw_data/hic/5kb_resolution_intrachromosomal/'+chromosome+'/'+chromosome+'_5kb.RAWobserved' + hic_norm_file = 'data/raw_data/hic/5kb_resolution_intrachromosomal/'+chromosome+'/'+chromosome+'_5kb.KRnorm' hic_is_vc = True - hic_threshold = 0.#2.788652e-03 * 0.01 #75% of contacts + hic_threshold = 0. #2.788652e-03 * 0.01 #75% of contacts #load HIC matrix HiC = load_hic(hic_file = hic_file, hic_norm_file = hic_norm_file,