refine_docs

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)
 
 

missingfiles.txt

@@ -0,0 +1,3 @@
+data/enhancers.gas.class.tsv
+data/gene_tss.uniq.tsv
+raw_data/hic/5kb_resolution_intrachromosomal

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

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,