/
classify_bedregions.py
147 lines (116 loc) · 5.23 KB
/
classify_bedregions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
import argparse
import os
import numpy as np
from keras import backend as K
from keras.layers import Conv2D
from pkg_resources import resource_filename
from janggu import Janggu
from janggu import Scorer
from janggu import inputlayer
from janggu import outputconv
from janggu.data import Bioseq
from janggu.data import Cover
from janggu.layers import DnaConv2D
from janggu.layers import LocalAveragePooling2D
from janggu.utils import ExportClustermap
from janggu.utils import ExportTsne
from janggu.utils import ExportTsv
np.random.seed(1234)
# Fetch parser arguments
PARSER = argparse.ArgumentParser(description='Command description.')
PARSER.add_argument('-path', dest='path',
default='tf_results',
help="Output directory for the examples.")
PARSER.add_argument('-order', dest='order', type=int,
default=1,
help="One-hot order.")
PARSER.add_argument('-epochs', dest='epochs', type=int,
default=100,
help="Number of epochs to train.")
args = PARSER.parse_args()
os.environ['JANGGU_OUTPUT'] = args.path
# load the dataset
# The pseudo genome represents just a concatenation of all sequences
# in sample.fa and sample2.fa. Therefore, the results should be almost
# identically to the models obtained from classify_fasta.py.
REFGENOME = resource_filename('janggu', 'resources/pseudo_genome.fa')
# ROI contains regions spanning positive and negative examples
ROI_TRAIN_FILE = resource_filename('janggu', 'resources/roi_train.bed')
ROI_TEST_FILE = resource_filename('janggu', 'resources/roi_test.bed')
# PEAK_FILE only contains positive examples
PEAK_FILE = resource_filename('janggu', 'resources/scores.bed')
# Training input and labels are purely defined genomic coordinates
DNA = Bioseq.create_from_refgenome('dna', refgenome=REFGENOME,
roi=ROI_TRAIN_FILE,
binsize=200,
order=args.order,
cache=True)
LABELS = Cover.create_from_bed('peaks', roi=ROI_TRAIN_FILE,
bedfiles=PEAK_FILE,
binsize=200,
resolution=200,
cache=True,
storage='sparse')
DNA_TEST = Bioseq.create_from_refgenome('dna', refgenome=REFGENOME,
roi=ROI_TEST_FILE,
binsize=200,
order=args.order)
LABELS_TEST = Cover.create_from_bed('peaks',
bedfiles=PEAK_FILE,
roi=ROI_TEST_FILE,
binsize=200,
resolution=200,
storage='sparse')
@inputlayer
@outputconv('sigmoid')
def double_stranded_model_dnaconv(inputs, inp, oup, params):
""" keras model for scanning both DNA strands.
A more elegant way of scanning both strands for motif occurrences
is achieved by the DnaConv2D layer wrapper, which internally
performs the convolution operation with the normal kernel weights
and the reverse complemented weights.
"""
with inputs.use('dna') as layer:
# the name in inputs.use() should be the same as the dataset name.
layer = DnaConv2D(Conv2D(params[0], (params[1], 1),
activation=params[2]))(layer)
output = LocalAveragePooling2D(window_size=layer.shape.as_list()[1],
name='motif')(layer)
return inputs, output
modeltemplate = double_stranded_model_dnaconv
K.clear_session()
# create a new model object
model = Janggu.create(template=modeltemplate,
modelparams=(30, 21, 'relu'),
inputs=DNA,
outputs=LABELS)
model.compile(optimizer='adadelta', loss='binary_crossentropy',
metrics=['acc'])
hist = model.fit(DNA, LABELS, epochs=100)
print('#' * 40)
print('loss: {}, acc: {}'.format(hist.history['loss'][-1],
hist.history['acc'][-1]))
print('#' * 40)
# clustering plots based on hidden features
heatmap_eval = Scorer('heatmap', exporter=ExportClustermap(z_score=1.))
tsne_eval = Scorer('tsne', exporter=ExportTsne())
# output the predictions as tables or json files
pred_tsv = Scorer('pred', exporter=ExportTsv(row_names=DNA_TEST.gindexer.chrs))
# do the evaluation on the independent test data
model.evaluate(DNA_TEST, LABELS_TEST, datatags=['test'],
callbacks=['auc', 'prc', 'roc', 'auprc'])
pred = model.predict(DNA_TEST)
cov_pred = Cover.create_from_array('BindingProba', pred, LABELS_TEST.gindexer)
print('Oct4 predictions scores should be greater than Mafk scores:')
print('Prediction score examples for Oct4')
for i in range(4):
print('{}.: {}'.format(i, pred[i]))
print('Prediction score examples for Mafk')
for i in range(1, 5):
print('{}.: {}'.format(i, pred[-i]))
model.predict(DNA_TEST, datatags=['test'],
callbacks=[pred_tsv],
layername='motif')
model.predict(DNA_TEST, datatags=['test'],
callbacks=[heatmap_eval, tsne_eval],
layername='motif')