add annotation for promoter regions #41

Author: sreichl

workflow/rules/quantification.smk

@@ -25,6 +25,7 @@ rule get_promoter_regions:
         config["gencode_gtf"],
     output:
         promoter_regions = os.path.join(result_path,"counts","promoter_regions.bed"),
+        promoter_annot = os.path.join(result_path,"counts","promoter_annotation.csv"),
     params:
         # cluster parameters
         partition=config.get("partition"),

workflow/scripts/get_promoter_regions.py

@@ -2,6 +2,7 @@
 
 #### libraries
 import pybedtools as bedtools
+import pandas as pd
 
 # extract promoter regions
 def get_promoter(feature, upstream, downstream, chrom_sizes):
@@ -26,7 +27,7 @@ def get_promoter(feature, upstream, downstream, chrom_sizes):
         start,
         end,
         gene_id,
-#         feature.attrs['gene_name'] if 'gene_name' in feature.attrs else feature.attrs['gene_id'],
+        feature.attrs['gene_name'] if 'gene_name' in feature.attrs else feature.attrs['gene_id'],
 #         '.',
 #         feature.strand
     ])
@@ -41,10 +42,12 @@ def get_promoter(feature, upstream, downstream, chrom_sizes):
 
 # output
 promoter_regions_path = snakemake.output["promoter_regions"]
+promoter_annot_path = snakemake.output["promoter_annot"]
 
 # parameters
 TSS_up = snakemake.config["proximal_size_up"]
 TSS_dn = snakemake.config["proximal_size_dn"]
+genome_fasta_path = snakemake.config["genome_fasta"]
 
 # load the genome annotation file using pybedtools
 gtf = bedtools.BedTool(gtf_file)
@@ -68,3 +71,14 @@ def get_promoter(feature, upstream, downstream, chrom_sizes):
 
 # save the promoter regions to a BED file
 promoters.saveas(promoter_regions_path)
+
+# calculate GC content and length for each region and save as annotation
+gc_content_length = promoters.nucleotide_content(fi=genome_fasta_path).to_dataframe()
+gc_content_length.columns = [col.split('_', 1)[-1].replace('at', 'AT').replace('gc', 'GC').replace('oth', 'otherBases') for col in gc_content_length.columns]
+gc_content_length = gc_content_length.add_prefix('bedtools_')
+gc_content_length.columns = ["chr", "start", "end", "gene", "gene_name"] + gc_content_length.columns[5:].tolist()
+gc_content_length.set_index("gene", inplace=True)
+gc_content_length.to_csv(promoter_annot_path)
+
+# load, remove last column (gene name) and save again as final promoter BED file for quantification
+bedtools.BedTool(promoter_regions_path).cut(range(0, 4)).saveas(promoter_regions_path)

workflow/scripts/map_consensus_tss.py

@@ -44,13 +44,10 @@ def map_region(x):
 annot_regions.set_index('peak_id', inplace=True)
 TSS_annot = annot_regions.loc[TSS_regions["peak_id"],:]
 TSS_annot.reset_index(inplace=True)
-TSS_annot = TSS_annot.sort_values(by="peak_id")
 TSS_annot.index = TSS_regions.index
-
 TSS_annot.to_csv(tss_annot_path)
 
 # save bed file of TSS regions
-TSS_annot.reset_index(inplace=True)
-TSS_bed_df = TSS_annot[["gencode_chr",  'gencode_start', 'gencode_end', 'homer_Nearest_Ensembl']]
-TSS_bed = pybedtools.BedTool.from_dataframe(TSS_bed_df)
+TSS_bed_df = TSS_annot.sort_values(by="peak_id")[["gencode_chr",  'gencode_start', 'gencode_end', 'homer_Nearest_Ensembl']]
+TSS_bed = bedtools.BedTool.from_dataframe(TSS_bed_df)
 TSS_bed.saveas(tss_bed_path)