adding sequencing effort limitation by random subsampling

Author: maxibor

adrsm/adrsm.py

@@ -5,6 +5,7 @@
 from adrsm import __version__
 import click
 import csv
+import random
 
 
 @click.command()
@@ -68,6 +69,14 @@
     show_default=True,
     help="Deamination substitution max frequency",
 )
+@click.option(
+    "-e",
+    "--effort",
+    default=100,
+    type=int,
+    show_default=True,
+    help="Sequencing effort, maximum number of reads to be generated",
+)
 @click.option(
     "-s",
     "--seed",
@@ -157,6 +166,7 @@ def main(
     geom_p,
     mind,
     maxd,
+    effort,
     seed,
     threads,
     output,
@@ -167,8 +177,9 @@ def main(
     fastq_list = []
     stat_dict = {}
     all_genomes = read_config(conffile)
+    all_reads = []
     for agenome in all_genomes.keys():
-        stat_and_run = ad.run_read_simulation_multi(
+        reads, stat_and_run = ad.run_read_simulation_multi(
             INFILE=agenome,
             COV=all_genomes[agenome]["cov"],
             READLEN=readlength,
@@ -185,14 +196,27 @@ def main(
             THEMIN=mind,
             THEMAX=maxd,
             PROCESS=threads,
-            FASTQ_OUT=output,
         )
         stat_dict[ad.get_basename(agenome)] = stat_and_run
+        all_reads.extend(reads)
 
+    if len(all_reads) > effort:
+        all_reads = random.sample(all_reads, effort)
+    ad.write_fastq_multi(all_reads, output)
     ad.write_stat(stat_dict=stat_dict, stat_out=stats)
-    print("\n-- ADRSM v" + __version__ + " finished generating this mock metagenome --")
     print(
-        "-- FASTQ files written to " + output + ".1.fastq and " + output + ".2.fastq --"
+        "\n-- ADRSM v"
+        + __version__
+        + " finished generating "
+        + str(len(all_reads))
+        + " reads for this mock metagenome --"
+    )
+    print(
+        "-- FASTQ files written to "
+        + output
+        + ".1.fastq.gz and "
+        + output
+        + ".2.fastq.gz --"
     )
     print("-- Statistic file written to " + stats + " --")
 

adrsm/lib/adrsmlib.py

@@ -188,7 +188,6 @@ def run_read_simulation_multi(
     THEMIN,
     THEMAX,
     PROCESS,
-    FASTQ_OUT,
 ):
     print("===================\n===================")
     print("Genome: ", INFILE)
@@ -269,8 +268,9 @@ def run_read_simulation_multi(
         read_length=READLEN,
         process=PROCESS,
     )
-    write_fastq_multi(fastq_list=result, outputfile=FASTQ_OUT)
-    return [nread * INSERLEN, INSERLEN, COV, DAMAGE]
+
+    # write_fastq_multi(fastq_list=result, outputfile=FASTQ_OUT)
+    return result, [nread * INSERLEN, INSERLEN, COV, DAMAGE]
 
 
 def specie_to_taxid(specie):