docs: def njsd_geneset(network, ref, query, gene_set, file, verbose=True):
      """Compute gene set-specified nJSD between reference and query expression profiles.
      Attribute;
          network (str): File path to a network file.
          ref (str): File path to a reference expression file.
          query (str): File path to a query expression file.
          geneset (str): File path to a gene set file.
      """
      graph, gene_set_total = util.parse_network(network)
      ref_gene_expression_dict = util.parse_gene_expression(ref, mean=True)
      query_gene_expression_dict = util.parse_gene_expression(query, mean=False)
      group_gene_set_dict = util.parse_gene_set(gene_set)
      maximally_ambiguous_gene_experession_dict = util.get_maximally_ambiguous_network(query_gene_expression_dict)
      gene_set_present = set(query_gene_expression_dict.keys())
  
      with open(file, 'w') as outFile:
          print('Gene_set_ID', 'nJSD_NT', 'nJSD_TA', 'tITH', sep='\t', file=outFile)
  
      for group, gene_set in group_gene_set_dict.items():
          gene_set_to_be_analyzed = gene_set.intersection(gene_set_present)
          # If no genes are available for the group, just ignore it.
          if len(gene_set_to_be_analyzed) == 0:
              logger.warning('%s has no genes available for analysis. Ignoring the group.' % group)
              continue
  
          # If every gene has a single neighbor, just ignore it.
          if all([graph.degree(gene) == 1 for gene in gene_set_to_be_analyzed]):
              logger.warning('%s has no genes with enough neighbors. Ignoring the group.' % group)
              continue
  
          normal_to_tumor_njsd = entropy.njsd(network=graph,
                                              ref_gene_expression_dict=ref_gene_expression_dict,
                                              query_gene_expression_dict=query_gene_expression_dict,
                                              gene_set=gene_set)
  
          tumor_to_ambiguous_njsd = entropy.njsd(network=graph,
                                                 ref_gene_expression_dict=maximally_ambiguous_gene_experession_dict,
                                                 query_gene_expression_dict=query_gene_expression_dict,
                                                 gene_set=gene_set)
  
          tITH = normal_to_tumor_njsd / (normal_to_tumor_njsd + tumor_to_ambiguous_njsd)
          with open(file, 'a') as outFile:
              print(group, normal_to_tumor_njsd, tumor_to_ambiguous_njsd, tITH, sep='\t', file=outFile)