Repository: vanheeringen-lab/ANANSE

Issue:
ANANSE network: KeyError: 'activity'
command:
```
ananse network -g GRCz11.fa -a GRCz11.annotation.bed -n 12 -b binding_24hpf.tsv -e 24hpf_rep1.tsv 24hpf_rep2.tsv -o network_24hpf.txt
```

input:
```
$head -n 2 binding_24hpf.tsv
factor  enhancer        binding
SOX15   1:5842-6042     0.11866

$head -n 2 24hpf_rep1.tsv
gene_name       tpm
fam162a 19.823413000000002
```

full error message:
```
2021-06-29 18:14:58 | INFO | Loading expression
2021-06-29 18:18:00 | INFO | creating expression dataframe
2021-06-29 18:18:00 | INFO | Aggregate binding
2021-06-29 18:18:00 | INFO | reading enhancers
2021-06-29 18:20:47 | INFO | Reading binding file...
2021-06-29 18:20:50 | INFO | Grouping by tf and target gene...
2021-06-29 18:20:50 | INFO | Done grouping...
2021-06-29 18:20:50 | INFO | Reading factor activity
Traceback (most recent call last):
  File "/vol/mbconda/siebrenf/envs/ananse/lib/python3.9/site-packages/pandas/core/indexes/base.py", line 3080, in get_loc
    return self._engine.get_loc(casted_key)
  File "pandas/_libs/index.pyx", line 70, in pandas._libs.index.IndexEngine.get_loc
  File "pandas/_libs/index.pyx", line 101, in pandas._libs.index.IndexEngine.get_loc
  File "pandas/_libs/hashtable_class_helper.pxi", line 4554, in pandas._libs.hashtable.PyObjectHashTable.get_item
  File "pandas/_libs/hashtable_class_helper.pxi", line 4562, in pandas._libs.hashtable.PyObjectHashTable.get_item
KeyError: 'activity'

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/vol/mbconda/siebrenf/envs/ananse/bin/ananse", line 318, in <module>
    args.func(args)
  File "/vol/mbconda/siebrenf/envs/ananse/lib/python3.9/site-packages/ananse/commands/network.py", line 28, in network
    b.run_network(
  File "/vol/mbconda/siebrenf/envs/ananse/lib/python3.9/site-packages/ananse/network.py", line 558, in run_network
    act["activity"] = minmax_scale(rankdata(act["activity"], method="min"))
  File "/vol/mbconda/siebrenf/envs/ananse/lib/python3.9/site-packages/pandas/core/frame.py", line 3024, in __getitem__
    indexer = self.columns.get_loc(key)
  File "/vol/mbconda/siebrenf/envs/ananse/lib/python3.9/site-packages/pandas/core/indexes/base.py", line 3082, in get_loc
    raise KeyError(key) from err
KeyError: 'activity'
```

[responsible code ](https://github.com/vanheeringen-lab/ANANSE/blob/cd1f3fef7fcf486329a782748da8a099362fcf4f/ananse/network.py#L555)

I can't figure out where "factor_activity.tsv" is made. It must exist, but it seems to be missing a column.

Discussion:
codeclimate[bot]: Code Climate has analyzed commit e3455307 and detected **0 issues** on this pull request.

The test coverage on the diff in this pull request is **91.6%** (**50%** is the threshold).

This pull request will bring the total coverage in the repository to **43.7%**.

View more on [Code Climate](https://codeclimate.com/github/vanheeringen-lab/ANANSE/pull/102).


Patch:
diff --git a/.gitignore b/.gitignore
index b669bb4..8623dbc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,6 +9,7 @@ dist/
 site/
 __pycache__
 **/.ipynb_checkpoints/
+dask-worker-space/
 *egg-info
 bld.bat
 build.sh
diff --git a/ananse/network.py b/ananse/network.py
index b7a0d12..b6f9e83 100644
--- a/ananse/network.py
+++ b/ananse/network.py
@@ -24,13 +24,10 @@ from tempfile import NamedTemporaryFile
 from dask.distributed import progress
 from loguru import logger
 
-from genomepy import Genome
 import pyranges as pr
 
-import ananse
-from ananse.utils import region_gene_overlap
-
 warnings.filterwarnings("ignore")
+PACKAGE_DIR = os.path.dirname(__file__)
 
 
 class Network(object):
@@ -42,25 +39,24 @@ class Network(object):
         include_promoter=False,
         include_enhancer=True,
     ):
-        """[infer cell type-specific gene regulatory network]
-
-        Arguments:
-            object {[type]} -- [description]
-
-        Keyword Arguments:
-            ncore {int} -- [Specifies the number of threads to use during analysis.] (default: {1})
-            genome {str} -- [The genome that is used for the gene annotation and the enhancer location.] (default: {"hg38"})
-            gene_bed {[type]} -- [Gene annotation for the genome specified with -g as a 12 column BED file.] (default: {None})
-            include_promoter {bool} -- [Include or exclude promoter peaks (<= TSS +/- 2kb) in network inference.] (default: {False})
-            include_enhancer {bool} -- [Include or exclude enhancer peaks (> TSS +/- 2kb) in network inference.] (default: {True})
+        """
+        infer cell type-specific gene regulatory network
 
-        Raises:
-            TypeError: [description]
+        Parameters
+        ----------
+            ncore : int
+                Specifies the number of threads to use during analysis. (default: 1)
+            genome : str
+                The genome that is used for the gene annotation and the enhancer location. (default: "hg38")
+            gene_bed : str, optional
+                Gene annotation for the genome specified with -g as a 12 column BED file. (default: None)
+            include_promoter : bool
+                Include or exclude promoter peaks (<= TSS +/- 2kb) in network inference. (default: False)
+            include_enhancer : bool
+                Include or exclude enhancer peaks (> TSS +/- 2kb) in network inference. (default: True)
         """
         self.ncore = ncore
         self.genome = genome
-        g = Genome(self.genome)
-        self.gsize = g.sizes_file
         self._tmp_files = []
 
         # # Motif information file
@@ -72,31 +68,27 @@ class Network(object):
         # self.motifs2factors = self.pfmfile.replace(".pfm", ".motif2factors.txt")
         # self.factortable = self.pfmfile.replace(".pfm", ".factortable.txt")
 
-        package_dir = os.path.dirname(ananse.__file__)
-
         # Gene information file
-        if self.genome == "hg38":
-            if gene_bed is None:
-                self.gene_bed = os.path.join(package_dir, "db", "hg38.genes.bed")
-            else:
-                self.gene_bed = gene_bed
-        elif self.genome == "hg19":
-            if gene_bed is None:
-                self.gene_bed = os.path.join(package_dir, "db", "hg19.genes.bed")
+        self.gene_bed = gene_bed
+        if gene_bed is None:
+            if self.genome in ["hg38", "hg19"]:
+                self.gene_bed = os.path.join(
+                    PACKAGE_DIR, "db", f"{self.genome}.genes.bed"
+                )
             else:
-                self.gene_bed = gene_bed
-        else:
-            if gene_bed is None:
                 raise TypeError("Please provide a gene bed file with -a argument.")
-            else:
-                self.gene_bed = gene_bed
+        if not os.path.exists(self.gene_bed):
+            raise FileNotFoundError(
+                f"Could not find the gene bed file {self.gene_bed}."
+            )
 
         # self.promoter = promoter
         self.include_promoter = include_promoter
 
         self.include_enhancer = include_enhancer
 
-    def unique_enhancers(self, fname, chrom=None):
+    @staticmethod
+    def unique_enhancers(fname):
         """Extract a list of unique enhancers.
 
         Parameters
@@ -104,14 +96,10 @@ class Network(object):
         fname : str
             File name of a tab-separated file that contains an 'enhancer' column.
 
-        chrom : str, optional
-            Only return enhancers on this chromosome.
-
         Returns
         -------
             PyRanges object with enhancers
         """
-        p = re.compile("[:-]")
         logger.info("reading enhancers")
 
         # Read enhancers from binding file
@@ -148,8 +136,8 @@ class Network(object):
         )
         return enhancers
 
+    @staticmethod
     def distance_weight(
-        self,
         include_promoter=False,
         include_enhancer=True,
         alpha=1e4,
@@ -174,7 +162,7 @@ class Network(object):
 
         Parameters
         ----------
-        include_promoer : bool, optional
+        include_promoter : bool, optional
             Include promoter regions. Default is False.
         include_enhancer : bool, optional
             Include enhancer regions, ie. regions that are distal to the
@@ -207,7 +195,7 @@ class Network(object):
 
         weight1 = pd.DataFrame(
             {
-                "weight": [promoter_weight for z in range(0, promoter_region + 1)],
+                "weight": [promoter_weight for _ in range(0, promoter_region + 1)],
                 "dist": range(0, promoter_region + 1),
             }
         )
@@ -216,7 +204,7 @@ class Network(object):
             {
                 "weight": [
                     enhancer_weight
-                    for z in range(promoter_region + 1, full_weight_region + 1)
+                    for _ in range(promoter_region + 1, full_weight_region + 1)
                 ],
                 "dist": range(promoter_region + 1, full_weight_region + 1),
             }
@@ -430,7 +418,7 @@ class Network(object):
         return fname
 
     def create_expression_network(
-        self, fin_expression, column="tpm", tfs=None, rank=True, bindingfile=None
+        self, fin_expression, column="tpm", tfs=None, factor_activity_file=None
     ):
         """Create a gene expression based network.
 
@@ -442,8 +430,8 @@ class Network(object):
         Parameters
         ----------
         fin_expression : str or list
-            Filename of file that contains gene expression data (TPM), or a
-            list of filenames. First column should contain the gene name.
+            One of more files that contains gene expression data.
+            First column should contain the gene names in HGNC symbols.
 
         column : str, optional
             Column name that contains gene expression, 'tpm' by default (case insensitive).
@@ -451,17 +439,14 @@ class Network(object):
         tfs : list, optional
             List of TF gene names. All TFs will be used by default.
 
-        rank : bool, optional
-            Rank expression levels before scaling.
-
-        bindingfile : str, optional
-            Filename with binding information.
+        factor_activity_file : str, optional
+            factor_activity.tsv created by ANANSE binding. Contains TF names.
 
         Returns
         -------
             Dask DataFrame with gene expression based values.
         """
-        # Convert it to a list if it's not a list of files, but a single file name
+        # Convert to a list of filename(s)
         if isinstance(fin_expression, str):
             fin_expression = [fin_expression]
 
@@ -480,17 +465,11 @@ class Network(object):
         expression[column] = np.log2(expression[column] + 1e-5)
 
         # Create the TF list, based on valid transcription factors
-        if bindingfile is None:
-            bindingfile = "/na/"
         if tfs is None:
-            activity_fname = bindingfile.replace("binding.tsv", "factor_activity.tsv")
-            if os.path.exists(activity_fname):
-                tfs = list(
-                    set(pd.read_table(activity_fname, index_col=0).index.tolist())
-                )
+            if factor_activity_file and os.path.exists(factor_activity_file):
+                tfs = list(set(pd.read_table(factor_activity_file, index_col=0).index))
             else:
-                package_dir = os.path.dirname(ananse.__file__)
-                tffile = os.path.join(package_dir, "db", "tfs.txt")
+                tffile = os.path.join(PACKAGE_DIR, "db", "tfs.txt")
                 tfs = pd.read_csv(tffile, header=None)[0].tolist()
 
         # Save TFs and targets as temporary files
@@ -539,7 +518,6 @@ class Network(object):
         binding,
         fin_expression=None,
         tfs=None,
-        corrfiles=None,
         outfile=None,
         up=1e5,
         down=1e5,
@@ -559,10 +537,8 @@ class Network(object):
         tfs : list, optional
             List of transcription factors to use, by default None, which means
             all TFs will be used.
-        corrfiles : [type], optional
-            Correlation files by default None. CURRENTLY UNUSED.
         outfile : str, optional
-            Output file.
+            Output file. If None, returns a dataframe.
         up : int, optional
             Upstream maximum distance, by default 100kb.
         down : int, optional
@@ -575,8 +551,9 @@ class Network(object):
             Region that will receive full weight, by default 5000."""
         # Expression base network
         logger.info("Loading expression")
+        activity_fname = os.path.join(os.path.dirname(binding), "factor_activity.tsv")
         df_expression = self.create_expression_network(
-            fin_expression, tfs=tfs, rank=True, bindingfile=binding
+            fin_expression, tfs=tfs, factor_activity_file=activity_fname
         )
 
         # Use a version of the binding network, either promoter-based, enhancer-based
@@ -594,7 +571,6 @@ class Network(object):
                 combine_function="sum",
             )
 
-            activity_fname = binding.replace("binding.tsv", "factor_activity.tsv")
             if os.path.exists(activity_fname):
                 logger.info("Reading factor activity")
                 act = pd.read_table(activity_fname, index_col=0)
@@ -636,11 +612,12 @@ class Network(object):
             result["prob"] = result[["tf_expression", "target_expression"]].mean(1)
             result = result.compute()
 
-        logger.info("Writing network")
-        dirname = os.path.dirname(outfile)
-        if dirname:
-            os.makedirs(dirname, exist_ok=True)
-        result[["tf_target", "prob"]].to_csv(outfile, sep="\t", index=False)
+        if outfile:
+            logger.info("Writing network")
+            os.makedirs(os.path.dirname(outfile), exist_ok=True)
+            result[["tf_target", "prob"]].to_csv(outfile, sep="\t", index=False)
+        else:
+            return result[["tf_target", "prob"]]
 
     def __del__(self):
         if not hasattr(self, "_tmp_files"):
@@ -649,3 +626,53 @@ class Network(object):
         for fname in self._tmp_files:
             if os.path.exists(fname):
                 os.unlink(fname)
+
+
+def region_gene_overlap(
+    region_pr,
+    gene_bed,
+    up=100_000,
+    down=100_000,
+):
+    """
+    Couple enhancers to genes.
+
+    Parameters
+    ----------
+    region_pr : PyRanges object
+        PyRanges object with enhancer regions.
+    gene_bed : str
+        gene_bed
+    up : int, optional
+        Upstream maximum distance, by default 100kb.
+    down : int, optional
+        Upstream maximum distance, by default 100kb.
+
+    Returns
+    -------
+    pandas.DataFrame
+        DataFrame with enhancer regions, gene names, distance and weight.
+    """
+    genes = pr.read_bed(gene_bed)
+    # Convert to DataFrame & we don't need intron/exon information
+    genes = genes.as_df().iloc[:, :6]
+
+    # Get the TSS only
+    genes.loc[genes["Strand"] == "+", "End"] = genes.loc[
+        genes["Strand"] == "+", "Start"
+    ]
+    genes.loc[genes["Strand"] == "-", "Start"] = genes.loc[
+        genes["Strand"] == "-", "End"
+    ]
+
+    # Extend up and down
+    genes.loc[genes["Strand"] == "+", "Start"] -= up
+    genes.loc[genes["Strand"] == "+", "End"] += down
+    genes.loc[genes["Strand"] == "-", "Start"] -= down
+    genes.loc[genes["Strand"] == "-", "End"] += up
+
+    # Perform the overlap
+    genes = pr.PyRanges(genes)
+    genes = genes.join(region_pr).as_df()
+
+    return genes
diff --git a/ananse/utils.py b/ananse/utils.py
index f64e470..71a6783 100644
--- a/ananse/utils.py
+++ b/ananse/utils.py
@@ -10,7 +10,6 @@ import warnings
 import genomepy.utils
 from pybedtools import BedTool
 import pysam
-import pyranges as pr
 
 
 def shhh_bedtool(func):
@@ -71,7 +70,7 @@ def samc(ncore):
 def bam_index(bam, force=True, ncore=1):
     if force or not os.path.exists(f"{bam}.bai"):
         index_parameters = [f"-@ {samc(ncore)}", bam]
-        pysam.index(*index_parameters)
+        pysam.index(*index_parameters)  # noqa
 
 
 def bam_sort(bam, ncore=1):
@@ -246,50 +245,3 @@ def check_input_factors(factors):
 
     factors = [line.strip() for line in open(fname)]
     return factors
-
-
-def region_gene_overlap(
-    region_pr,
-    gene_bed,
-    up=100_000,
-    down=100_000,
-):
-    """Couple enhancers to genes.
-
-    Parameters
-    ----------
-    pr : PyRanges object
-        PyRanges object with enhancer regions.
-    up : int, optional
-        Upstream maximum distance, by default 100kb.
-    down : int, optional
-        Upstream maximum distabce, by default 100kb.
-
-    Returns
-    -------
-    pandas.DataFrame
-        DataFrame with enhancer regions, gene names, distance and weight.
-    """
-    genes = pr.read_bed(gene_bed)
-    # Convert to DataFrame & we don't need intron/exon information
-    genes = genes.as_df().iloc[:, :6]
-
-    # Get the TSS only
-    genes.loc[genes["Strand"] == "+", "End"] = genes.loc[
-        genes["Strand"] == "+", "Start"
-    ]
-    genes.loc[genes["Strand"] == "-", "Start"] = genes.loc[
-        genes["Strand"] == "-", "End"
-    ]
-
-    # Extend up and down
-    genes.loc[genes["Strand"] == "+", "Start"] -= up
-    genes.loc[genes["Strand"] == "+", "End"] += down
-    genes.loc[genes["Strand"] == "-", "Start"] -= down
-    genes.loc[genes["Strand"] == "-", "End"] += up
-
-    # Perform the overlap
-    genes = pr.PyRanges(genes)
-    genes = genes.join(region_pr).as_df()
-
-    return genes
diff --git a/requirements.yaml b/requirements.yaml
index 80d395f..6c30812 100644
--- a/requirements.yaml
+++ b/requirements.yaml
@@ -31,7 +31,6 @@ dependencies:
   - scikit-learn
   - scipy >=0.9.0
   - seaborn
-  - six
   - sklearn-contrib-lightning
   - statsmodels
   - tqdm
diff --git a/setup.py b/setup.py
index cbfc930..51b872a 100644
--- a/setup.py
+++ b/setup.py
@@ -43,7 +43,11 @@ setup(
         "gimmemotifs >=0.15.1",
         "loguru",
         "networkx",
+        "numpy",
         "openpyxl",
+        "pandas",
+        "scipy",
+        "scikit-learn",
         "tables",
         "genomepy >= 0.9.3",
         "pyranges",


Test patch:
diff --git a/tests/continuous_integration/test_05_network.py b/tests/continuous_integration/test_05_network.py
index 4acdaf6..2feb15a 100644
--- a/tests/continuous_integration/test_05_network.py
+++ b/tests/continuous_integration/test_05_network.py
@@ -1,4 +1,3 @@
-import os
 from collections import namedtuple
 from tempfile import NamedTemporaryFile
 
@@ -8,7 +7,6 @@ import pandas as pd
 
 from ananse.network import Network
 from ananse.commands import network
-from .test_02_utils import write_file
 
 
 @pytest.fixture
@@ -18,11 +16,7 @@ def binding_fname():
 
 @pytest.fixture
 def network_obj():
-    genome = "tests/data/genome.fa"
-    if not os.path.exists(genome):
-        write_file(genome, [">chr1", "N"])
-
-    return Network(genome=genome, gene_bed="ananse/db/hg38.genes.bed")
+    return Network(genome="", gene_bed="ananse/db/hg38.genes.bed")
 
 
 def test_unique_enhancer(network_obj, binding_fname):
@@ -76,19 +70,17 @@ def test_command():
             "genome annotation include_promoter include_enhancer binding fin_expression outfile ncore",
         )
         args = Args(
-            "hg38",
-            None,
-            True,
-            True,
-            "tests/data/network/binding.tsv.gz",
-            "tests/data/network/heart_expression.txt",
-            fname,
-            2,
+            genome="hg38",
+            annotation=None,
+            include_promoter=True,
+            include_enhancer=True,
+            binding="tests/data/network/binding.tsv.gz",
+            fin_expression="tests/data/network/heart_expression.txt",
+            outfile=fname,
+            ncore=2,
         )
         network(args)
 
-        df = pd.read_table(fname)
-        assert df.shape[0] == 30690
-        assert df.shape[1] == 2
-        assert df.columns[0] == "tf_target"
-        assert df.columns[1] == "prob"
+        df = pd.read_table(fname, sep="\t")
+        assert df.shape[0] == 68820  # 30690
+        assert list(df.columns).__eq__(["tf_target", "prob"])