Skip to content
Snippets Groups Projects
Commit 9239ba7e authored by Thomas Faraut's avatar Thomas Faraut
Browse files

delly-ready pipeline

parent b23f51b5
No related branches found
No related tags found
No related merge requests found
......@@ -155,7 +155,7 @@ rule parsing:
unpack(snakemake_utils.getToolOuputFile),
reference = REFERENCE
output:
parseoutput = "{tool}/{tool}_{chrom}_{svtype}_parsed.vcf.gz"
parseoutput = "parse/{tool}/{tool}_{chrom}_{svtype}_parsed.vcf.gz"
params:
filterbed = FILTERBED
threads:
......
......@@ -10,54 +10,62 @@ from subprocess import run
import pysam
from extractDiscSplitReads import extractDiscSplitReads
from insert_size_distro import insertsizes
def eprint(*args, **kwargs):
''' A simple to sys.stderr printer wrapper '''
print(*args, file=sys.stderr, **kwargs)
def abs_file_path(f):
def existing_file_abs_path(f):
try:
my_real_path = Path(f).resolve()
Path(f).resolve()
except FileNotFoundError:
print("File %s is not an existing and valid file" % f)
exit(1)
else:
my_abs_path = os.path.abspath(f)
return my_abs_path
def file_abs_path(f):
return os.path.abspath(f)
def get_bamfiles(bamlist):
with open(bamlist, "r") as fin:
bamfiles = [abs_file_path(f.strip()) for f in fin]
bamfiles = [existing_file_abs_path(f.strip()) for f in fin]
return bamfiles
def get_all_chromosomes(bamfile):
with pysam.AlignmentFile(bamfile, "rb") as samfile:
header = samfile.header
chromosomes = [contig['SN'] for contig in samfile.header['SQ']]
return chromosomes
def get_unwanted_chroms(chroms_all, chromosome):
if chromosome not in chroms_all:
eprint("Selected chromosome is not a valid chromosome")
chroms_all.remove(chromosome)
return chroms_all
def append_chroms_to_exclude(exclusion_file, unwanted_chromosomes, odir="."):
new_exclusion = os.path.join(odir,"excluded.tsv")
new_exclusion = os.path.join(odir, "excluded.tsv")
copyfile(exclusion_file, new_exclusion)
with open(new_exclusion, "a") as fout:
for chrom in unwanted_chromosomes:
fout.write(chrom+"\n")
return new_exclusion
def runDelly(args):
bamlist = args.bamlist
genome = args.genome
excluded = abs_file_path(args.excluded)
chrom = args.chrom
svtype = args.svtype
output = args.output
excluded = existing_file_abs_path(args.excluded)
output = file_abs_path(args.output)
bamfiles = get_bamfiles(bamlist)
......@@ -85,17 +93,16 @@ def runDelly(args):
# construct delly string
delly_filter_str = "delly filter "
delly_filter_str += "-f germline "
delly_filter_str += "-m 100 "
delly_filter_str += "-t %s " % svtype
delly_filter_str += "-o %s " % output
delly_filter_str += "delly_raw.bcf "
completed = run(delly_filter_str, shell=True)
#delly filter -f germline -m 100 -t ${svtype} -o ${outputfile} delly_raw.bcf
os.chdir(oldir)
#rmtree(tempdir)
rmtree(tempdir)
#return completed.returncode
return completed.returncode
def parse_arguments():
parser = argparse.ArgumentParser(
......
......@@ -61,7 +61,7 @@ class SnakemakeUtils:
inputs = []
for tool in self.tools:
if wildcards.svtype in self.tools_sv_types[tool]:
inputs.append("{tool}/{tool}_{chrom}_{svtype}_parsed.vcf.gz".format(
inputs.append("parse/{tool}/{tool}_{chrom}_{svtype}_parsed.vcf.gz".format(
tool=tool, chrom=wildcards.chrom, svtype=wildcards.svtype
))
return inputs
......@@ -2,7 +2,6 @@ lumpy_sv = ['DEL', 'INV', 'DUP']
rule lumpy:
input:
bams = expand("extraction/{{chrom}}/{sample}.bam", sample=samples),
bamlist = "bamlist.txt",
excluded = "exclusion/excluded.bed"
output:
......@@ -16,4 +15,4 @@ rule lumpy:
stderr = "logs/lumpy/{chrom}.e"
shell:
"lumpy.py -b {input.bamlist} -c {wildcards.chrom}"
" -r {params.readlength -t {threads} | bgzip -c > {output}"
" -r {params.readlength} -t {threads} | bgzip -c > {output}"
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment