Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

change cli to remove inputs as options #185

Merged
merged 6 commits into from
Jan 28, 2025
Merged

change cli to remove inputs as options #185

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
488db5d
add compression for output fasta
pdimens Jan 27, 2025
fa1a936
add gz to fasta
pdimens Jan 27, 2025
5674b0b
change cli to remove inputs as options
pdimens Jan 28, 2025
416d78a
update rule to match new cli
pdimens Jan 28, 2025
c8517e0
update summary for new api
pdimens Jan 28, 2025
e91ddab
rm redundant required= from argparse [skip cli]
pdimens Jan 28, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 4 additions & 5 deletions harpy/bin/inline_to_haplotag.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,13 +10,13 @@
parser = argparse.ArgumentParser(
prog = 'inline_to_haplotag.py',
description = 'Moves inline linked read barcodes to read headers (OX:Z) and converts them into haplotag ACBD format (BX:Z). Barcodes must all be the same length.',
usage = f"inline_to_haplotag.py -f <forward.fq.gz> -r <reverse.fq.gz> -b <barcodes.txt> -p <prefix>",
usage = f"inline_to_haplotag.py -b <barcodes.txt> -p <prefix> FORWARD.fq.gz REVERSE.fq.gz",
exit_on_error = False
)
parser.add_argument("-f", "--forward", required = True, type = str, help = "Forward reads of paired-end FASTQ file pair (gzipped)")
parser.add_argument("-r", "--reverse", required = True, type = str, help = "Reverse reads of paired-end FASTQ file pair (gzipped)")
parser.add_argument("-p", "--prefix", required = True, type = str, help = "Prefix for outfile files (e.g. <prefix>.R1.fq.gz)")
parser.add_argument("-b", "--barcodes", required = True, type=str, help="File listing the linked-read barcodes to convert to haplotag format, one barcode per line")
parser.add_argument("-b", "--barcodes", required = True, type=str, help="Barcode conversion key file with format: ATCG<tab>ACBD")
parser.add_argument("forward", required = True, type = str, help = "Forward reads of paired-end FASTQ file pair (gzipped)")
parser.add_argument("reverse", required = True, type = str, help = "Reverse reads of paired-end FASTQ file pair (gzipped)")
if len(sys.argv) == 1:
parser.print_help(sys.stderr)
sys.exit(1)
Expand Down Expand Up @@ -112,7 +112,6 @@ def process_record(fw_entry, rv_entry, barcode_database, bc_len):
sys.exit(1)

insert_key_value(bc_db, ATCG, ACBD)
#bc_dict[ATCG] = ACBD
lengths.add(len(ATCG))
if len(lengths) > 1:
sys.stderr.write("Can only search sequences for barcodes of a single length, but multiple barcode legnths detected: " + ",".join([str(i) for i in lengths]))
Expand Down
4 changes: 2 additions & 2 deletions harpy/snakefiles/simulate_linkedreads.smk
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ rule demultiplex_barcodes:
container:
None
shell:
"inline_to_haplotag.py -f {input.fw} -r {input.rv} -b {input.barcodes} -p {params} 2> {log}"
"inline_to_haplotag.py -b {input.barcodes} -p {params} {input.fw} {input.rv} 2> {log}"

rule workflow_summary:
default_target: True
Expand Down Expand Up @@ -202,7 +202,7 @@ rule workflow_summary:
haplosim += f"\tHaploSim.pl -g genome1,genome2 -a dwgsimreads1,dwgsimreads2 -l {params.lrbc_len} -p {params.lrsproj_dir}/linked_molecules/lrsim -b BARCODES -i {params.lrsoutdist} -s {params.lrsdist_sd} -x {params.lrsn_pairs} -f {params.lrsmol_len} -t {params.lrsparts} -m {params.lrsmols_per} -z THREADS {params.lrsstatic}"
summary.append(haplosim)
bxconvert = "Inline barcodes were converted in haplotag BX:Z tags using:\n"
bxconvert += "\tinline_to_haplotag.py -f <forward.fq.gz> -r <reverse.fq.gz> -b <barcodes.txt> -p <prefix>"
bxconvert += "\tinline_to_haplotag.py -b <barcodes.txt> -p <prefix> forward.fq.gz reverse.fq.gz"
summary.append(bxconvert)
sm = "The Snakemake workflow was called via command line:\n"
sm += f"\t{config['workflow_call']}"
Expand Down
33 changes: 20 additions & 13 deletions harpy/snakefiles/simulate_snpindel.smk
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ rule simulate_haploid:
in_vcfs,
geno = genome
output:
f"{outdir}/{outprefix}.simseq.genome.fa",
temp(f"{outdir}/{outprefix}.simseq.genome.fa"),
f"{outdir}/{outprefix}.refseq2simseq.SNP.vcf" if snp else [],
f"{outdir}/{outprefix}.refseq2simseq.INDEL.vcf" if indel else [],
f"{outdir}/{outprefix}.refseq2simseq.map.txt"
Expand All @@ -115,14 +115,17 @@ rule rename_haploid:
indelvcf = f"{outdir}/{outprefix}.refseq2simseq.INDEL.vcf" if indel else [],
mapfile = f"{outdir}/{outprefix}.refseq2simseq.map.txt"
output:
fasta = f"{outdir}/{outprefix}.fasta",
fasta = f"{outdir}/{outprefix}.fasta.gz",
snpvcf = f"{outdir}/{outprefix}.snp.vcf" if snp else [],
indelvcf = f"{outdir}/{outprefix}.indel.vcf" if indel else [],
mapfile = f"{outdir}/{outprefix}.map"
run:
for i,j in zip(input, output):
if i:
os.rename(i,j)
shell(f"bgzip -c {input.fasta} > {output.fasta}")
os.rename(input.mapfile, output.mapfile)
if input.snpvcf:
os.rename(input.snpvcf, output.snpvcf)
if input.indelvcf:
os.rename(input.indelvcf, output.indelvcf)

rule diploid_snps:
input:
Expand Down Expand Up @@ -160,8 +163,8 @@ rule simulate_diploid:
indel_hap = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.indel.vcf" if indel else [],
geno = genome
output:
f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.simseq.genome.fa",
f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.refseq2simseq.map.txt",
temp(f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.simseq.genome.fa"),
temp(f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.refseq2simseq.map.txt"),
temp(f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.refseq2simseq.INDEL.vcf") if indel else [],
temp(f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.refseq2simseq.SNP.vcf") if snp else []
log:
Expand All @@ -180,18 +183,22 @@ rule rename_diploid:
fasta= f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.simseq.genome.fa",
mapfile = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.refseq2simseq.map.txt",
output:
fasta = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.fasta",
fasta = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.fasta.gz",
mapfile = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.map"
run:
for i,j in zip(input, output):
os.rename(i,j)
container:
None
shell:
"""
bgzip -c {input.fasta} > {output.fasta}
cp {input.mapfile} {output.mapfile}
"""

rule workflow_summary:
default_target: True
input:
f"{outdir}/{outprefix}.fasta",
f"{outdir}/{outprefix}.fasta.gz",
collect(f"{outdir}/{outprefix}" + ".{var}.vcf", var = variants),
collect(f"{outdir}/haplotype_" + "{n}/" + outprefix + ".hap{n}.fasta", n = [1,2]) if heterozygosity > 0 and not only_vcf else [],
collect(f"{outdir}/haplotype_" + "{n}/" + outprefix + ".hap{n}.fasta.gz", n = [1,2]) if heterozygosity > 0 and not only_vcf else [],
collect(f"{outdir}/haplotype_" + "{n}/" + outprefix + ".hap{n}" + ".{var}.vcf", n = [1,2], var = variants) if heterozygosity > 0 else [],
collect(f"{outdir}/haplotype_" + "{n}/" + outprefix + ".hap{n}" + ".map", n = [1,2]) if heterozygosity > 0 else []
params:
Expand Down
29 changes: 17 additions & 12 deletions harpy/snakefiles/simulate_variants.smk
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ rule simulate_haploid:
vcf_correct if vcf else [],
geno = genome
output:
f"{outdir}/{outprefix}.simseq.genome.fa",
temp(f"{outdir}/{outprefix}.simseq.genome.fa"),
f"{outdir}/{outprefix}.refseq2simseq.{simuG_variant}.vcf",
f"{outdir}/{outprefix}.refseq2simseq.map.txt"
log:
Expand All @@ -77,12 +77,13 @@ rule rename_haploid:
vcf = f"{outdir}/{outprefix}.refseq2simseq.{simuG_variant}.vcf",
mapfile = f"{outdir}/{outprefix}.refseq2simseq.map.txt"
output:
fasta = f"{outdir}/{outprefix}.fasta",
fasta = f"{outdir}/{outprefix}.fasta.gz",
vcf = f"{outdir}/{outprefix}.{variant}.vcf",
mapfile = f"{outdir}/{outprefix}.{variant}.map"
run:
for i,j in zip(input, output):
os.rename(i,j)
shell(f"bgzip -c {input.fasta} > {output.fasta}")
os.rename(input.mapfile, output.mapfile)
os.rename(input.vcf, output.vcf)

rule diploid_variants:
input:
Expand Down Expand Up @@ -112,8 +113,8 @@ rule simulate_diploid:
hap = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.{variant}.vcf",
geno = genome
output:
f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.simseq.genome.fa",
f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.refseq2simseq.map.txt",
temp(f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.simseq.genome.fa"),
temp(f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.refseq2simseq.map.txt"),
temp(f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.refseq2simseq.{simuG_variant}.vcf")
log:
f"{outdir}/logs/{outprefix}.hap{{haplotype}}.log"
Expand All @@ -130,18 +131,22 @@ rule rename_diploid:
fasta = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.simseq.genome.fa",
mapfile = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.refseq2simseq.map.txt"
output:
fasta = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.fasta",
fasta = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.fasta.gz",
mapfile = f"{outdir}/haplotype_{{haplotype}}/{outprefix}.hap{{haplotype}}.{variant}.map"
run:
for i,j in zip(input, output):
os.rename(i,j)
container:
None
shell:
"""
bgzip -c {input.fasta} > {output.fasta}
cp {input.mapfile} {output.mapfile}
"""

rule workflow_summary:
default_target: True
input:
f"{outdir}/{outprefix}.fasta",
f"{outdir}/{outprefix}.fasta.gz",
f"{outdir}/{outprefix}.{variant}.vcf",
collect(f"{outdir}/haplotype_" + "{n}" + f"/{outprefix}.hap" + "{n}.fasta", n = [1,2]) if heterozygosity > 0 and not only_vcf else [],
collect(f"{outdir}/haplotype_" + "{n}" + f"/{outprefix}.hap" + "{n}.fasta.gz", n = [1,2]) if heterozygosity > 0 and not only_vcf else [],
collect(f"{outdir}/haplotype_" + "{n}" + f"/{outprefix}.hap" + "{n}" + f".{variant}.vcf", n = [1,2]) if heterozygosity > 0 else [],
collect(f"{outdir}/haplotype_" + "{n}" + f"/{outprefix}.hap" + "{n}" + f".{variant}.map", n = [1,2]) if heterozygosity > 0 else []
params:
Expand Down
Loading