From 35cab6fe81bc3518a493153fdfef4c58fb66f4ec Mon Sep 17 00:00:00 2001 From: Max Schubach Date: Tue, 15 Nov 2022 10:49:28 +0100 Subject: [PATCH 1/7] initial commit --- negative_training_sampler/cli.py | 12 +++- .../fasta_gc_calculator.py | 62 +++++++++++++++++++ .../negative_training_balancer.py | 34 +++++++--- negative_training_sampler/utils.py | 10 ++- 4 files changed, 104 insertions(+), 14 deletions(-) create mode 100644 negative_training_sampler/fasta_gc_calculator.py 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..620d8d5 --- /dev/null +++ b/negative_training_sampler/fasta_gc_calculator.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python + +"""calculates gc content for coordinates in a .bed file""" + +import dask.dataframe as dd +import pybedtools + +BEDCOLS = ["chrom", "chromStart", "chromEnd", + "name", "score", "strand", + "thickStart", "thickEnd", "itemRGB", + "blockCount", "blockSizes", "blockStarts"] + + +def generate_colnames(df, 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 = [] + for i in range(df.columns.str.contains("usercol").sum()-label_num): + colnames.append(BEDCOLS[i]) + for i in range(label_num): + colnames.append("label_{}".format(i+1)) + colnames.append("gc") + colnames.append("num_N") + return colnames + +def get_gc(fasta_file, 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] + """ + bed_df = pybedtools.BedTool(label_file) + bed_gc = bed_df.nuc(reference_file) + gc_df = dd.read_table(bed_gc.fn) + gc_df = gc_df.loc[:, gc_df.columns.str.contains("usercol|gc|num_N")] + colnames = generate_colnames(gc_df, label_num) + gc_df.columns = colnames + 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) + return gc_df + + +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..bc60850 100644 --- a/negative_training_sampler/negative_training_balancer.py +++ b/negative_training_sampler/negative_training_balancer.py @@ -13,6 +13,7 @@ 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 @@ -20,6 +21,7 @@ def balance_trainingdata(label_file, reference_file, genome_file, output_file, + fasta_file, precision, label_num, bgzip, @@ -37,6 +39,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] @@ -61,29 +64,38 @@ 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 # pylint: disable=pointless-statement logging.info("---------------------\ncalculating GC content...\n---------------------") cl_gc = get_gc(label_file, reference_file, label_num, precision) + if useFastaAsPositive(): + positive_sample = get_fasta_gc(fasta_file, precision) + cl_gc = combine_samples(cl_gc, positive_sample, sort=False)) logging.info("---------------------\nextracting positive samples...\n---------------------") - positive_sample = get_positive(cl_gc) + if not useFastaAsPositive(): + positive_sample = get_positive(cl_gc) 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 = (cl_gc.apply(get_negative, seed, meta=dts)).compute() + else: + negative_sample = (cl_gc.groupby(["chrom"], group_keys=False).apply(get_negative, + seed, + meta=dts)).compute() logging.info("---------------------\nloading contigs...\n---------------------") @@ -93,14 +105,18 @@ def balance_trainingdata(label_file, logging.info("---------------------\ncleaning samples\n---------------------") - positive_sample_cleaned = clean_sample(positive_sample, contigs) + if not useFastaAsPositive(): + 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) + if useFastaAsPositive(): + sample_df = negative_sample_cleaned.sort_values(by=["chrom", "chromStart"]) + else: + 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) diff --git a/negative_training_sampler/utils.py b/negative_training_sampler/utils.py index 34f4d24..5b473c7 100644 --- a/negative_training_sampler/utils.py +++ b/negative_training_sampler/utils.py @@ -5,7 +5,7 @@ import pandas as pd -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 +16,9 @@ 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) From e43d69b48da64245061ff2c2f4cdf8226e0519f2 Mon Sep 17 00:00:00 2001 From: sroener <40714954+sroener@users.noreply.github.com> Date: Tue, 15 Nov 2022 11:59:53 +0100 Subject: [PATCH 2/7] refactor: fix typo and add recommendation for conda env --- README.md | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) 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 From 50b7aca3639f4f07683c09ab6976e4bf759473f9 Mon Sep 17 00:00:00 2001 From: sroener <40714954+sroener@users.noreply.github.com> Date: Tue, 15 Nov 2022 12:02:14 +0100 Subject: [PATCH 3/7] feat: add bedtools explicitly for using newer versions --- environment.yml | 1 + 1 file changed, 1 insertion(+) 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 From 56f6012b611e348c0b397985b876e0e062ce3742 Mon Sep 17 00:00:00 2001 From: sroener <40714954+sroener@users.noreply.github.com> Date: Tue, 15 Nov 2022 12:05:43 +0100 Subject: [PATCH 4/7] docs: add comments, improve documentation of get_gc --- .../bed_gc_calculator.py | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) 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) From d20825b4d0d4711365e7423506018684d2964749 Mon Sep 17 00:00:00 2001 From: sroener <40714954+sroener@users.noreply.github.com> Date: Tue, 15 Nov 2022 12:10:09 +0100 Subject: [PATCH 5/7] refactor: apply black code style --- .../negative_training_balancer.py | 68 +++++++++++-------- 1 file changed, 41 insertions(+), 27 deletions(-) diff --git a/negative_training_sampler/negative_training_balancer.py b/negative_training_sampler/negative_training_balancer.py index bc60850..d9bf9f6 100644 --- a/negative_training_sampler/negative_training_balancer.py +++ b/negative_training_sampler/negative_training_balancer.py @@ -17,19 +17,21 @@ from negative_training_sampler.utils import combine_samples -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'): +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. @@ -53,7 +55,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" @@ -69,39 +71,49 @@ def useFastaAsPositive(): logging.info("---------------------\nstarting workers...\n---------------------") - client = Client(n_workers=cores, - threads_per_worker=1, - memory_limit=memory_per_core, - dashboard_address=None) + 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) if useFastaAsPositive(): positive_sample = get_fasta_gc(fasta_file, precision) - cl_gc = combine_samples(cl_gc, positive_sample, sort=False)) + cl_gc = combine_samples(cl_gc, positive_sample, sort=False) - logging.info("---------------------\nextracting positive samples...\n---------------------") + logging.info( + "---------------------\nextracting positive samples...\n---------------------" + ) if not useFastaAsPositive(): 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) if useFastaAsPositive(): negative_sample = (cl_gc.apply(get_negative, seed, meta=dts)).compute() else: - negative_sample = (cl_gc.groupby(["chrom"], group_keys=False).apply(get_negative, - seed, - meta=dts)).compute() + 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---------------------") @@ -123,6 +135,8 @@ def useFastaAsPositive(): else: write_to_stdout(sample_df, precision) - logging.info("---------------------\nshutting down workers...\n---------------------") + logging.info( + "---------------------\nshutting down workers...\n---------------------" + ) client.close() From d0de11e2347e7965fcfc58d6d4216a8e12c186e0 Mon Sep 17 00:00:00 2001 From: Max Schubach Date: Tue, 15 Nov 2022 15:22:18 +0100 Subject: [PATCH 6/7] feat: implemented fasta gc positive sample --- .../fasta_gc_calculator.py | 69 +++++++++++++------ .../negative_training_balancer.py | 21 +++--- negative_training_sampler/utils.py | 15 ++++ 3 files changed, 73 insertions(+), 32 deletions(-) diff --git a/negative_training_sampler/fasta_gc_calculator.py b/negative_training_sampler/fasta_gc_calculator.py index 620d8d5..35d5432 100644 --- a/negative_training_sampler/fasta_gc_calculator.py +++ b/negative_training_sampler/fasta_gc_calculator.py @@ -3,15 +3,13 @@ """calculates gc content for coordinates in a .bed file""" import dask.dataframe as dd -import pybedtools +import dask.bag as db +import pandas as pd -BEDCOLS = ["chrom", "chromStart", "chromEnd", - "name", "score", "strand", - "thickStart", "thickEnd", "itemRGB", - "blockCount", "blockSizes", "blockStarts"] +COLUMNS = ["chrom", "chromStart", "chromEnd"] -def generate_colnames(df, label_num): +def generate_colnames(label_num): """Generates column names for an input dataframe. Arguments: @@ -21,16 +19,47 @@ def generate_colnames(df, label_num): Returns: [list] -- [list of generated column names] """ - colnames = [] - for i in range(df.columns.str.contains("usercol").sum()-label_num): - colnames.append(BEDCOLS[i]) + colnames = [] + COLUMNS for i in range(label_num): colnames.append("label_{}".format(i+1)) colnames.append("gc") - colnames.append("num_N") return colnames -def get_gc(fasta_file, precision): +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: @@ -41,16 +70,14 @@ def get_gc(fasta_file, precision): Returns: [dataframe] -- [Dataframe containing viable entries and their respective gc content] """ - bed_df = pybedtools.BedTool(label_file) - bed_gc = bed_df.nuc(reference_file) - gc_df = dd.read_table(bed_gc.fn) - gc_df = gc_df.loc[:, gc_df.columns.str.contains("usercol|gc|num_N")] - colnames = generate_colnames(gc_df, label_num) - gc_df.columns = colnames - 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) - return gc_df + colnames = generate_colnames(label_num) + 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(): diff --git a/negative_training_sampler/negative_training_balancer.py b/negative_training_sampler/negative_training_balancer.py index d9bf9f6..d68f6c6 100644 --- a/negative_training_sampler/negative_training_balancer.py +++ b/negative_training_sampler/negative_training_balancer.py @@ -15,6 +15,7 @@ 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 +from negative_training_sampler.utils import combine_dataframe def balance_trainingdata( @@ -85,15 +86,14 @@ def useFastaAsPositive(): cl_gc = get_gc(label_file, reference_file, label_num, precision) if useFastaAsPositive(): - positive_sample = get_fasta_gc(fasta_file, precision) - cl_gc = combine_samples(cl_gc, positive_sample, sort=False) + positive_sample = get_fasta_gc(fasta_file, label_num, precision) + cl_gc = combine_dataframe(cl_gc, positive_sample) logging.info( "---------------------\nextracting positive samples...\n---------------------" ) - - if not useFastaAsPositive(): - positive_sample = get_positive(cl_gc) + + positive_sample = get_positive(cl_gc) logging.info( "---------------------\nbalancing negative sample set...\n---------------------" @@ -101,7 +101,7 @@ def useFastaAsPositive(): dts = dict(cl_gc.dtypes) if useFastaAsPositive(): - negative_sample = (cl_gc.apply(get_negative, seed, meta=dts)).compute() + negative_sample = get_negative(cl_gc.compute(), seed) else: negative_sample = ( cl_gc.groupby(["chrom"], group_keys=False).apply( @@ -117,16 +117,15 @@ def useFastaAsPositive(): logging.info("---------------------\ncleaning samples\n---------------------") - if not useFastaAsPositive(): + 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---------------------") - if useFastaAsPositive(): - sample_df = negative_sample_cleaned.sort_values(by=["chrom", "chromStart"]) - else: - sample_df = combine_samples(positive_sample_cleaned, negative_sample_cleaned) + sample_df = combine_samples(positive_sample_cleaned, negative_sample_cleaned) # print(sample_df.head()) diff --git a/negative_training_sampler/utils.py b/negative_training_sampler/utils.py index 5b473c7..819767b 100644 --- a/negative_training_sampler/utils.py +++ b/negative_training_sampler/utils.py @@ -3,6 +3,7 @@ """Utility functions""" import pandas as pd +import dask.dataframe as dd def combine_samples(positive_sample_cleaned, negative_sample_cleaned, sort=True): @@ -22,3 +23,17 @@ def combine_samples(positive_sample_cleaned, negative_sample_cleaned, sort=True) 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 From d373d00f22b72ffdc0069affe37c32bc95109fbb Mon Sep 17 00:00:00 2001 From: Max Schubach Date: Tue, 15 Nov 2022 15:36:36 +0100 Subject: [PATCH 7/7] adding option for gz --- negative_training_sampler/fasta_gc_calculator.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/negative_training_sampler/fasta_gc_calculator.py b/negative_training_sampler/fasta_gc_calculator.py index 35d5432..cd4f8a8 100644 --- a/negative_training_sampler/fasta_gc_calculator.py +++ b/negative_training_sampler/fasta_gc_calculator.py @@ -5,6 +5,7 @@ import dask.dataframe as dd import dask.bag as db import pandas as pd +import gzip COLUMNS = ["chrom", "chromStart", "chromEnd"] @@ -71,7 +72,11 @@ def get_gc(fasta_file, label_num, precision): [dataframe] -- [Dataframe containing viable entries and their respective gc content] """ colnames = generate_colnames(label_num) - file_handle = open(fasta_file, 'r') + + 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)