Commit 90670433 authored by Beatrice Tan's avatar Beatrice Tan

Added circos zoomed plot.

parent 35d2dcf8
......@@ -14,7 +14,7 @@ include: "rules/GenePrioritization.smk"
#Rules to compare different inputs.
include: "rules/ComparisonRegions.smk"
include: "rules/Circos.yaml"
include: "rules/Circos.smk"
include: "rules/SampleSizes.smk"
include: "rules/UseControl.smk"
......@@ -42,8 +42,9 @@ onerror:
rule all:
"""Define desired output from pipeline."""
input:
os.path.join(config["gisticdir"], "gistic2")
#"Reports/Results.html"
"Samplesize/Report.txt"
#"Reports/Overlap_plots/12.png"
#expand("Reports/Overlap_plots/{region}.png", region=range(22))
rule help:
......
from Circos import InputCircos
from Circos import InputCircos, bed_to_circos
rule circos_input:
"""Make input files for making a circos diagram."""
......@@ -60,3 +60,31 @@ rule add_legend_circos:
offset = ((c_w - l_w), 0)
circos.paste(legend, offset)
circos.save(output[0])
rule circos_input_zoom:
"""Make input files for making a circos diagram."""
input:
bed="Reports/Overlap_known_genes.bed"
output:
gistic="Reports/Circos/Zoom/GISTIC.txt",
rubic="Reports/Circos/Zoom/RUBIC.txt",
genes="Reports/Circos/Zoom/Genes.txt",
run:
bed_to_circos(input.bed, output.rubic, output.gistic, output.genes)
rule make_circos_zoom:
"""Compare locations of known genes, recurrent regions from RUBIC and recurrent regions from GISTIC2."""
input:
gistic="Reports/Circos/Zoom/GISTIC.txt",
rubic="Reports/Circos/Zoom/RUBIC.txt",
genes="Reports/Circos/Zoom/Genes.txt",
output:
plots="Reports/Overlap_plots/12.png"
params:
workflow.basedir + "/scripts/circos/circos_zoom.conf",
chrom='hs12'
conda:
workflow.basedir + "/envs/circos.yaml"
shell:
"circos -conf {params[0]} -outputfile {output[0]} -param gistic_file={input.gistic} -param rubic_file={input.rubic} \
-param gene_file={input.genes} -param chrom={params.chrom}"
......@@ -55,24 +55,3 @@ rule bed_known_genes:
workflow.basedir + "/envs/bedtools.yaml"
shell:
"bedtools intersect -a {input.known} -b {input.gistic} {input.rubic} -names GISTIC RUBIC -wo > {output}"
def get_regions(bed_file):
plot_names = []
with open(bed_file, 'r') as bed:
bed.readline()
for line in bed:
chrom, start = line.split("\t")[0:2]
plot_names.append(chrom + "." + start)
return plot_names
#rule compare_regions:
# """Compare locations of known genes, recurrent regions from RUBIC and recurrent regions from GISTIC2."""
# input:
# overlap="Reports/Overlap_regions.bed"
# params:
# known=os.path.join(workflow.basedir, config["prev_found_genes"]) if config["prev_found_genes"].startswith("input_files") else config["prev_found_genes"]
# output:
# plots=expand("Reports/Overlap_plots/{region}.png", region=get_regions(input.overlap))
# shell:
# "R {workflow.basedir}/scripts/plot_regions.R"
......@@ -14,7 +14,7 @@ rule seg_subsets:
rule run_gistic_subsets:
"""Run GISTIC2 for the segmentation files with different subsets."""
input:
gistic_directory=config["gisticdir"],
gistic_directory=os.path.join(config["gisticdir"], "gistic2"),
seg="Samplesize/Input/Size{rand_nr}.Rep{rep_nr}.txt"
output:
"Samplesize/GISTIC/Size{rand_nr}.Rep{rep_nr}/all_lesions.conf_" + config["gistic_precision"] + ".txt"
......@@ -50,9 +50,9 @@ rule run_rubic_subsets:
rule report_sizes:
"""Report the difference when using different sample sizes."""
input:
#gistic=expand("Samplesize/GISTIC/Size{rand_nr}.Rep{rep_nr}/all_lesions.conf_" + config["gistic_precision"] + ".txt", rand_nr=config["sizes"], rep_nr=config["repeats"]),
rubic_gains=expand("Samplesize/RUBIC/Size{rand_nr}.Rep{rep_nr}/gains.txt", rand_nr=config["sizes"], rep_nr=config["repeats"]),
rubic_losses=expand("Samplesize/RUBIC/Size{rand_nr}.Rep{rep_nr}/losses.txt", rand_nr=config["sizes"], rep_nr=config["repeats"])
gistic=expand("Samplesize/GISTIC/Size{rand_nr}.Rep{rep_nr}/all_lesions.conf_" + config["gistic_precision"] + ".txt", rand_nr=config["sizes"], rep_nr=config["repeats"]),
#rubic_gains=expand("Samplesize/RUBIC/Size{rand_nr}.Rep{rep_nr}/gains.txt", rand_nr=config["sizes"], rep_nr=config["repeats"]),
#rubic_losses=expand("Samplesize/RUBIC/Size{rand_nr}.Rep{rep_nr}/losses.txt", rand_nr=config["sizes"], rep_nr=config["repeats"])
output:
report="Samplesize/Report.txt",
plots="Samplesize/Plots/"
......
......@@ -45,3 +45,19 @@ class InputCircos:
qval = str(max([float(left_q), float(right_q)]))
out_line = ["hs" + chrom, start, end, qval]
out.write(" ".join(out_line) + "\n")
def bed_to_circos(bed_file, rubic_file, gistic_file, gene_file):
rubic = open(rubic_file, 'w')
gistic = open(gistic_file, 'w')
genes = open(gene_file, 'w')
with open(bed_file, 'r') as bed:
for line in bed:
gene_chrom, gene_start, gene_end, gene_name, tool_name, tool_chrom, tool_start, tool_end, \
amp_name, overlap_bp = line.split("\t")
chrom = 'hs' + gene_chrom.strip("chr")
genes.write(" ".join([chrom, gene_start, gene_end]) + "\n")
if tool_name == 'GISTIC':
gistic.write(" ".join([chrom, tool_start, tool_end]) + "\n")
else:
rubic.write(" ".join([chrom, tool_start, tool_end]) + "\n")
rubic.close(), gistic.close(), genes.close()
<<include etc/colors_fonts_patterns.conf>>
<<include ideogram.conf>>
<<include ticks.conf>>
<image>
<<include etc/image.conf>>
</image>
karyotype = data/karyotype/karyotype.human.txt
chromosomes_units = 1000000
chromosomes = conf(chrom)
chromosomes_display_default = no
<highlights>
z = 5
<highlight>
file = conf(gistic_file)
r0 = 0.9r
r1 = 0.7r
fill_color = lgrey
</highlight>
<highlight>
file = conf(rubic_file)
r0 = 0.7r
r1 = 0.5r
fill_color = lyellow
</highlight>
<highlight>
file = conf(gene_file)
r0 = 1.0r
r1 = 0.9r
fill_color = black
</highlight>
</highlights>
<<include etc/housekeeping.conf>>
......@@ -8,9 +8,9 @@ from snakemake.shell import shell
#Convert file locations to absolute paths
segments = os.path.abspath(snakemake.input.seg)
gistic_dir = os.path.abspath(snakemake.input.gistic_directory).split("gistic2")[0]
gistic_dir = os.path.abspath(snakemake.input.gistic_directory).split("/gistic2")[0]
outfolder = os.path.abspath(snakemake.output[0]).split("all_lesions")[0]
print(outfolder)
#Select reference file
ref = snakemake.params.get("ref", "")
ref_file = snakemake.params.get("ref_file", "")
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment