feat: use fasta as positive

README.md

@@ -20,6 +20,12 @@ dependencies:
 
 ### Setup
 
+It is recommended to setup a conda environment prior to installing the package. By doing so, one has more control over the bedtools version used by pybedtools.
+
+```bash
+env create -n negative_training_sampler --file environment.yml
+```
+
 Installing the package:
 
 ```bash
@@ -33,13 +39,13 @@ The package needs at least a minimal bed file with positive (1) and negative (0)
 General use:
 
 ```bash
-negative_input_sampler -i LABEL_FILE -r REFERENCE_FILE -g GENOME_FILE -o OUTPUT_FILE
+negative_training_sampler -i LABEL_FILE -r REFERENCE_FILE -g GENOME_FILE -o OUTPUT_FILE
 ```
 
 More advanced use:
 
 ```bash
-negative_input_sampler -i LABEL_FILE -r REFERENCE_FILE -g GENOME_FILE -o OUTPUT_FILE --cores INT --memory [INT]GB
+negative_training_sampler -i LABEL_FILE -r REFERENCE_FILE -g GENOME_FILE -o OUTPUT_FILE --cores INT --memory [INT]GB
 ```
 
 ### help

environment.yml

@@ -6,6 +6,7 @@ channels:
 dependencies:
   - click
   - pandas
+  - bedtools
   - pybedtools
   - dask
 

negative_training_sampler/bed_gc_calculator.py

@@ -40,13 +40,39 @@ def get_gc(label_file, reference_file, label_num, precision):
 
     Returns:
         [dataframe] -- [Dataframe containing viable entries and their respective gc content]
+
+
+    Example: 
+
+    For an input label file like this:
+
+    chr1	191200	191500	0.0
+    chr1	205950	206250	0.0
+    chr1	296600	296900	0.0
+    chr1	354950	355250	0.0
+    chr1	356950	357250	0.0
+
+    The output of the function will look like this:
+
+    chrom  chromStart  chromEnd  label_1         gc
+    0  chr1      191200    191500      0.0  54.669998
+    1  chr1      205950    206250      0.0  42.669998
+    2  chr1      296600    296900      0.0  35.330002
+    3  chr1      354950    355250      0.0  43.669998
+    4  chr1      356950    357250      0.0  49.669998
     """
+    # loads bed file
     bed_df = pybedtools.BedTool(label_file)
+    # calculates nucleotide content -> https://daler.github.io/pybedtools/autodocs/pybedtools.bedtool.BedTool.nucleotide_content.html#pybedtools.bedtool.BedTool.nucleotide_content
     bed_gc = bed_df.nuc(reference_file)
+    # reads output table of nuc as table -> bed_gc.fn contains the filename of the temporary file created by bedtools
     gc_df = dd.read_table(bed_gc.fn)
+    #get usercolumns (all columns provided in the bedfile), GC content and number of Ns
     gc_df = gc_df.loc[:, gc_df.columns.str.contains("usercol|gc|num_N")]
+    # generate column names used later on 
     colnames = generate_colnames(gc_df, label_num)
     gc_df.columns = colnames
+    # drop entries if they contain Ns and drop the column afterwards
     gc_df = gc_df.loc[gc_df.num_N == 0].drop("num_N", axis=1)
     gc_df["gc"] = gc_df["gc"].astype("float32")*100
     gc_df["gc"] = gc_df["gc"].round(precision)

negative_training_sampler/cli.py

@@ -4,6 +4,7 @@
 
 from negative_training_sampler.negative_training_balancer import balance_trainingdata
 
+
 @click.command()
 @click.option("-i",
               "--label-file",
@@ -23,6 +24,12 @@
               required=True,
               type=click.Path(exists=True, readable=True),
               help="Input genome file of reference")
+@click.option("-f",
+              "--fasta-file",
+              'fasta_file',
+              required=False,
+              type=click.Path(exists=True, readable=True),
+              help="Use fasta sequences as positive label insteald of label bed file")
 @click.option("-o",
               "--output_file",
               'output_file',
@@ -60,8 +67,8 @@
 @click.option("--memory",
               default="2GB",
               help="amount of memory per core (e.g. 2 cores * 2GB = 4GB)\ndefault: 2GB")
-def cli(label_file, reference_file, genome_file, output_file, precision, label_num,
-        bgzip, log_file, verbose, seed, cores, memory): # pylint: disable=no-value-for-parameter
+def cli(label_file, reference_file, genome_file, fasta_file, output_file, precision, label_num,
+        bgzip, log_file, verbose, seed, cores, memory):  # pylint: disable=no-value-for-parameter
     '''
     A simple script that takes a tsv file with positive and negative labels
     and a reference file. Generates negative samples with the same GC distribution
@@ -72,6 +79,7 @@ def cli(label_file, reference_file, genome_file, output_file, precision, label_n
                          reference_file=reference_file,
                          genome_file=genome_file,
                          precision=precision,
+                         fasta_file=fasta_file,
                          output_file=output_file,
                          label_num=label_num,
                          bgzip=bgzip,

negative_training_sampler/fasta_gc_calculator.py

@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+
+"""calculates gc content for coordinates in a .bed file"""
+
+import dask.dataframe as dd
+import dask.bag as db
+import pandas as pd
+import gzip
+
+COLUMNS = ["chrom", "chromStart", "chromEnd"]
+
+
+def generate_colnames(label_num):
+    """Generates column names for an input dataframe.
+
+    Arguments:
+        df {dataframe} -- [dataframe for which column names are generated]
+        label_num {int} -- [number of provided labels]
+
+    Returns:
+        [list] -- [list of generated column names]
+    """
+    colnames = [] + COLUMNS
+    for i in range(label_num):
+        colnames.append("label_{}".format(i+1))
+    colnames.append("gc")
+    return colnames
+
+def read_fasta(filehandle):
+    """
+    Reads fasta entries from filehandle
+    supports multiline fasta
+    """
+    count = 0
+    seq = ""
+    id = ""
+    for line in filehandle:
+        line = line.rstrip("\n\r")
+        if ( (count == 0) and line.startswith(">")): # Read identifier
+            id = line[1:]
+            count+=1
+        elif count == 1:        # read sequence
+            seq = line
+            count+=1
+        elif count == 2:  # multiple case:  a) more sequence (fasta), c) next sequence identifier (fasta)
+            if line.startswith(">") : # case c)
+                yield seq
+                id, seq  = None, None
+                count = 1
+                id = line[1:]
+            else: # case b)
+                seq = seq + line
+        else:
+          sys.stderr.write("Unexpected line:" + str(line.strip()) + "\n")
+          count = 0
+    if id and seq:
+        yield seq
+    return
+
+def compute_gc(seq):
+    return (seq.count("G")+seq.count("C")+seq.count("c")+seq.count("g"))/(len(seq)-seq.count("n")-seq.count("N"))
+
+def get_gc(fasta_file, label_num, precision):
+    """Calculates gc content for all viable entries in an input dataframe.
+
+    Arguments:
+        fasta_file {dataframe} -- [Dataframe containing genomic regions labeled
+                                   as positive(1) or negative(0)]
+        genome_file {str} -- [Path to a refrence genome in FASTA format]
+
+    Returns:
+        [dataframe] -- [Dataframe containing viable entries and their respective gc content]
+    """
+    colnames = generate_colnames(label_num)
+
+    if fasta_file.endswith(".gz"):
+        file_handle = gzip.open(fasta_file, 'rt')
+    else:
+        file_handle = open(fasta_file, 'r')
+    seq_db = db.from_sequence([i for i in read_fasta(file_handle)])
+    seq_dd = dd.from_pandas(pd.DataFrame({"gc":db.map(compute_gc,seq_db).compute(),"label_1":1},columns=colnames), npartitions=1)
+
+    seq_dd["gc"] = seq_dd["gc"].astype("float32")*100
+    seq_dd["gc"] = seq_dd["gc"].round(precision)
+    return seq_dd
+
+
+def main():
+    """[summary]
+    """
+
+
+if __name__ == "__main__":
+    main()

negative_training_sampler/negative_training_balancer.py

@@ -13,21 +13,26 @@
 from negative_training_sampler.io import write_to_file
 from negative_training_sampler.io import write_to_stdout
 from negative_training_sampler.io import load_contigs
+from negative_training_sampler.fasta_gc_calculator import get_gc as get_fasta_gc
 from negative_training_sampler.utils import combine_samples
-
-
-def balance_trainingdata(label_file,
-                         reference_file,
-                         genome_file,
-                         output_file,
-                         precision,
-                         label_num,
-                         bgzip,
-                         log_file,
-                         verbose,
-                         seed,
-                         cores=1,
-                         memory_per_core='2GB'):
+from negative_training_sampler.utils import combine_dataframe
+
+
+def balance_trainingdata(
+    label_file,
+    reference_file,
+    genome_file,
+    output_file,
+    fasta_file,
+    precision,
+    label_num,
+    bgzip,
+    log_file,
+    verbose,
+    seed,
+    cores=1,
+    memory_per_core="2GB",
+):
     """
     Function that calculates the GC content for positive and negative labeled genomic regions and
     balances their number based on GC content per chromosome.
@@ -37,6 +42,7 @@ def balance_trainingdata(label_file,
                                     labeled as positive(1) or negative(0)]
         reference_file {str}    -- [Path to a reference genome in FASTA format]
         genome_file {str}       -- [Path to the genome file of the reference]
+        fast_file {str}         -- [Path to a fasta file to use as positve labels. If not None.]
         output_file {str}       -- [Name of the output file. File will be in .bed format]
         precision {int}         -- [Precision of decimals when computing the attributes like GC content]
         label_num {int}         -- [Number of provided label columns]
@@ -50,7 +56,7 @@ def balance_trainingdata(label_file,
     """
 
     loglevel = logging.INFO
-    logformat = '%(message)s'
+    logformat = "%(message)s"
     if verbose:
         loglevel = logging.DEBUG
         logformat = "%(asctime)s: %(levelname)s - %(message)s"
@@ -61,52 +67,75 @@ def balance_trainingdata(label_file,
     else:
         logging.basicConfig(stream=sys.stderr, level=loglevel, format=logformat)
 
+    def useFastaAsPositive():
+        return fasta_file is not None
+
     logging.info("---------------------\nstarting workers...\n---------------------")
 
-    client = Client(n_workers=cores,
-                    threads_per_worker=1,
-                    memory_limit=memory_per_core,
-                    dashboard_address=None)
-    client # pylint: disable=pointless-statement
+    client = Client(
+        n_workers=cores,
+        threads_per_worker=1,
+        memory_limit=memory_per_core,
+        dashboard_address=None,
+    )
+    client  # pylint: disable=pointless-statement
 
-    logging.info("---------------------\ncalculating GC content...\n---------------------")
+    logging.info(
+        "---------------------\ncalculating GC content...\n---------------------"
+    )
 
     cl_gc = get_gc(label_file, reference_file, label_num, precision)
-
-    logging.info("---------------------\nextracting positive samples...\n---------------------")
-
+    if useFastaAsPositive():
+        positive_sample = get_fasta_gc(fasta_file, label_num, precision)
+        cl_gc = combine_dataframe(cl_gc, positive_sample)
+
+    logging.info(
+        "---------------------\nextracting positive samples...\n---------------------"
+    )
+
     positive_sample = get_positive(cl_gc)
 
-    logging.info("---------------------\nbalancing negative sample set...\n---------------------")
+    logging.info(
+        "---------------------\nbalancing negative sample set...\n---------------------"
+    )
 
     dts = dict(cl_gc.dtypes)
-    negative_sample = (cl_gc.groupby(["chrom"], group_keys=False).apply(get_negative,
-                                                                        seed,
-                                                                        meta=dts)
-                       ).compute()
+    if useFastaAsPositive():
+        negative_sample = get_negative(cl_gc.compute(), seed)
+    else:
+        negative_sample = (
+            cl_gc.groupby(["chrom"], group_keys=False).apply(
+                get_negative, seed, meta=dts
+            )
+        ).compute()
 
     logging.info("---------------------\nloading contigs...\n---------------------")
 
     contigs = load_contigs(genome_file)
 
-#    print(contigs)
+    #    print(contigs)
 
     logging.info("---------------------\ncleaning samples\n---------------------")
 
-    positive_sample_cleaned = clean_sample(positive_sample, contigs)
+    if useFastaAsPositive():
+        positive_sample_cleaned = positive_sample
+    else:
+        positive_sample_cleaned = clean_sample(positive_sample, contigs)
     negative_sample_cleaned = clean_sample(negative_sample, contigs)
 
     logging.info("---------------------\nsaving results\n---------------------")
 
     sample_df = combine_samples(positive_sample_cleaned, negative_sample_cleaned)
 
-    #print(sample_df.head())
+    # print(sample_df.head())
 
     if output_file:
         write_to_file(sample_df, output_file, bgzip)
     else:
         write_to_stdout(sample_df, precision)
 
-    logging.info("---------------------\nshutting down workers...\n---------------------")
+    logging.info(
+        "---------------------\nshutting down workers...\n---------------------"
+    )
 
     client.close()