<jupyter_start><jupyter_text>Housekeeping genes have significantly fewer TATA boxes than variable genes    now I need to rerun analyses using gat enrichment If binding sites you're mapping are small, need to get the mapability genome containing all regions that are uniquely mappable with reads of 24 bases. https://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg38&g=mappabilitySee https://gat.readthedocs.io/en/latest/tutorialGenomicAnnotation.html Downloaded TATA_boxes.bed and TATA_boxes.fps (both the same, different formats) from EPDUsed the following search parameters for download: FindM Genome Assembly : <NAME> (Feb 2011 TAIR10/araTha1)Series : EPDnew, the Arabidopsis Curated Promoter DatabaseSample : TSS from EPDnew rel 004Repeat masking: off5' border: -50     3' border: 0Search mode: forwardSelection mode : all matches In the end I didn't need these files, can use existing tatabox files for the specific genes of interest (responsive_housekeeping_TATA_box_positive.bed)Copied the chromsizes.chr to data/EPD_promoter_analysis/TATA and converted it into a BED file for the workspace.<jupyter_code>#create a bed file containing all 100 constitutive/responsive promoters with the fourth column annotating whether it's constitutive or responsive

promoters_no_random = promoters.copy()
#drop randCont rows
promoters_no_random = promoters_filtered[~(promoters.gene_type == 'randCont')]
promoters_no_random
promoterbedfile = '../../data/FIMO/responsivepromoters.bed'
promoters = pd.read_table(promoterbedfile, sep='\t', header=None)
cols = ['chr', 'start', 'stop', 'promoter_AGI', 'score', 'strand', 'source', 'feature_name', 'dot2', 'attributes']
promoters.columns = cols
merged = pd.merge(promoters,promoters_no_random, on='promoter_AGI')
merged
merged_reordered = merged[['chr','start','stop','gene_type', 'strand', 'source', 'attributes','promoter_AGI']]
sorted_motifs = merged_reordered.sort_values(['chr','start'])
bed = BedTool.from_dataframe(sorted_motifs).saveas('../../data/EPD_promoter_analysis/TATA/promoters_norandom.bed')
def add_chr_linestart(input_location,output_location):
    """this function removes characters from the start of each line in the input file and sends modified lines to output"""
    output = open(output_location, 'w') #make output file with write capability
    #open input file
    with open(input_location, 'r') as infile:  
        #iterate over lines in file
        for line in infile:
            line = line.strip() # removes hidden characters/spaces
            if line[0].isdigit():
                                 
                line = 'chr' + line   #prepend chr to the beginning of line if starts with a digit
            output.write(line + '\n') #output to new file
    output.close()
add_chr_linestart('../../data/EPD_promoter_analysis/TATA/promoters_norandom.bed', '../../data/EPD_promoter_analysis/TATA/promoters_norandom_renamed.bed')

# #In bash I ran this:
# gat-run.py --ignore-segment-tracks --segments=../../data/EPD_promoter_analysis/responsive_housekeeping_TATA_box_positive.bed `#TATA box annotations` \
#     --annotations=../../data/EPD_promoter_analysis/TATA/promoters_norandom.bed `#100 constitutive/responsive promoter annotations` \
#     --workspace=../../data/EPD_promoter_analysis/TATA/chromsizes.bed `#Arabidopsis chromosome bed file` \
#     --num-samples=1000 --log=../../data/EPD_promoter_analysis/TATA/gat.log > ../../data/EPD_promoter_analysis/TATA/gat_TATA.out
# # note,  --num-threads=7 is currently broken`
# #test run
# gat-run.py --ignore-segment-tracks --segments=../../data/EPD_promoter_analysis/responsive_housekeeping_TATA_box_positive.bed `#TATA box annotations` \
#     --annotations=../../data/EPD_promoter_analysis/TATA/promoters_norandom_renamed.bed `#100 constitutive/responsive promoter annotations` \
#     --workspace=../../data/EPD_promoter_analysis/TATA/chromsizes.bed `#Arabidopsis chromosome bed file` \
#     --num-samples=1000 --log=../../data/EPD_promoter_analysis/TATA/gat.log > ../../data/EPD_promoter_analysis/TATA/gat_TATA.out<jupyter_output><empty_output><jupyter_text>Calculate distance of TATA box from TSS<jupyter_code>cols = ['chrTATA', 'startTATA', 'stopTATA', 'gene_IDTATA','number','strandTATA','TATA_present','promoter_AGI']
TATA.columns = cols
TATA
#merge TATA bed with promoters 
sorted_motifs

TATA_distance = pd.merge(TATA,sorted_motifs, how='inner', on='promoter_AGI')
TATA_distance
#calculate distance between TATA and TSS
TATA_distance.loc[TATA_distance.strand =='+', 'TATAdistance(bp)'] = TATA_distance.startTATA - TATA_distance.stop
TATA_distance.loc[TATA_distance.strand =='-', 'TATAdistance(bp)'] = TATA_distance.start - TATA_distance.startTATA
TATA_distance<jupyter_output><empty_output><jupyter_text>Create distribution plotNote:The y axis is a density, not a probability. The normalized histogram does not show a probability mass function, where the sum the bar heights equals 1; the normalization ensures that the sum of the bar heights times the bar widths equals 1. This is what ensures that the normalized histogram is comparable to the kernel density estimate, which is normalized so that the area under the curve is equal to 1.<jupyter_code>dist_plot = TATA_distance['TATAdistance(bp)']
#create figure with no transparency
dist_plot_fig = sns.distplot(dist_plot).get_figure()
dist_plot_fig.savefig('../../data/plots/TATAbox/TATA_distance_from_extracted_promoters.pdf', format='pdf')

TATA_distance['TATAdistance(bp)']
#Make TATA box segment the actual size - I will set all to 15 bp
TATA_15bp = TATA.copy()
TATA_15bp
#Make TATA box segment the actual size - I will set all to 15 bp
TATA_15bp.loc[TATA_15bp.strandTATA =='+', 'stopTATA'] = TATA_15bp.stopTATA + 14
TATA_15bp.loc[TATA_15bp.strandTATA =='-', 'startTATA'] = TATA_15bp.startTATA - 14
TATA_15bp
#make into bed file
sorted_TATA = TATA_15bp.sort_values(['chrTATA','startTATA'])
bed = BedTool.from_dataframe(sorted_TATA).saveas('../../data/EPD_promoter_analysis/TATA/TATA_15bp.bed')
#extend promoter 3' end by 661 bp (to furthest registered TATA box)
responsive_constitutive_promoters_extended = sorted_motifs.copy()
responsive_constitutive_promoters_extended.loc[responsive_constitutive_promoters_extended.strand =='+', 'stop'] = responsive_constitutive_promoters_extended.stop + 675
responsive_constitutive_promoters_extended.loc[responsive_constitutive_promoters_extended.strand =='-', 'start'] = responsive_constitutive_promoters_extended.start - 675

sorted_proms = responsive_constitutive_promoters_extended.sort_values(['chr','start'])
bed = BedTool.from_dataframe(sorted_proms).saveas('../../data/EPD_promoter_analysis/TATA/responsive_constitutive_promoters_extended.bed')
#add chr to chromosome name
add_chr_linestart('../../data/EPD_promoter_analysis/TATA/responsive_constitutive_promoters_extended.bed', '../../data/EPD_promoter_analysis/TATA/responsive_constitutive_promoters_extended_renamed.bed')
#rerun analysis using nonbidirectional promoters
nonbidirectional_proms_file = '../../data/FIMO/nonbidirectional_proms.bed'
nonbidirectional_proms = pd.read_table(nonbidirectional_proms_file, sep='\t', header=None)
cols3 = ['chr', 'start', 'stop','promoter_AGI','dot1', 'strand','source_bi', 'type','dot2', 'attributes']
nonbidirectional_proms.columns = cols3

nonbidir_const_var_proms = pd.merge(sorted_motifs, nonbidirectional_proms[['promoter_AGI','source_bi']], how='left', on='promoter_AGI')
nonbidir_const_var_proms = nonbidir_const_var_proms[~nonbidir_const_var_proms['source_bi'].isnull()]
nonbidir_const_var_proms
#number of nonbidirectional housekeeping genes
len(nonbidir_const_var_proms[nonbidir_const_var_proms.gene_type == 'housekeeping'])
#number of nonbidirectional variable genes
len(nonbidir_const_var_proms[nonbidir_const_var_proms.gene_type == 'highVar'])
# gat-run.py --ignore-segment-tracks --segments=../../data/EPD_promoter_analysis/TATA/TATA_15bp.bed `#TATA box annotations` \
#     --annotations=../../data/EPD_promoter_analysis/TATA/responsive_constitutive_promoters_extended_renamed.bed `#100 constitutive/responsive promoter annotations` \
#     --workspace=../../data/EPD_promoter_analysis/TATA/chromsizes.bed `#Arabidopsis chromosome bed file` \
#     --num-samples=1000 --log=../../data/EPD_promoter_analysis/TATA/gat.log > ../../data/EPD_promoter_analysis/TATA/gat_TATA.out
#Create file with only the variable promoters
extended_promoters_file = '../../data/EPD_promoter_analysis/TATA/responsive_constitutive_promoters_extended_renamed.bed'
extended_promoters = pd.read_table(extended_promoters_file, sep='\t', header=None)
#make a new gat workspace file with all promoters (first 3 columns)
bed = BedTool.from_dataframe(extended_promoters[[0,1,2]]).saveas('../../data/EPD_promoter_analysis/TATA/responsive_constitutive_promoters_extended_workspace.bed')
#select only variable promoters
variable_promoters_extended = extended_promoters[extended_promoters[3] == 'highVar']
sorted_variable = variable_promoters_extended.sort_values([0,1])
bed = BedTool.from_dataframe(sorted_variable).saveas('../../data/EPD_promoter_analysis/TATA/variable_promoters_extended.bed')
#make a constitutive only file for completness sake
constitutive_promoters_extended = extended_promoters[extended_promoters[3] == 'housekeeping']
sorted_constitutive = constitutive_promoters_extended.sort_values([0,1])
bed = BedTool.from_dataframe(sorted_constitutive).saveas('../../data/EPD_promoter_analysis/TATA/constitutive_promoters_extended.bed')
log2fold = pd.read_csv('../../data/EPD_promoter_analysis/TATA/TATAlogfold.csv', header=0)
log2fold
#rename
log2fold.Gene_type.replace('Variable','variable', inplace=True)
log2fold.Gene_type.replace('Constitutive','constitutive', inplace=True)
#set style to ticks
sns.set(style="ticks", color_codes=True)
#bar chart, 95% confidence intervals
plot = sns.barplot(x="Gene_type", y="Log2-fold", data=log2fold)
plot.axhline(0, color='black')
plt.xlabel("Gene type")
plt.ylabel("Log2-fold enrichment over background").get_figure().savefig('../../data/plots/TATAbox/log2fold.pdf', format='pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12131440484463417, 'fraction_numerical': 0.01161520897448625, 'mean_word_length': 6.030006978367062, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 17, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '521951', 'license': 'MIT', 'n_tokens_mistral': 3503, 'n_tokens_neox': 3250, 'n_words': 666, 'path': 'src/plotting/TATA_enrichmentold.ipynb', 'repo_name': 'Switham1/PromoterArchitecture', 'size': '10078'}