Purpose of the notebook is to generate reports via saved images

import os
import sys
sys.path.append("/oak/stanford/groups/akundaje/manyu/C2H2_ZNF_project/utils/modisco_utils")
sys.path.append(os.path.abspath("/oak/stanford/groups/akundaje/manyu/softwares/tfmodisco_tf_models/src/"))
sys.path.append(os.path.abspath("/oak/stanford/groups/akundaje/manyu/softwares/tfmodisco_tf_models/notebooks/reports/"))
from tfmodisco.run_tfmodisco import import_shap_scores, import_tfmodisco_results
from motif.read_motifs import pfm_info_content, trim_motif_by_ic
import motif.moods as moods
import plot.viz_sequence as viz_sequence
from util import figure_to_vdom_image, import_peak_table
import h5py
import pandas as pd
import numpy as np
import modisco
import sklearn.decomposition
import umap
import matplotlib.pyplot as plt
import matplotlib.font_manager as font_manager
import vdom.helpers as vdomh
from IPython.display import display
import tqdm
import io
import base64
import urllib
import deepdish as dd
import shutil
from vdom.helpers import h1, p, img, div, b

Main class to generate reports in HTML form using the saved images and h5 files

from shutil import copyfile
def get_read_handle(h5_file):
    '''Copy to temp with a hashed name and give the read handle'''
    name = h5_file.split('/')[-1]
    #copy_file = '/users/manyu/temp/{}_{}'.format(str(hash(h5_file)),name)
    copy_file = '/srv/scratch/manyu/tmp/{}_{}'.format(str(hash(h5_file)),name)
    
    if os.path.exists(copy_file):
        os.remove(copy_file)
    copyfile(h5_file,copy_file)
    f = h5py.File(copy_file,'r')
    return(f)


def close_h5file_read_handle(h5_file,read_handle):
    '''Removes the copied h5 file in temp and the read handle to avoid bloating'''
    read_handle.close()
    name = h5_file.split('/')[-1]
    #copy_file = '/users/manyu/temp/{}_{}'.format(str(hash(h5_file)),name)
    copy_file = '/srv/scratch/manyu/tmp/{}_{}'.format(str(hash(h5_file)),name)
    if os.path.exists(copy_file):
        os.remove(copy_file)
    
    
    
  

    

class ReportGenerator():
    def __init__(self,tf,dataset):
        self.dataset = dataset
        self.tf = tf
        # Path to saved modisco results on /oak. Images and h5 files
        self.reports_base = '/oak/stanford/groups/akundaje/manyu/C2H2_ZNF_project/train_profile_models_2020/modisco_reports'
        self.results_dir = '{}/{}/{}'.format(self.reports_base,dataset,tf)
        
        ## Paths to the saved modisco report results to copy/copied on Mitra: Images and h5 files
        self.modisco_reports_results_mitra = '/srv/www/kundaje/manyu/C2H2_ZNF_project/modisco_report_saved_information'
        self.mitra_results_dataset_dir = '{}/{}'.format(self.modisco_reports_results_mitra,self.dataset)
        if not os.path.exists(self.mitra_results_dataset_dir):
            os.makedirs(self.mitra_results_dataset_dir)
        self.save_results_mitra_dir  = '{}/{}'.format(self.mitra_results_dataset_dir,self.tf)
        
        self.main_patterns_dir = '{}/main_patterns'.format(self.save_results_mitra_dir)
        self.main_patterns = '{}/all_motifs.h5'.format(self.main_patterns_dir)
        
        self.subpatterns_dir  = '{}/sub_patterns'.format(self.save_results_mitra_dir)
        self.subpatterns = '{}/all_motif_subclusters.h5'.format(self.subpatterns_dir)
        
        
        ## Paths to the modisco reports html on lab cluster
        self.modisco_reports_base = '/srv/www/kundaje/manyu/C2H2_ZNF_project/modisco_reports'
        self.modisco_reports_dataset_base = '{}/{}'.format(self.modisco_reports_base,self.dataset)
        if not os.path.exists(self.modisco_reports_dataset_base):
            os.makedirs(self.modisco_reports_dataset_base)
        self.modisco_report_filename = '{}/{}.html'.format(self.modisco_reports_dataset_base,self.tf)
        
        ## Path to webpages and weblinks to images(needs weblinks to load; not local links)
        self.mitra_web_base = 'http://mitra.stanford.edu/kundaje/manyu/C2H2_ZNF_project/'
        self.mitra_web_results_base = '{}/modisco_report_saved_information/{}/{}'.format(self.mitra_web_base,self.dataset,self.tf)
        self.mitra_web_results_main_patterns = '{}/main_patterns'.format(self.mitra_web_results_base)
        self.mitra_web_results_sub_patterns = '{}/sub_patterns'.format(self.mitra_web_results_base)
        self.mitra_web_results_profiles = '{}/profiles_dir'.format(self.mitra_web_results_base)
        self.mitra_web_results_b1h_aln = '{}/b1h_aligments'.format(self.mitra_web_results_base)
        
        ##Do we need this here?
        #self.mitra_web_reports_base = '{}/modisco_reports'.format(self.mitra_web_base)
        
        
    
    def copy_to_mitra(self):
        assert(os.path.exists(self.results_dir))
        if os.path.exists(self.save_results_mitra_dir):
            shutil.rmtree(self.save_results_mitra_dir)
        shutil.copytree(self.results_dir,self.save_results_mitra_dir)
        print('Finished copying files to Mitra')
        permissions_cmd = 'chmod -R 777 {}'.format(self.save_results_mitra_dir)
        os.system(permissions_cmd)
        
        
        
    def show_subcluster_table(self):
      
        assert(os.path.exists(self.main_patterns))
        
        ## Opening a read handle for h5. Remmeber to close it also
        try:
            patterns_read_handle = get_read_handle(self.main_patterns)
        except Exception as e:
            print(e)
        
        try:
            sub_patterns_read_handle = get_read_handle(self.subpatterns)
        except Exception as e:
            print(e)

        num_patterns  = len(patterns_read_handle.keys())

     
        
        for pattern_i in range(num_patterns):
            
            display(vdomh.h3("Pattern {}/{}".format(str(pattern_i + 1), num_patterns)))
        
            

            # Sort this out:
            colgroup = vdomh.colgroup(
                vdomh.col(style={"width": "5%"}),
                vdomh.col(style={"width": "5%"}),
                vdomh.col(style={"width": "50%"}),
                vdomh.col(style={"width": "40%"})
            )
            header = vdomh.thead(
                vdomh.tr(
                    vdomh.th("Subpattern", style={"text-align": "center"}),
                    vdomh.th("Seqlets", style={"text-align": "center"}),
                    vdomh.th("Embeddings", style={"text-align": "center"}),
                    vdomh.th("CWM", style={"text-align": "center"}),
                    vdomh.th("PFM", style={"text-align": "center"})
                )
            )
            
            
            ## Fill in the paths and seqlets of the main pattern
            total_seqlets = patterns_read_handle[str(pattern_i)]['n_seqlets'].value
            emb_agg_path = '{}/{}_subcluster_agg.png'.format(self.mitra_web_results_sub_patterns,str(pattern_i))
            pwm_path = '{}/{}_pfm_full.png'.format(self.mitra_web_results_main_patterns,str(pattern_i))
            cwm_path = '{}/{}_cwm_full.png'.format(self.mitra_web_results_main_patterns,str(pattern_i))
            
            
            table_rows = [vdomh.tr(
                vdomh.td("Agg."),
                vdomh.td(str(total_seqlets)),
                vdomh.td(img(src=emb_agg_path)),
                vdomh.td(img(src=cwm_path)),
                vdomh.td(img(src=pwm_path))
            )]

            
            total_sub_pats = len(sub_patterns_read_handle[str(pattern_i)].keys())
            for subpat_i in range(total_sub_pats):
                emb_path = '{}/{}_subcluster_{}.png'.format(self.mitra_web_results_sub_patterns,str(pattern_i),str(subpat_i))
            
                total_seqlets = sub_patterns_read_handle[str(pattern_i)]['subcluster_{}'.format(str(subpat_i))]['n_seqlets'].value
                pwm_path = '{}/{}_subcluster_{}_pfm_trimmed.png'.format(self.mitra_web_results_sub_patterns,str(pattern_i),str(subpat_i))
                cwm_path = '{}/{}_subcluster_{}_cwm_trimmed.png'.format(self.mitra_web_results_sub_patterns,str(pattern_i),str(subpat_i))
                table_rows.append(vdomh.tr(
                vdomh.td('{}'.format(str(subpat_i))),
                vdomh.td(str(total_seqlets)),
                vdomh.td(img(src=emb_path)),
                vdomh.td(img(src=cwm_path)),
                vdomh.td(img(src=pwm_path))
            ))
            table = vdomh.table(header,vdomh.tbody(*table_rows))
            
            display(table)
            
        #Remove unnecessary temp files
        close_h5file_read_handle(self.main_patterns,patterns_read_handle)
        close_h5file_read_handle(self.subpatterns,sub_patterns_read_handle)
    
    
    def show_footprints_and_peak_dist(self):
        try:
            patterns_read_handle = get_read_handle(self.main_patterns)
        except Exception as e:
            print(e)
        
        num_patterns  = len(patterns_read_handle.keys())
        for pattern_i in range(num_patterns):
            
            display(vdomh.h3("Pattern {}/{} Footprint".format(str(pattern_i + 1), num_patterns)))
            footprint_plot = '{}/profiles_pattern_{}.png'.format(self.mitra_web_results_profiles,str(pattern_i))
            motif_distplot  = '{}/motif_distances_pattern_{}.png'.format(self.mitra_web_results_profiles,str(pattern_i))
            display(img(src=footprint_plot))
            display(vdomh.h3("Pattern {}/{} Motif Distance Distribution".format(str(pattern_i + 1), num_patterns)))
            display(img(src=motif_distplot))
    
    
    def show_b1h_alignments(self):
        try:
            patterns_read_handle = get_read_handle(self.main_patterns)
        except Exception as e:
            print(e)
        


        num_patterns  = len(patterns_read_handle.keys())
        header = vdomh.thead(
            vdomh.tr(
                vdomh.th("Pattern", style={"text-align": "center"}),
                vdomh.th("N-Seqlets", style={"text-align": "center"}),
                vdomh.th("Alignment", style={"text-align": "center"})

            )
        )

         
        display(vdomh.h2("B1H alignments of all patterns"))
        table_rows = []
        for pattern_i in range(num_patterns):
            
            
        
            total_seqlets = patterns_read_handle[str(pattern_i)]['n_seqlets'].value
            alignment_img = '{}/pattern_{}_aligned.png'.format(self.mitra_web_results_b1h_aln,str(pattern_i))
            

            table_rows.append(vdomh.tr(
            vdomh.td('{}'.format(str(pattern_i+1))),
            vdomh.td('{}'.format(str(total_seqlets))),
            vdomh.td(img(src=alignment_img)),

            ))
        close_h5file_read_handle(self.main_patterns,patterns_read_handle)
        table = vdomh.table(header,vdomh.tbody(*table_rows))
        display(table)

# Define parameters/fetch arguments
dataset = os.environ["DATASET"]
tf = os.environ["TF"]

# tf = 'ENCSR904JEY'
# dataset = 'ENCODE_new'

vdomh.h1('Running {} in dataset: {}'.format(tf,dataset))

Running ENCSR701AQS in dataset: ENCODE_new

a = ReportGenerator(tf=tf,dataset=dataset)

Show motifs and submotifs and relevant embeddings

a.show_subcluster_table()

Pattern 1/30

/users/manyu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:127: H5pyDeprecationWarning: dataset.value has been deprecated. Use dataset[()] instead.
/users/manyu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:146: H5pyDeprecationWarning: dataset.value has been deprecated. Use dataset[()] instead.

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	1162
0	252
1	205
2	197
3	185
4	164
5	159

Pattern 2/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	1074
0	124
1	121
2	114
3	113
4	106
5	93
6	78
7	64
8	62
9	56
10	50
11	48
12	45

Pattern 3/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	571
0	78
1	75
2	65
3	51
4	47
5	44
6	44
7	43
8	43
9	40
10	35
11	4
12	2

Pattern 4/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	528
0	78
1	69
2	66
3	60
4	52
5	52
6	49
7	40
8	40
9	22

Pattern 5/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	373
0	54
1	53
2	50
3	47
4	39
5	37
6	30
7	24
8	18
9	17
10	4

Pattern 6/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	322
0	80
1	79
2	66
3	52
4	31
5	11
6	3

Pattern 7/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	265
0	56
1	50
2	49
3	40
4	40
5	28
6	2

Pattern 8/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	260
0	49
1	36
2	31
3	28
4	27
5	25
6	22
7	22
8	20

Pattern 9/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	201
0	44
1	44
2	34
3	27
4	25
5	19
6	5
7	3

Pattern 10/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	197
0	68
1	67
2	28
3	26
4	8

Pattern 11/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	201
0	42
1	38
2	32
3	28
4	25
5	22
6	7
7	5
8	2

Pattern 12/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	169
0	50
1	38
2	28
3	23
4	19
5	11

Pattern 13/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	115
0	30
1	25
2	20
3	17
4	13
5	7
6	3

Pattern 14/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	105
0	26
1	24
2	15
3	15
4	13
5	12

Pattern 15/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	101
0	25
1	22
2	19
3	14
4	13
5	8

Pattern 16/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	103
0	32
1	25
2	21
3	13
4	12

Pattern 17/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	99
0	30
1	28
2	25
3	13
4	3

Pattern 18/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	93
0	20
1	20
2	18
3	16
4	12
5	4
6	3

Pattern 19/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	97
0	24
1	18
2	18
3	17
4	12
5	8

Pattern 20/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	79
0	29
1	28
2	22

Pattern 21/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	73
0	23
1	18
2	17
3	10
4	5

Pattern 22/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	80
0	28
1	26
2	20
3	4
4	2

Pattern 23/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	117
0	33
1	25
2	24
3	18
4	17

Pattern 24/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	52
0	21
1	15
2	9
3	7

Pattern 25/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	64
0	16
1	15
2	15
3	9
4	9

Pattern 26/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	60
0	23
1	16
2	16
3	3
4	2

Pattern 27/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	51
0	23
1	17
2	9
3	2

Pattern 28/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	55
0	18
1	16
2	11
3	10

Pattern 29/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	51
0	17
1	10
2	9
3	8
4	7

Pattern 30/30

Subpattern	Seqlets	Embeddings	CWM	PFM
Agg.	45
0	18
1	17
2	10

Show footprints and motif distributions within peaks

a.show_footprints_and_peak_dist()

Pattern 1/30 Footprint

Pattern 1/30 Motif Distance Distribution

Pattern 2/30 Footprint

Pattern 2/30 Motif Distance Distribution

Pattern 3/30 Footprint

Pattern 3/30 Motif Distance Distribution

Pattern 4/30 Footprint

Pattern 4/30 Motif Distance Distribution

Pattern 5/30 Footprint

Pattern 5/30 Motif Distance Distribution

Pattern 6/30 Footprint

Pattern 6/30 Motif Distance Distribution

Pattern 7/30 Footprint

Pattern 7/30 Motif Distance Distribution

Pattern 8/30 Footprint

Pattern 8/30 Motif Distance Distribution

Pattern 9/30 Footprint

Pattern 9/30 Motif Distance Distribution

Pattern 10/30 Footprint

Pattern 10/30 Motif Distance Distribution

Pattern 11/30 Footprint

Pattern 11/30 Motif Distance Distribution

Pattern 12/30 Footprint

Pattern 12/30 Motif Distance Distribution

Pattern 13/30 Footprint

Pattern 13/30 Motif Distance Distribution

Pattern 14/30 Footprint

Pattern 14/30 Motif Distance Distribution

Pattern 15/30 Footprint

Pattern 15/30 Motif Distance Distribution

Pattern 16/30 Footprint

Pattern 16/30 Motif Distance Distribution

Pattern 17/30 Footprint

Pattern 17/30 Motif Distance Distribution

Pattern 18/30 Footprint

Pattern 18/30 Motif Distance Distribution

Pattern 19/30 Footprint

Pattern 19/30 Motif Distance Distribution

Pattern 20/30 Footprint

Pattern 20/30 Motif Distance Distribution

Pattern 21/30 Footprint

Pattern 21/30 Motif Distance Distribution

Pattern 22/30 Footprint

Pattern 22/30 Motif Distance Distribution

Pattern 23/30 Footprint

Pattern 23/30 Motif Distance Distribution

Pattern 24/30 Footprint

Pattern 24/30 Motif Distance Distribution

Pattern 25/30 Footprint

Pattern 25/30 Motif Distance Distribution

Pattern 26/30 Footprint

Pattern 26/30 Motif Distance Distribution

Pattern 27/30 Footprint

Pattern 27/30 Motif Distance Distribution

Pattern 28/30 Footprint

Pattern 28/30 Motif Distance Distribution

Pattern 29/30 Footprint

Pattern 29/30 Motif Distance Distribution

Pattern 30/30 Footprint

Pattern 30/30 Motif Distance Distribution

Show B1H alignments of main motifs

a.show_b1h_alignments()

B1H alignments of all patterns

/users/manyu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:207: H5pyDeprecationWarning: dataset.value has been deprecated. Use dataset[()] instead.

Pattern	N-Seqlets	Alignment
1	1162
2	1074
3	571
4	528
5	373
6	322
7	265
8	260
9	201
10	197
11	201
12	169
13	115
14	105
15	101
16	103
17	99
18	93
19	97
20	79
21	73
22	80
23	117
24	52
25	64
26	60
27	51
28	55
29	51
30	45