diff --git a/README.md b/README.md index 9ac94ea..ccd35cb 100644 --- a/README.md +++ b/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 diff --git a/environment.yml b/environment.yml index e709b16..71e966a 100644 --- a/environment.yml +++ b/environment.yml @@ -6,6 +6,7 @@ channels: dependencies: - click - pandas + - bedtools - pybedtools - dask diff --git a/negative_training_sampler/bed_gc_calculator.py b/negative_training_sampler/bed_gc_calculator.py index 22c3db2..f51cb67 100644 --- a/negative_training_sampler/bed_gc_calculator.py +++ b/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) diff --git a/negative_training_sampler/cli.py b/negative_training_sampler/cli.py index 7e94438..9957a43 100644 --- a/negative_training_sampler/cli.py +++ b/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, diff --git a/negative_training_sampler/fasta_gc_calculator.py b/negative_training_sampler/fasta_gc_calculator.py new file mode 100644 index 0000000..cd4f8a8 --- /dev/null +++ b/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() diff --git a/negative_training_sampler/negative_training_balancer.py b/negative_training_sampler/negative_training_balancer.py index 06ee094..d68f6c6 100644 --- a/negative_training_sampler/negative_training_balancer.py +++ b/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() diff --git a/negative_training_sampler/utils.py b/negative_training_sampler/utils.py index 34f4d24..819767b 100644 --- a/negative_training_sampler/utils.py +++ b/negative_training_sampler/utils.py @@ -3,9 +3,10 @@ """Utility functions""" import pandas as pd +import dask.dataframe as dd -def combine_samples(positive_sample_cleaned, negative_sample_cleaned): +def combine_samples(positive_sample_cleaned, negative_sample_cleaned, sort=True): """Combines positive an negative samples into on dataframe. Arguments: @@ -16,5 +17,23 @@ def combine_samples(positive_sample_cleaned, negative_sample_cleaned): [dataframe] -- Dataframe containing positive and negative labeled samples """ - return pd.concat([positive_sample_cleaned, negative_sample_cleaned], - axis=0).sort_values(by=["chrom", "chromStart"]) + if sort: + return pd.concat([positive_sample_cleaned, negative_sample_cleaned], + axis=0).sort_values(by=["chrom", "chromStart"]) + else: + return pd.concat([positive_sample_cleaned, negative_sample_cleaned], + axis=0) + + +def combine_dataframe(sample_1, sample_2): + """Combines two dask dataframes. + + Arguments: + sample_1 {dask dataframe} -- Dataframe containing sample 1 + sample_2 {dask dataframe} -- Dataframe containing sample 2 + + Returns: + [dataframe] -- Dask Dataframe containing sample 1 and sample 2 + """ + + return dd.concat([sample_1,sample_2]) \ No newline at end of file