<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def njsd(network, ref_gene_expression_dict, query_gene_expression_dict, gene_set):

 """Calculate Jensen-Shannon divergence between query and reference gene expression profile.
 """

    gene_jsd_dict = dict()
    reference_genes = ref_gene_expression_dict.keys()
    assert len(reference_genes) != 'Reference gene expression profile should have > 0 genes.'

    for gene in gene_set:
        if gene not in network.nodes:
            continue

        neighbors = find_neighbors(network, gene)

        query_expression_vec = get_neighbor_expression_vector(neighbors, query_gene_expression_dict)
        ref_expression_vec = get_neighbor_expression_vector(neighbors, ref_gene_expression_dict)
        assert len(query_expression_vec) == len(ref_expression_vec), 'Topology of reference network and query network differs. Please check.'

        # A gene which has non-expressed neighbors is ignored.
        if np.sum(query_expression_vec) == 0 and np.sum(ref_expression_vec) == 0:
            continue

        query_p_vec = exp2prob(query_expression_vec)
        ref_p_vec = exp2prob(ref_expression_vec)

        gene_jsd_dict[gene] = jsd(query_p_vec, ref_p_vec)

    return np.mean(list(gene_jsd_dict.values()))