Reformatting the output

Author: weichan

config/config_CD19_BBz.yml

@@ -36,6 +36,8 @@ annotate_ucsc_exons: "/home/weichan/permanent/Projects/VIS/dev/UCSC/UCSC_GENCODE
 annotate_ucsc_introns: "/home/weichan/permanent/Projects/VIS/dev/UCSC/UCSC_GENCODE_V44_Introns"
 #promoters
 annotate_ucsc_promoter: "/home/weichan/permanent/Projects/VIS/dev/UCSC/UCSC_Promoters_EPDnew"
+#enhancer
+annotate_ucsc_enhancer: "/home/weichan/permanent/Projects/VIS/dev/UCSC/UCSC_VISTA_Enhancer_02_2025.bed"
 #sedb
 #sedb_cd4: "/home/weichan/permanent/Projects/VIS/dev/SEdb/sedb_cd4_closest_gene.bed"
 #sedb_cd8: "/home/weichan/permanent/Projects/VIS/dev/SEdb/sedb_cd8_closest_gene.bed"
@@ -61,6 +63,6 @@ downstream: "rules/downstream.smk"
 #epigenetics: "rules/epigenetics.smk"
 #variants rule collection WIP
 #variants: "rules/variants.smk"
-cmod: "rules/cmod.smk"
+base_modifications: "rules/base_modifications.smk"
 plot_functional_genomics: "rules/plot_functional_genomics.smk"
 threads: 10

workflow/Snakefile

@@ -20,23 +20,15 @@ cwd=os.getcwd()
 
 import VIS_helper_functions as vhf #custom functions to make snakemake pipeline leaner
 
-#inmport rules
+#inmport default rules
 include: config["detection"]
 include: config["quality_control"]
 
-#analysis-specific rules
-#include: config["downstream"]
-
-#development
-#include: config["epigenetics"]
-#include: config["variants"]
-#include: config["development"]
-
-
 #conditional rule all based on defined rules
 conditional_output = list()
 
 if "functional_genomics" in config:
+	import VIS_functional_genomics_helper_functions as vhf_fg
 	include: config["functional_genomics"]
 	conditional_output.append(expand(f"{outdir}/final/functional_genomics/Functional_distances_to_Insertions_{{sample}}.bed", sample=SAMPLES))
 	conditional_output.append(expand(f"{outdir}/final/functional_genomics/Functional_distances_to_Insertions_{{sample}}.bed", sample=SAMPLES))
@@ -45,17 +37,15 @@ if "functional_genomics" in config:
 									sample=SAMPLES))
 
 if "plot_functional_genomics" in config:
-
-	import VIS_plot_functional_annotation_helper_functions as vhf_fa
+	import VIS_plot_functional_genomics_helper_functions as vhf_pfg
 	include: config["plot_functional_genomics"]
-
 	conditional_output.append(expand(f"{outdir}/final/functional_genomics/Plot_Distance_to_Genes_{fragmentsize}_{{sample}}.png", sample=SAMPLES))
 	conditional_output.append(expand(f"{outdir}/final/functional_genomics/Insertion_Scoring_{{sample}}.svg", sample=SAMPLES))
 
-if "cmod" in config:
-	include: config["cmod"]
-	conditional_output.append(expand(f"{outdir}/intermediate/cmod/Isolated_Reads_{{sample}}.tsv", sample=SAMPLES))
-	conditional_output.append(expand(f"{outdir}/intermediate/cmod/Calls_Isolated_Reads_{{sample}}.tsv", sample=SAMPLES))
+if "base_modifications" in config:
+	include: config["base_modifications"]
+	conditional_output.append(expand(f"{outdir}/final/base_modifications/Isolated_Reads_{{sample}}.tsv", sample=SAMPLES))
+	conditional_output.append(expand(f"{outdir}/final/base_modifications/Calls_Isolated_Reads_{{sample}}.tsv", sample=SAMPLES))
 
 rule all:
 	input:

workflow/rules/cmod.smk workflow/rules/base_modifications.smkworkflow/rules/cmod.smk renamed to workflow/rules/base_modifications.smk

@@ -1,6 +1,6 @@
 ######
 ######
-###### C-modifications of reads with isnertions
+###### Base modifications of reads with isnertions
 ######
 ######
 
@@ -11,9 +11,9 @@ rule readnames_final:
 	input:
 		final_insertions=f"{outdir}/final/localization/ExactInsertions_{{sample}}.bed"
 	log:
-		log=f"{outdir}/intermediate/log/cmod/readnames_final/{{sample}}.log"
+		log=f"{outdir}/intermediate/log/base_modifications/readnames_final/{{sample}}.log"
 	output:
-		readnames=f"{outdir}/intermediate/cmod/final_insertion_readnames_{{sample}}.txt"
+		readnames=temp(f"{outdir}/intermediate/base_modifications/final_insertion_readnames_{{sample}}.txt")
 	conda:
 		"../envs/VIS_dummy_env.yml"
 	shell:
@@ -23,15 +23,14 @@ rule readnames_final:
 		) > {log.log} 2>&1
 		"""
 
-
-rule only_keep_valid_insertions:
+rule only_keep_filtered_insertions:
 	input:
 		isobam=f"{outdir}/intermediate/mapping/Precut_{{sample}}_sorted.bam",
 		readnames=f"{outdir}/intermediate/cmod/final_insertion_readnames_{{sample}}.txt"
 	log:
-		log=f"{outdir}/intermediate/log/cmod/filter/{{sample}}.log"
+		log=f"{outdir}/intermediate/log/base_modifications/only_keep_filtered_insertions/{{sample}}.log"
 	output:
-		bam=f"{outdir}/intermediate/cmod/Final_Isolated_Reads_{{sample}}.bam"
+		bam=f"{outdir}/intermediate/base_modifications/Final_Isolated_Reads_{{sample}}.bam"
 	conda:
 		"../envs/VIS_samtools_env.yml"
 	shell:
@@ -42,40 +41,40 @@ rule only_keep_valid_insertions:
 		) > {log.log} 2>&1
 		"""
 
-
 rule modkit:
 	input:
-		isobam=f"{outdir}/intermediate/cmod/Final_Isolated_Reads_{{sample}}.bam"
+		isobam=f"{outdir}/intermediate/base_modifications/Final_Isolated_Reads_{{sample}}.bam"
 	log:
-		log=f"{outdir}/intermediate/log/cmod/modkit/{{sample}}.log"
+		log=f"{outdir}/intermediate/log/base_modifications/modkit/{{sample}}.log"
 	output:
-		tsv=f"{outdir}/intermediate/cmod/Isolated_Reads_{{sample}}.tsv"
+		tsv=f"{outdir}/final/base_modifications/Isolated_Reads_{{sample}}.tsv"
 	conda:
 		"../envs/VIS_modkit_env.yml"
+	threads: config["threads"]
 	shell:
 		"""
 		(
-		modkit extract full -t 20 {input.isobam} {output.tsv}
+		modkit extract full -t {threads} {input.isobam} {output.tsv}
 		) > {log.log} 2>&1
 		"""
 
 rule call_modkit:
 	input:
-		isobam=f"{outdir}/intermediate/cmod/Final_Isolated_Reads_{{sample}}.bam"
+		isobam=f"{outdir}/intermediate/base_modifications/Final_Isolated_Reads_{{sample}}.bam"
 	log:
-		log=f"{outdir}/intermediate/log/cmod/Calls_modkit/{{sample}}.log"
+		log=f"{outdir}/intermediate/log/base_modifications/call_modkit/{{sample}}.log"
 	output:
-		tsv=f"{outdir}/intermediate/cmod/Calls_Isolated_Reads_{{sample}}.tsv"
+		tsv=f"{outdir}/final/base_modifications/Calls_Isolated_Reads_{{sample}}.tsv"
 	conda:
 		"../envs/VIS_modkit_env.yml"
+	threads: config["threads"]
 	shell:
 		"""
 		(
-		modkit extract calls -t 20 {input.isobam} {output.tsv}
+		modkit extract calls -t {threads} {input.isobam} {output.tsv}
 		) > {log.log} 2>&1
 		"""
 
-
 '''
 rule specific_methylartist:
 	input:
@@ -93,7 +92,7 @@ rule specific_methylartist:
 		methylartist locus -b {input.bam} -i chr17:31124037-31180287 -n C -r {input.ref} -l 31154037-31159287 -o {output.plot}
 		) > {log.log} 2>&1
 		"""
-'''
+
 
 rule inserted_seq:
 	input:
@@ -109,3 +108,4 @@ rule inserted_seq:
 		except Exception as e:
 			with open(log.log, "a") as log_file:
 					log_file.write(f"Error: {str(e)}\n")
+'''

workflow/rules/detection.smk

@@ -114,7 +114,6 @@ rule fasta_insertion_reads:
 		) > {log.log} 2>&1
 		"""
 
-
 ######
 ######
 ###### "Clean" BAM: Cut-out fasta to BAM via Mapping to reference 
@@ -262,6 +261,7 @@ rule split_fasta:
 ###### Insertion preparation: Fragmentation 
 ######
 ######
+
 rule prepare_insertion:
 	input:
 		config["insertion_fasta"] #insertion fasta sequence
@@ -484,7 +484,23 @@ rule calculate_exact_insertion_coordinates:
 	    except Exception as e:
 	        with open(log.log, "a") as log_file:
                     log_file.write(f"Error: {str(e)}\n")
-                
+
+rule insertion_points:
+    input:
+        f"{outdir}/intermediate/localization/ExactInsertions_{{sample}}.bed"
+    output:
+        f"{outdir}/final/localization/InsertionPoints_{{sample}}.bed"
+    log:
+        f"{outdir}/intermediate/log/detection/insertion_points/{{sample}}.log"
+    conda:
+        "../envs/VIS_dummy_env.yml"
+    shell:
+        """
+        (
+        awk '{{OFS="\t"; $3 = $2 + 1; print $0}}' {input} > {output}
+        ) > {log} 2>&1
+        """
+               
 rule collect_outputs:
 	input:
 		coordinates=f"{outdir}/intermediate/localization/ExactInsertions_{{sample}}.bed",

workflow/rules/functional_genomics.smk

@@ -4,13 +4,13 @@
 rule sort_insertion_file:
     input:
         exact=f"{outdir}/intermediate/localization/ExactInsertions_{{sample}}.bed",
-        point=f"{outdir}/intermediate/localization/annotation/temp_Insertions_{{sample}}.bed"
+        point=f"{outdir}/final/localization/InsertionPoints_{{sample}}.bed"
     log:
         log=f"{outdir}/intermediate/log/functional_genomics/sort_insertion_file/{{sample}}.log",
-        log2=f"{outdir}/intermediate/log/functional_genomics/sort_insertion_file/temp_{{sample}}.log"
+        log2=f"{outdir}/intermediate/log/functional_genomics/sort_insertion_file/points_{{sample}}.log"
     output:
         exact=temp(f"{outdir}/intermediate/localization/Sorted_ExactInsertions_{{sample}}.bed"),
-        point=f"{outdir}/intermediate/localization/temp_Sorted_ExactInsertions_{{sample}}.bed"
+        point=f"{outdir}/intermediate/localization/Sorted_InsertionPoints_{{sample}}.bed"
     conda:
         "../envs/VIS_dummy_env.yml"
     shell:
@@ -26,7 +26,7 @@ rule sort_insertion_file:
 rule calc_distance_to_elements:
     input:
         #insertions=f"{outdir}/intermediate/localization/Sorted_ExactInsertions_{{sample}}.bed",
-        insertions=f"{outdir}/intermediate/localization/temp_Sorted_ExactInsertions_{{sample}}.bed",
+        insertions=f"{outdir}/intermediate/localization/Sorted_InsertionPoints_{{sample}}.bed",
         ref=config["ref_genome_ctrl"]
     params:
         annotation_files={k: v for k, v in config.items() if k.startswith("annotate_")}
@@ -37,26 +37,10 @@ rule calc_distance_to_elements:
     run:
         vhf.calculate_element_distance(input.insertions, output[0], log.log, params.annotation_files)
 
-rule modify_insertions:
-    input:
-        f"{outdir}/intermediate/localization/ExactInsertions_{{sample}}.bed"
-    output:
-        temp(f"{outdir}/intermediate/localization/annotation/temp_Insertions_{{sample}}.bed")
-    log:
-        f"{outdir}/intermediate/log/functional_genomics/modify_insertions/{{sample}}.log"
-    conda:
-        "../envs/VIS_dummy_env.yml"
-    shell:
-        """
-        (
-        awk '{{OFS="\t"; $3 = $2 + 1; print $0}}' {input} > {output}
-        ) > {log} 2>&1
-        """
-
 #kinda 'hacky' solution for flexibility: Now it does not matter if one annotation or 20 are defined. 
 rule annotation_overlap_insertion:
     input:
-        insertions_bed=f"{outdir}/intermediate/localization/temp_Sorted_ExactInsertions_{{sample}}.bed"
+        insertions_bed=f"{outdir}/intermediate/localization/Sorted_InsertionPoints_{{sample}}.bed"
     params:
         annotation_files={k.replace("annotate_", ""): v for k, v in config.items() if k.startswith("annotate_")}
     log:
@@ -65,5 +49,5 @@ rule annotation_overlap_insertion:
         **{k.replace("annotate_", ""): f"{outdir}/intermediate/localization/annotation/Annotation_{k.replace('annotate_', '')}_Insertions_{{sample}}.bed"
            for k in config if k.startswith("annotate_")}
     run:
-        vhf.run_bedtools_intersect(input.insertions_bed, output, log.log, params.annotation_files)
+        vhf_fg.run_bedtools_intersect(input.insertions_bed, output, log.log, params.annotation_files)
 

workflow/rules/plot_functional_genomics.smk

@@ -11,7 +11,7 @@ rule plot_distance_to_elements:
 		scatter=report(f"{outdir}/final/functional_genomics/Plot_Distance_to_Genes_{fragmentsize}_{{sample}}.png"),
 	run:
 	    try:
-	        vhf_fa.plot_element_distance(input.distancetable, params.distances, params.threshold, output.scatter, log.log1)
+	        vhf_pfg.plot_element_distance(input.distancetable, params.distances, params.threshold, output.scatter, log.log1)
 	    except Exception as e:
 	        with open(log.log1, "a") as log_file:
                     log_file.write(f"Error: {str(e)}\n")
@@ -26,7 +26,7 @@ rule plot_scoring:
         data=(f"{outdir}/final/functional_genomics/Insertion_Scoring_Data_{{sample}}.txt")
     run:
         try:
-            vhf_fa.scoring_insertions(input[0], output.plot, output.data, log.log)
+            vhf_pfg.scoring_insertions(input[0], output.plot, output.data, log.log)
         except Exception as e:
             with open(log.log, "a") as log_file:
                 log_file.write(f"Error: {str(e)}\n")

workflow/scripts/VIS_functional_genomics_helper_functions.py

@@ -0,0 +1,68 @@
+#!/usr/bin/env python3
+
+import pandas as pd
+import pybedtools
+
+#wrapper
+from VIS_helper_functions import redirect_logging
+
+#functional
+@redirect_logging(logfile_param="logfile")
+def calculate_element_distance(insertions_bed, output_bed, logfile, annotation_files):
+    """
+    Calculates distances between insertion sites and genomic annotations using bedtools closest.
+    
+    Parameters:
+    - insertions_bed (str): Path to the BED file containing the insertion sites.
+    - output_bed (str): Path to save the output BED file.
+    - annotation_files (dict): Dict with config_entry as keys and pathways as values.
+    """
+    
+    # At least one annotation input necessary
+    if not annotation_files:
+        raise ValueError("At least one annotation file must be provided.")
+    
+    # Create DataFrame for combined annotations
+    combined_df = pd.DataFrame()
+
+    for tag, file in annotation_files.items():
+        try:
+            df = pd.read_csv(file, sep="\t", header=None, usecols=[0, 1, 2, 3, 4, 5])
+            if df.iloc[0, 0].startswith("chr"): 
+                df["source"] = tag  # Add source column
+                print(f"Loaded {tag}: {df.head()}")
+                combined_df = pd.concat([combined_df, df], ignore_index=True)
+        except:
+            print(f"Error reading {file})")
+            print("BED files are expected to follow BED6 format convention. If more than 6 columns are provided, the first 6 will be used.")
+            continue
+    
+    # Convert combined DataFrame bed object
+    combined_bed = pybedtools.BedTool.from_dataframe(combined_df)
+    sorted_annotations = combined_bed.sort()
+
+    insertions = pybedtools.BedTool(insertions_bed)
+    
+    #bedtools closest operation
+    closest = insertions.closest(sorted_annotations, D="a", filenames=True)
+
+    print(type(closest))
+    print(closest)
+    closest.saveas(output_bed)
+    print(f"Distances calculated and saved to {output_bed}")
+
+@redirect_logging(logfile_param="logfile")
+def run_bedtools_intersect(insertions_bed, output_files, logfile, annotations):
+    """
+    Runs bedtools intersect for each annotation using pybedtools.
+    """
+    insertions = pybedtools.BedTool(insertions_bed)
+    
+    for key, annotation_file in annotations.items():
+        output_file = output_files[key]
+        annotation = pybedtools.BedTool(annotation_file)
+        
+        intersected = insertions.intersect(annotation, wb=True)
+        
+        intersected.saveas(output_file)
+        print(f"Completed: {key}")