diff --git a/config/example.config.yml b/config/example.config.yml index 3f439f4..7bf9f05 100644 --- a/config/example.config.yml +++ b/config/example.config.yml @@ -34,7 +34,13 @@ plotting: smoothing: True # bool; apply Savitzky-Golay filter or not rolling: True # bool; smooth by substracting the rolling median background_norm: True # bool, normalize with random background or not + edge_norm: False # bool, normalize with the mean of both outer 100 bp + win_len: 21 #Savitzky-Golay filter window length (length of the filter window i.e. the number of coefficients) + poly: 2 #Savitzky-Golay filter polyorder (order of the polynomial used to fit the samples) +### Midpoint coverage + +midpoint_coverage: False # bool, compute midpoint coverage score additionally to WPS ## unsupervised diff --git a/snakefile_WPS.smk b/snakefile_WPS.smk index 8491bfc..cd5f671 100644 --- a/snakefile_WPS.smk +++ b/snakefile_WPS.smk @@ -48,6 +48,15 @@ def get_STARTS_ref(sample): GENOME=genomes, ) +def get_MP_ref(sample): + ref_samples = samples["ref_samples"][sample].split(",") + genomes = samples["genome_build"][ref_samples].values.tolist() + return expand( + "results/intermediate/{{ID}}/table/{GENOME}/target/{{target_region}}--{ref_SAMPLE}_MIDPOINT.csv.gz", + zip, + ref_SAMPLE=ref_samples, + GENOME=genomes, + ) def get_WPS_background_ref(sample): ref_samples = samples["ref_samples"][sample].split(",") @@ -81,6 +90,15 @@ def get_STARTS_background_ref(sample): GENOME=genomes, ) +def get_MP_background_ref(sample): + ref_samples = samples["ref_samples"][sample].split(",") + genomes = samples["genome_build"][ref_samples].values.tolist() + return expand( + "results/intermediate/{{ID}}/table/{GENOME}/background/{{target_region}}--{ref_SAMPLE}_MIDPOINT.background.csv.gz", + zip, + ref_SAMPLE=ref_samples, + GENOME=genomes, + ) def get_length(input): if exists(input): @@ -92,6 +110,18 @@ def get_length(input): #print(f"{input} does not exist: using length of {length}") return length +midpoint_file=list() +if config["midpoint_coverage"] == True: + midpoint_file.append(expand(expand( + "results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_MIDPOINT.csv.gz", + zip, + SAMPLE=samples["sample"], + ID=samples["ID"], + GENOME=samples["genome_build"], + allow_missing=True, + ),target_region=regions["target"], + )) + rule all: input: @@ -152,6 +182,7 @@ rule all: ), target_region=regions["target"], ), + midpoint_file rule add_flanks: @@ -207,76 +238,168 @@ rule generate_random_background: gzip -c > {output} """ - -rule extract_counts: - input: - target="results/intermediate/{ID}/regions/{GENOME}/target_region/{target_region}_blacklist-excluded.bed.gz", - BAMFILE=lambda wildcards: samples["path"][wildcards.SAMPLE], - output: - WPS="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_WPS.csv.gz", - COV="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_COV.csv.gz", - STARTS="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_STARTS.csv.gz", - params: - minRL=config["minRL"], - maxRL=config["maxRL"], - out_pre="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_%s.csv.gz", - conda: - "workflow/envs/cfDNA.yml" - shell: - """ - workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py \ - --minInsert={params.minRL} \ - --maxInsert={params.maxRL} \ - -i {input.target} \ - -o {params.out_pre} {input.BAMFILE} - """ - - -rule extract_counts_background: - input: - background="results/intermediate/{ID}/regions/{GENOME}/background/{target_region}_background_regions.bed.gz", - BAMFILE=lambda wildcards: samples["path"][wildcards.SAMPLE], - output: - WPS="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_WPS.background.csv.gz", - COV="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_COV.background.csv.gz", - STARTS="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_STARTS.background.csv.gz", - params: - minRL=config["minRL"], - maxRL=config["maxRL"], - out_pre="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_%s.background.csv.gz", - conda: - "workflow/envs/cfDNA.yml" - shell: - """ - workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py \ - --minInsert={params.minRL} \ - --maxInsert={params.maxRL} \ - -i {input.background} \ - -o {params.out_pre} {input.BAMFILE} - """ - - -rule plot_overlays: - input: - WPS=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/target/{{target_region}}--{{SAMPLE}}_WPS.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), - WPS_ref=lambda wildcards: get_WPS_ref(wildcards.SAMPLE), - COV=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/target/{{target_region}}--{{SAMPLE}}_COV.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), - COV_ref=lambda wildcards: get_COV_ref(wildcards.SAMPLE), - WPS_back=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/background/{{target_region}}--{{SAMPLE}}_WPS.background.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), - WPS_back_ref=lambda wildcards: get_WPS_background_ref(wildcards.SAMPLE), - COV_back=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/background/{{target_region}}--{{SAMPLE}}_COV.background.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), - COV_back_ref=lambda wildcards: get_COV_background_ref(wildcards.SAMPLE), - output: - "results/plots/overlays/{ID}/{target_region}--{SAMPLE}_overlays.pdf", - params: - target="{target_region}", - sample="{SAMPLE}", - ref_IDs=lambda wildcards: samples["ref_samples"][wildcards.SAMPLE].split(","), - overlay_mode = config["plotting"]["overlay_mode"], - smoothing = config["plotting"]["smoothing"], - rolling = config["plotting"]["rolling"], - background_norm = config["plotting"]["background_norm"] - conda: - "workflow/envs/overlays.yml" - script: - "workflow/scripts/WPS/overlays.py" +if config["midpoint_coverage"] == True: + rule extract_counts: + input: + target="results/intermediate/{ID}/regions/{GENOME}/target_region/{target_region}_blacklist-excluded.bed.gz", + BAMFILE=lambda wildcards: samples["path"][wildcards.SAMPLE], + output: + WPS="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_WPS.csv.gz", + COV="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_COV.csv.gz", + STARTS="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_STARTS.csv.gz", + MIDPOINT="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_MIDPOINT.csv.gz", + params: + minRL=config["minRL"], + maxRL=config["maxRL"], + out_pre="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_%s.csv.gz", + conda: + "workflow/envs/cfDNA.yml" + shell: + """ + workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py \ + --minInsert={params.minRL} \ + --maxInsert={params.maxRL} \ + -i {input.target} \ + --compute_midpoint_coverage=True \ + -o {params.out_pre} {input.BAMFILE} \ + """ + + + rule extract_counts_background: + input: + background="results/intermediate/{ID}/regions/{GENOME}/background/{target_region}_background_regions.bed.gz", + BAMFILE=lambda wildcards: samples["path"][wildcards.SAMPLE], + output: + WPS="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_WPS.background.csv.gz", + COV="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_COV.background.csv.gz", + STARTS="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_STARTS.background.csv.gz", + MIDPOINT="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_MIDPOINT.background.csv.gz", + params: + minRL=config["minRL"], + maxRL=config["maxRL"], + out_pre="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_%s.background.csv.gz", + conda: + "workflow/envs/cfDNA.yml" + shell: + """ + workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py \ + --minInsert={params.minRL} \ + --maxInsert={params.maxRL} \ + -i {input.background} \ + --compute_midpoint_coverage=True \ + -o {params.out_pre} {input.BAMFILE} + """ + + + rule plot_overlays: + input: + WPS=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/target/{{target_region}}--{{SAMPLE}}_WPS.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + WPS_ref=lambda wildcards: get_WPS_ref(wildcards.SAMPLE), + COV=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/target/{{target_region}}--{{SAMPLE}}_COV.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + COV_ref=lambda wildcards: get_COV_ref(wildcards.SAMPLE), + WPS_back=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/background/{{target_region}}--{{SAMPLE}}_WPS.background.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + WPS_back_ref=lambda wildcards: get_WPS_background_ref(wildcards.SAMPLE), + COV_back=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/background/{{target_region}}--{{SAMPLE}}_COV.background.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + COV_back_ref=lambda wildcards: get_COV_background_ref(wildcards.SAMPLE), + MP=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/target/{{target_region}}--{{SAMPLE}}_MIDPOINT.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + MP_ref=lambda wildcards: get_MP_ref(wildcards.SAMPLE), + MP_back=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/background/{{target_region}}--{{SAMPLE}}_MIDPOINT.background.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + MP_back_ref=lambda wildcards: get_MP_background_ref(wildcards.SAMPLE), + output: + "results/plots/overlays/{ID}/{target_region}--{SAMPLE}_overlays.pdf", + params: + target="{target_region}", + sample="{SAMPLE}", + ref_IDs=lambda wildcards: samples["ref_samples"][wildcards.SAMPLE].split(","), + overlay_mode = config["plotting"]["overlay_mode"], + smoothing = config["plotting"]["smoothing"], + rolling = config["plotting"]["rolling"], + background_norm = config["plotting"]["background_norm"], + edge_norm = config["plotting"]["edge_norm"], + win_len=config["plotting"]["win_len"], + poly=config["plotting"]["poly"], + compute_midpoint_coverage=config["midpoint_coverage"], + conda: + "workflow/envs/overlays.yml" + script: + "workflow/scripts/WPS/overlays.py" + +else: + rule extract_counts: + input: + target="results/intermediate/{ID}/regions/{GENOME}/target_region/{target_region}_blacklist-excluded.bed.gz", + BAMFILE=lambda wildcards: samples["path"][wildcards.SAMPLE], + output: + WPS="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_WPS.csv.gz", + COV="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_COV.csv.gz", + STARTS="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_STARTS.csv.gz", + params: + minRL=config["minRL"], + maxRL=config["maxRL"], + out_pre="results/intermediate/{ID}/table/{GENOME}/target/{target_region}--{SAMPLE}_%s.csv.gz", + conda: + "workflow/envs/cfDNA.yml" + shell: + """ + workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py \ + --minInsert={params.minRL} \ + --maxInsert={params.maxRL} \ + -i {input.target} \ + --compute_midpoint_coverage=False \ + -o {params.out_pre} {input.BAMFILE} + """ + + + rule extract_counts_background: + input: + background="results/intermediate/{ID}/regions/{GENOME}/background/{target_region}_background_regions.bed.gz", + BAMFILE=lambda wildcards: samples["path"][wildcards.SAMPLE], + output: + WPS="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_WPS.background.csv.gz", + COV="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_COV.background.csv.gz", + STARTS="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_STARTS.background.csv.gz", + params: + minRL=config["minRL"], + maxRL=config["maxRL"], + out_pre="results/intermediate/{ID}/table/{GENOME}/background/{target_region}--{SAMPLE}_%s.background.csv.gz", + conda: + "workflow/envs/cfDNA.yml" + shell: + """ + workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py \ + --minInsert={params.minRL} \ + --maxInsert={params.maxRL} \ + -i {input.background} \ + --compute_midpoint_coverage=False \ + -o {params.out_pre} {input.BAMFILE} + """ + + + rule plot_overlays: + input: + WPS=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/target/{{target_region}}--{{SAMPLE}}_WPS.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + WPS_ref=lambda wildcards: get_WPS_ref(wildcards.SAMPLE), + COV=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/target/{{target_region}}--{{SAMPLE}}_COV.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + COV_ref=lambda wildcards: get_COV_ref(wildcards.SAMPLE), + WPS_back=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/background/{{target_region}}--{{SAMPLE}}_WPS.background.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + WPS_back_ref=lambda wildcards: get_WPS_background_ref(wildcards.SAMPLE), + COV_back=lambda wc: "results/intermediate/{{ID}}/table/{GENOME}/background/{{target_region}}--{{SAMPLE}}_COV.background.csv.gz".format(GENOME=samples["genome_build"].loc[samples["sample"] == wc.SAMPLE].values[0]), + COV_back_ref=lambda wildcards: get_COV_background_ref(wildcards.SAMPLE), + output: + "results/plots/overlays/{ID}/{target_region}--{SAMPLE}_overlays.pdf", + params: + target="{target_region}", + sample="{SAMPLE}", + ref_IDs=lambda wildcards: samples["ref_samples"][wildcards.SAMPLE].split(","), + overlay_mode = config["plotting"]["overlay_mode"], + smoothing = config["plotting"]["smoothing"], + rolling = config["plotting"]["rolling"], + background_norm = config["plotting"]["background_norm"], + edge_norm = config["plotting"]["edge_norm"], + win_len=config["plotting"]["win_len"], + poly=config["plotting"]["poly"], + compute_midpoint_coverage=config["midpoint_coverage"], + conda: + "workflow/envs/overlays.yml" + script: + "workflow/scripts/WPS/overlays.py" diff --git a/workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py b/workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py index bebdb1a..966393d 100755 --- a/workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py +++ b/workflow/scripts/WPS/extractFromBAM_RegionBed_WPS_Cov.py @@ -79,7 +79,11 @@ def parseRegion(regionstr): parser.add_argument("--downsample", dest="downsample", help="Ratio to down sample reads (default OFF)",default=None,type=float) parser.add_argument("--onefile", dest="onefile", help="Print as single output to stdout (default OFF)",default=False,action="store_true") parser.add_argument("-v","--verbose", dest="verbose", help="Turn debug output on",default=False,action="store_true") +parser.add_argument("-a","--compute_midpoint_coverage", dest="compute_midpoint_coverage", help="compute midpoint coverage score additionally to WPS", default="False") options = parser.parse_args() +size_range_min = 100 +size_range_max = 200 +map_q = 20 minInsSize,maxInsSize = None,None if options.minInsSize > 0 and options.maxInsSize > 0 and options.minInsSize < options.maxInsSize: @@ -90,7 +94,10 @@ def parseRegion(regionstr): options.outfile = options.outfile.strip("""\'""") outfiles = {} if not options.onefile: - outfiles = { 'WPS':gzip.open(options.outfile%"WPS",'wt'), 'COV':gzip.open(options.outfile%"COV",'wt'), 'STARTS':gzip.open(options.outfile%"STARTS",'wt') } + if options.compute_midpoint_coverage == "True": + outfiles = { 'WPS':gzip.open(options.outfile%"WPS",'wt'), 'COV':gzip.open(options.outfile%"COV",'wt'), 'STARTS':gzip.open(options.outfile%"STARTS",'wt'), 'MIDPOINT':gzip.open(options.outfile%"MIDPOINT",'wt') } + else: + outfiles = { 'WPS':gzip.open(options.outfile%"WPS",'wt'), 'COV':gzip.open(options.outfile%"COV",'wt'), 'STARTS':gzip.open(options.outfile%"STARTS",'wt') } if maxInsSize: edge_extension = maxInsSize @@ -126,7 +133,7 @@ def parseRegion(regionstr): if options.verbose: sys.stderr.write("Processing region: %s:%s-%s %s %s %s\n"%(chrom,start,end,cid,score,strand)) - posRange = defaultdict(lambda:[0,0]) + posRange = defaultdict(lambda:[0,0,0]) filteredReads = Intersecter() for bamfile in options.files: @@ -141,6 +148,22 @@ def parseRegion(regionstr): break if options.verbose: sys.stderr.write("Retrieving reads...\n") for read in input_file.fetch(prefix+chrom,max(regionStart-edge_extension-1,0),regionEnd+edge_extension+1): + #Midpoint coverage implementation + if read.is_paired==True and read.is_duplicate==False and read.is_qcfail==False and read.mapq>=map_q: + rstart = min(read.pos,read.pnext)+1 # 1-based + lseq = abs(read.isize) + rend = rstart+lseq-1 # end included + if abs(read.isize)>=size_range_min and abs(read.isize)<=size_range_max: + tag = 1 + # if gc_correct: + # try: + # tag=round(dict(read.tags)["YC"],2) + # except KeyError as e: + # tag=1 + midpoint = int(math.floor((rstart+rend)/2)) + if midpoint in range(rstart,rend+1): + posRange[midpoint][2] += tag + if read.is_duplicate or read.is_qcfail or read.is_unmapped: continue if isSoftClipped(read.cigar): continue @@ -200,6 +223,7 @@ def parseRegion(regionstr): wps_list = [] cov_list = [] starts_list = [] + mp_list = [] for pos in range(regionStart,regionEnd+1): rstart,rend = pos-protection,pos+protection @@ -211,7 +235,7 @@ def parseRegion(regionstr): ecount += 1.0 else: gcount += 1.0 - covCount,startCount = posRange[pos] + covCount,startCount,mpValue = posRange[pos] cov_sites += covCount wpsValue = gcount-(bcount+ecount) #if (options.method != "WPSv1") and (2*gcount+bcount+ecount > 1): @@ -228,11 +252,16 @@ def parseRegion(regionstr): #elif (options.method != "WPSv1"): # wpsValue = 0.0 if options.onefile: - outLines.append("%s\t%d\t%.4f\t%.4f\t%.4f\n"%(chrom,pos,covCount,startCount,wpsValue)) + if options.compute_midpoint_coverage == "True": + outLines.append("%s\t%d\t%.4f\t%.4f\t%.4f\n"%(chrom,pos,covCount,startCount,wpsValue,mpValue)) + else: + outLines.append("%s\t%d\t%.4f\t%.4f\t%.4f\n"%(chrom,pos,covCount,startCount,wpsValue)) else: wps_list.append(wpsValue) cov_list.append(covCount) starts_list.append(startCount) + if options.compute_midpoint_coverage == "True": + mp_list.append(mpValue) if options.onefile: if strand == "-": outLines = outLines[::-1] @@ -246,6 +275,10 @@ def parseRegion(regionstr): outfiles['COV'].write(cid+","+",".join(map(lambda x:str(round(x,5)).replace(".0",""),cov_list))+"\n") outfiles['STARTS'].write(cid+","+",".join(map(lambda x:str(round(x,5)).replace(".0",""),starts_list))+"\n") + if options.compute_midpoint_coverage == "True": + if strand == "-": mp_list = mp_list[::-1] + outfiles['MIDPOINT'].write(cid+","+",".join(map(lambda x:str(round(x,5)),mp_list))+"\n") + if not options.onefile: for name,filestream in outfiles.items(): filestream.close() diff --git a/workflow/scripts/WPS/overlays.py b/workflow/scripts/WPS/overlays.py index f58b8a8..2df278c 100755 --- a/workflow/scripts/WPS/overlays.py +++ b/workflow/scripts/WPS/overlays.py @@ -36,8 +36,19 @@ smoothing = snakemake.params["smoothing"] rolling = snakemake.params["rolling"] background_norm = snakemake.params["background_norm"] +edge_norm = snakemake.params["edge_norm"] flank_edge = 500 +win_len = snakemake.params["win_len"] +poly = snakemake.params["poly"] +compute_midpoint_coverage = snakemake.params["compute_midpoint_coverage"] +if compute_midpoint_coverage == True: + MP = snakemake.input["MP"] + MP_refs = snakemake.input["MP_ref"] + MP_back = snakemake.input["MP_back"] + MP_back_refs = snakemake.input["MP_back_ref"] + + # def functions @@ -120,9 +131,9 @@ def add_sample(path_a: str, path_b: str, overlay_mode:str = "mean",smoothing:boo if smoothing: if rolling: - sample = sample.apply(lambda x:savgol_filter(x,window_length=21, polyorder=2)) - sample.apply(lambda x:savgol_filter(x,window_length=21, polyorder=2)).rolling(1000, center=True).median() + sample = sample.apply(lambda x:savgol_filter(x,window_length=win_len, polyorder=poly)) - sample.apply(lambda x:savgol_filter(x,window_length=21, polyorder=2)).rolling(1000, center=True).median() else: - sample = sample.apply(lambda x:savgol_filter(x,window_length=21, polyorder=2)) + sample = sample.apply(lambda x:savgol_filter(x,window_length=win_len, polyorder=poly)) else: if rolling: sample = sample - sample.rolling(1000, center=True).median() @@ -139,49 +150,131 @@ def add_sample(path_a: str, path_b: str, overlay_mode:str = "mean",smoothing:boo # average over all regions per sample and substract the trimmed mean to normalise -av_WPS = pd.DataFrame(add_sample(WPS, WPS_back,overlay_mode,smoothing,rolling,background_norm)) +data_sample = pd.DataFrame(add_sample(WPS, WPS_back,overlay_mode,smoothing,rolling,background_norm)) sys.stderr.write("WPS [%s]: %s %s\n"%(sample_ID,COV, COV_back)) +if(edge_norm == True): + edge_mean = pd.concat([data_sample.head(100), data_sample.tail(100)]).mean() + av_WPS = pd.DataFrame(data_sample/edge_mean) +else: + av_WPS = pd.DataFrame(data_sample) av_WPS.columns = av_WPS.columns.astype(str) av_WPS.columns.values[-1] = sample_ID for (ref_ID, WPS_ref, WPS_back_ref) in zip(ref_IDs, WPS_refs, WPS_back_refs): sys.stderr.write("WPS [%s]: %s %s\n"%(ref_ID, WPS_ref, WPS_back_ref)) - av_WPS[ref_ID] = add_sample(WPS_ref, WPS_back_ref,overlay_mode,smoothing,rolling,background_norm)["value"] + data_reference = add_sample(WPS_ref, WPS_back_ref,overlay_mode,smoothing,rolling,background_norm)["value"] + if(edge_norm == True): + edge_mean = pd.concat([data_reference.head(100), data_reference.tail(100)]).mean() + av_WPS[ref_ID] = data_reference/edge_mean + else: + av_WPS[ref_ID] = data_reference -av_COV = pd.DataFrame(add_sample(COV, COV_back,overlay_mode,smoothing,rolling, background_norm)) +data_sample = pd.DataFrame(add_sample(COV, COV_back,overlay_mode,smoothing,rolling,background_norm)) sys.stderr.write("COV [%s]: %s %s\n"%(sample_ID, COV, COV_back)) +if(edge_norm == True): + edge_mean = pd.concat([data_sample.head(100), data_sample.tail(100)]).mean() + av_COV = pd.DataFrame(data_sample/edge_mean) +else: + av_COV = pd.DataFrame(data_sample) av_COV.columns = av_COV.columns.astype(str) av_COV.columns.values[-1] = sample_ID for (ref_ID, COV_ref, COV_back_ref) in zip(ref_IDs, COV_refs, COV_back_refs): sys.stderr.write("COV [%s]: %s %s\n"%(ref_ID, COV_ref, COV_back_ref)) - av_COV[ref_ID] = add_sample(COV_ref, COV_back_ref,overlay_mode,smoothing,rolling,background_norm)["value"] + data_reference = add_sample(COV_ref, COV_back_ref,overlay_mode,smoothing,rolling,background_norm)["value"] + if(edge_norm == True): + edge_mean = pd.concat([data_reference.head(100), data_reference.tail(100)]).mean() + av_COV[ref_ID] = data_reference/edge_mean + else: + av_COV[ref_ID] = data_reference + +if compute_midpoint_coverage == True: + data_sample = pd.DataFrame(add_sample(MP, MP_back,overlay_mode,smoothing,rolling,background_norm)) + sys.stderr.write("MP [%s]: %s %s\n"%(sample_ID, MP, MP_back)) + if(edge_norm == True): + edge_mean = pd.concat([data_sample.head(100), data_sample.tail(100)]).mean() + av_MP = pd.DataFrame(data_sample/edge_mean) + else: + av_MP = pd.DataFrame(data_sample) + av_MP.columns = av_MP.columns.astype(str) + av_MP.columns.values[-1] = sample_ID + for (ref_ID, MP_ref, MP_back_ref) in zip(ref_IDs, MP_refs, MP_back_refs): + sys.stderr.write("MP [%s]: %s %s\n"%(ref_ID, MP_ref, MP_back_ref)) + data_reference = add_sample(MP_ref, MP_back_ref,overlay_mode,smoothing,rolling,background_norm)["value"] + if(edge_norm == True): + edge_mean = pd.concat([data_reference.head(100), data_reference.tail(100)]).mean() + av_MP[ref_ID] = data_reference/edge_mean + else: + av_MP[ref_ID] = data_reference # create line plots and save to a single pdf -if overlay_mode.lower() == "confidence": - av_WPS_long=av_WPS.reset_index().melt(id_vars=["position", "sample_nr"], value_name="score", var_name="phenotype").sort_values(by="position") - av_COV_long=av_COV.reset_index().melt(id_vars=["position", "sample_nr"], value_name="score", var_name="phenotype").sort_values(by="position") - Fig_WPS = sns.lineplot(data=av_WPS_long, x="position", y="score", hue="phenotype",) - plt.suptitle("adjusted WPS: {target} target regions") - plt.xlabel("Position relative to target site") - plt.ylabel("normalized WPS") - plt.close() - Fig_Cov = sns.lineplot(data=av_COV_long, x="position", y="score", hue="phenotype",) - plt.suptitle(f"adjusted read coverage: {target} target regions") - plt.xlabel("Position relative to target site") - plt.ylabel("normalized read coverage") - plt.close() +if compute_midpoint_coverage == True: + if overlay_mode.lower() == "confidence": + av_WPS_long=av_WPS.reset_index().melt(id_vars=["position", "sample_nr"], value_name="score", var_name="phenotype").sort_values(by="position") + av_COV_long=av_COV.reset_index().melt(id_vars=["position", "sample_nr"], value_name="score", var_name="phenotype").sort_values(by="position") + av_MP_long=av_MP.reset_index().melt(id_vars=["position", "sample_nr"], value_name="score", var_name="phenotype").sort_values(by="position") + Fig_WPS = sns.lineplot(data=av_WPS_long, x="position", y="score", hue="phenotype",) + plt.suptitle("adjusted WPS: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized WPS") + plt.close() + Fig_Cov = sns.lineplot(data=av_COV_long, x="position", y="score", hue="phenotype",) + plt.suptitle(f"adjusted read coverage: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized read coverage") + plt.close() + Fig_MP = sns.lineplot(data=av_MP_long, x="position", y="score", hue="phenotype",) + plt.suptitle(f"adjusted midpoint coverage: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized midpoint coverage") + plt.close() + else: + Fig_WPS = sns.lineplot(data=av_WPS) + plt.suptitle(f"adjusted WPS: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized WPS") + plt.close() + Fig_Cov = sns.lineplot(data=av_COV) + plt.suptitle(f"adjusted read coverage: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized read coverage") + plt.close() + Fig_MP = sns.lineplot(data=av_MP) + plt.suptitle(f"adjusted midpoint coverage: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized midpoint coverage") + plt.close() + + with PdfPages(outfile) as pdf: + pdf.savefig(Fig_WPS.get_figure()) + pdf.savefig(Fig_Cov.get_figure()) + pdf.savefig(Fig_MP.get_figure()) + else: - Fig_WPS = sns.lineplot(data=av_WPS) - plt.suptitle(f"adjusted WPS: {target} target regions") - plt.xlabel("Position relative to target site") - plt.ylabel("normalized WPS") - plt.close() - Fig_Cov = sns.lineplot(data=av_COV) - plt.suptitle(f"adjusted read coverage: {target} target regions") - plt.xlabel("Position relative to target site") - plt.ylabel("normalized read coverage") - plt.close() - -with PdfPages(outfile) as pdf: - pdf.savefig(Fig_WPS.get_figure()) - pdf.savefig(Fig_Cov.get_figure()) + if overlay_mode.lower() == "confidence": + av_WPS_long=av_WPS.reset_index().melt(id_vars=["position", "sample_nr"], value_name="score", var_name="phenotype").sort_values(by="position") + av_COV_long=av_COV.reset_index().melt(id_vars=["position", "sample_nr"], value_name="score", var_name="phenotype").sort_values(by="position") + Fig_WPS = sns.lineplot(data=av_WPS_long, x="position", y="score", hue="phenotype",) + plt.suptitle("adjusted WPS: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized WPS") + plt.close() + Fig_Cov = sns.lineplot(data=av_COV_long, x="position", y="score", hue="phenotype",) + plt.suptitle(f"adjusted read coverage: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized read coverage") + plt.close() + else: + Fig_WPS = sns.lineplot(data=av_WPS) + plt.suptitle(f"adjusted WPS: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized WPS") + plt.close() + Fig_Cov = sns.lineplot(data=av_COV) + plt.suptitle(f"adjusted read coverage: {target} target regions") + plt.xlabel("Position relative to target site") + plt.ylabel("normalized read coverage") + plt.close() + + with PdfPages(outfile) as pdf: + pdf.savefig(Fig_WPS.get_figure()) + pdf.savefig(Fig_Cov.get_figure())