<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Initializing filters on known motifs
In the scenario where data is scarse, it is often useful to initialize the filters of the first convolutional layer to some known position weights matrices (PWM's). That way, the model already starts with a parameter configuration much closer to the 'right' one.
Concise provides access to 2 PWM databases
Step1:  Let's choose PUM2 PWM (RBP in Human)
Step2:  Visualization - PWM class
The PWM class provides a method plotPWM to visualize the PWM.
Step3:  We can select the PWM with id 129.
Step4:  Initialize the conv filters with PWM values
Step5:  ci.PSSMKernelInitializer will set the filters of the first convolutional layer to the values of the position-specific scoring matrix (PSSM)
Step6:  Test-set performance
Step7:  Filter visualization
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
# RBP PWM's
from concise.data import attract
dfa = attract.get_metadata()
dfa
# TF PWM's
from concise.data import encode
dfe = encode.get_metadata()
dfe
# TF PWM's
from concise.data import hocomoco
dfh = hocomoco.get_metadata()
dfh
Explanation: Initializing filters on known motifs
In the scenario where data is scarse, it is often useful to initialize the filters of the first convolutional layer to some known position weights matrices (PWM's). That way, the model already starts with a parameter configuration much closer to the 'right' one.
Concise provides access to 2 PWM databases:
transcription factors from ENCODE (2067 PWMs)
transcription factors from HOCOMOCO v10 (640 PWMs)
rna-binding proteins from ATtrACT (1583 PWMs)
Find the motif of interest
Each PWM database is provided as a module under concise.data. It provides two functions:
concise.data.&lt;db&gt;.get_metadata() - returns a pandas.DataFrame with metadata information about each PWM 
concise.data.&lt;db&gt;.get_pwm_list() - given a list of PWM ids, return a list with concise.utils.pwm.PWM instances
Metadata tables
End of explanation
dfa_pum2 = dfa[dfa.Gene_name.str.match("PUM2") & \
               dfa.Organism.str.match("Homo_sapiens")]
dfa_pum2
Explanation: Let's choose PUM2 PWM (RBP in Human):
End of explanation
# Visualize the PUM2 Motifs from different experiments
from concise.utils.pwm import PWM
dfa_pum2_uniq = dfa_pum2[["Experiment_description", "PWM_id"]].drop_duplicates()
pwm_list = attract.get_pwm_list(dfa_pum2_uniq.PWM_id)
for i, pwm in enumerate(pwm_list):
    print("PWM_id:", pwm.name, "; Experiment_description:", dfa_pum2_uniq.Experiment_description.iloc[i])
    pwm.plotPWM(figsize=(3,1))
Explanation: Visualization - PWM class
The PWM class provides a method plotPWM to visualize the PWM.
End of explanation
pwm_list = [pwm for pwm in pwm_list if pwm.name == "129"]
pwm_list
Explanation: We can select the PWM with id 129.
End of explanation
import concise.layers as cl
import keras.layers as kl
import concise.initializers as ci
import concise.regularizers as cr
from keras.callbacks import EarlyStopping
from concise.preprocessing import encodeDNA
from keras.models import Model, load_model
from keras.optimizers import Adam
# get the data
def load(split="train", st=None):
    dt = pd.read_csv("../data/RBP/PUM2_{0}.csv".format(split))
    # DNA/RNA sequence
    xseq = encodeDNA(dt.seq) # list of sequences -> np.ndarray
    # response variable
    y = dt.binding_site.as_matrix().reshape((-1, 1)).astype("float")
    return {"seq": xseq}, y
train, valid, test = load("train"), load("valid"), load("test")
# deduce sequence length
seq_length = train[0]["seq"].shape[1]
# define the model
def model(train, filters=1, kernel_size=9, pwm_list=None, lr=0.001):
    seq_length = train[0]["seq"].shape[1]
    if pwm_list is None:
        kinit = "glorot_uniform"
        binit = "zeros"
    else:
        kinit = ci.PSSMKernelInitializer(pwm_list, add_noise_before_Pwm2Pssm=True)
        binit = "zeros"
        
    # sequence
    in_dna = cl.InputDNA(seq_length=seq_length, name="seq")
    x = cl.ConvDNA(filters=filters, 
                   kernel_size=kernel_size, 
                   activation="relu",
                   kernel_initializer=kinit,
                   bias_initializer=binit,
                   name="conv1")(in_dna)
    x = kl.AveragePooling1D(pool_size=4)(x)
    x = kl.Flatten()(x)
    
    x = kl.Dense(units=1)(x)
    m = Model(in_dna, x)
    m.compile(Adam(lr=lr), loss="binary_crossentropy", metrics=["acc"])
    return m
Explanation: Initialize the conv filters with PWM values
End of explanation
# create two models: with and without PWM initialization
m_rand_init = model(train, filters=3, pwm_list=None) # random initialization
m_pwm_init = model(train, filters=3, pwm_list=pwm_list) # motif initialization
print("Random initialization:")
m_rand_init.get_layer("conv1").plot_weights(figsize=(3, 5));
print("Known PWM initialization:")
m_pwm_init.get_layer("conv1").plot_weights(figsize=(3, 5));
# train the models
m_rand_init.fit(train[0], train[1], epochs=50, validation_data=valid, 
                verbose=0,
                callbacks=[EarlyStopping(patience=5)])
m_pwm_init.fit(train[0], train[1], epochs=50, validation_data=valid, 
                verbose=0,
                callbacks=[EarlyStopping(patience=5)]);
Explanation: ci.PSSMKernelInitializer will set the filters of the first convolutional layer to the values of the position-specific scoring matrix (PSSM):
$$ pssm_{ij} = log \frac{pwm_{ij}}{b_j} \;,$$
where $b_j$ is the background probability of observing base $j$.
We add gaussian noise to each individual filter. Let's visualize the filters:
End of explanation
import concise.eval_metrics as cem
# performance on the test-set
# Random initialization
print("Random intiailzation auPR:", cem.auprc(test[1], m_rand_init.predict(test[0])))
# PWM initialization
print("Known PWM initialization auPR:", cem.auprc(test[1], m_pwm_init.predict(test[0])))
Explanation: Test-set performance
End of explanation
m_rand_init.get_layer("conv1").plot_weights(plot_type="motif_pwm_info", figsize=(3, 5));
m_pwm_init.get_layer("conv1").plot_weights(plot_type="motif_pwm_info", figsize=(3, 5));
Explanation: Filter visualization
End of explanation
<END_TASK>