<jupyter_start><jupyter_text>NORMALIZATION 1. get coverage in bin mode2. identify highly covered bins3. apply RLE to get normalizing scaling factors<jupyter_code>def getCoverage(BamPathsList, bs=8000, p=20, e=True, sample_names=None):
    
    cr = crpb.CountReadsPerBin(BamPathsList, 
                               binLength=bs,
                               stepSize = bs,
                               numberOfProcessors=p, 
                               extendReads=e)
    coverage = cr.run()
    
    if sample_names == None:
        sample_names = [os.path.basename(bam) for bam in BamPathsList]
    
    coverage_df = pd.DataFrame(coverage, columns=sample_names)
    
    return coverage_df


def getRLEFactors(coverage_df, broad=False, plot=True):
    
    idx_dict = dict()
    
    for sample in list(coverage_df):
        
        sorted_df = coverage_df[sample].sort_values(ascending=True)
        cumsum_df = sorted_df.cumsum()
        prop_df = cumsum_df.divide(sorted_df.sum())

        X = (np.array(range(prop_df.shape[0]+1))[1:]) / prop_df.shape[0]
        delta = np.fabs(prop_df.subtract(X, axis=0))
        elbow_X = delta.idxmax()
        elbow_Y = prop_df.loc[elbow_X].max()
        
        if broad == False:
            selected_idx = prop_df.loc[prop_df > elbow_Y].index
        else:
            selected_idx = prop_df.loc[prop_df > 0].index
        
        idx_dict[sample] = set(selected_idx)
        
        if plot == True:
            plt.figure()
            plt.plot(X,prop_df)
            plt.plot(X,delta)
            plt.axhline(elbow_Y)
            plt.axvline(X[np.fabs(prop_df.subtract(X, axis=0)).values.argmax()])
    
    common_loci = set.intersection(*[idx_dict[sample] for sample in idx_dict])
    
    enrich_df = coverage_df.loc[common_loci]   

    return 1/rel.RLE(enrich_df)[0] 
SNAKEPIPE_OUTDIR="/data/manke/group/ferrari/PhD_project/reference_datasets/autoRELACS_dataset/output_snakePipes_Drosophila_merged"
BamPathsList_H3K4me3 = glob.glob(os.path.join(SNAKEPIPE_OUTDIR,"filtered_bam/*H3K4me3*.bam"))
c_H3K4me3 = getCoverage(BamPathsList_H3K4me3)
getRLEFactors(c_H3K4me3)/100
SNAKEPIPE_OUTDIR="/data/manke/group/ferrari/PhD_project/reference_datasets/autoRELACS_dataset/output_snakePipes_Drosophila_merged"
BamPathsList_H3K27ac = glob.glob(os.path.join(SNAKEPIPE_OUTDIR,"filtered_bam/*H3K27ac*.bam"))
c_H3K27ac = getCoverage(BamPathsList_H3K27ac)
getRLEFactors(c_H3K27ac)/100
SNAKEPIPE_OUTDIR="/data/manke/group/ferrari/PhD_project/reference_datasets/autoRELACS_dataset/output_snakePipes_Drosophila_merged"
BamPathsList_H3K27me3 = glob.glob(os.path.join(SNAKEPIPE_OUTDIR,"filtered_bam/*H3K27me3*.bam"))
c_H3K27me3 = getCoverage(BamPathsList_H3K27me3)
getRLEFactors(c_H3K27me3, broad=True)/100<jupyter_output><empty_output><jupyter_text>CREATE RLE SCALED COVERAGE TRACKS AND PLOT METAPROFILE ON THE MERGED PEAK SET FOR EACH MARK  H3K4me3<jupyter_code>%%bash

module load deeptools

SNAKEPIPE_OUTDIR="/data/manke/group/ferrari/PhD_project/reference_datasets/autoRELACS_dataset/output_snakePipes_Drosophila_merged"

bamCoverage -b $SNAKEPIPE_OUTDIR/filtered_bam/Biomek_H3K4me3.filtered.bam \
-o ../intermediate_files/Biomek_H3K4me3.filtered.RLEscaled.bw \
--scaleFactor 0.010824 -p 20 -e 

bamCoverage -b $SNAKEPIPE_OUTDIR/filtered_bam/manual_H3K4me3.filtered.bam \
-o ../intermediate_files/manual_H3K4me3.filtered.RLEscaled.bw \
--scaleFactor 0.009239 -p 20 -e 
%%bash

module load deeptools

computeMatrix reference-point -R ../intermediate_files/Merged_Peaks_H3K4me3.bed \
-S ../intermediate_files/Biomek_H3K4me3.filtered.RLEscaled.bw \
../intermediate_files/manual_H3K4me3.filtered.RLEscaled.bw \
-o ../intermediate_files/H3K4me3_merged_plotProfile.mat.gz \
--referencePoint center --missingDataAsZero -p 20 -b 2000 -a 2000
%%bash

module load deeptools

plotProfile -m ../intermediate_files/H3K4me3_merged_plotProfile.mat.gz \
-o ../figures/H3K4me3_merged_plotProfile.pdf --perGroup --averageType median \
--plotHeight 7 --plotWidth 8 --colors darkred coral<jupyter_output><empty_output><jupyter_text>H3K27ac<jupyter_code>%%bash

module load deeptools

SNAKEPIPE_OUTDIR="/data/manke/group/ferrari/PhD_project/reference_datasets/autoRELACS_dataset/output_snakePipes_Drosophila_merged"

bamCoverage -b $SNAKEPIPE_OUTDIR/filtered_bam/Biomek_H3K27ac.filtered.bam \
-o ../intermediate_files/Biomek_H3K27ac.filtered.RLEscaled.bw \
--scaleFactor 0.014424 -p 20 -e 

bamCoverage -b $SNAKEPIPE_OUTDIR/filtered_bam/manual_H3K27ac.filtered.bam \
-o ../intermediate_files/manual_H3K27ac.filtered.RLEscaled.bw \
--scaleFactor 0.006933 -p 20 -e 
%%bash

module load deeptools

computeMatrix reference-point -R ../intermediate_files/Merged_Peaks_H3K27ac.bed \
-S ../intermediate_files/Biomek_H3K27ac.filtered.RLEscaled.bw \
../intermediate_files/manual_H3K27ac.filtered.RLEscaled.bw \
-o ../intermediate_files/H3K27ac_merged_plotProfile.mat.gz \
--referencePoint center --missingDataAsZero -p 20 -b 2000 -a 2000
%%bash

module load deeptools

plotProfile -m ../intermediate_files/H3K27ac_merged_plotProfile.mat.gz \
-o ../figures/H3K27ac_merged_plotProfile.pdf --perGroup --averageType median \
--plotHeight 7 --plotWidth 8 --colors darkgreen limegreen<jupyter_output><empty_output><jupyter_text>H3K27me3<jupyter_code>%%bash

module load deeptools

SNAKEPIPE_OUTDIR="/data/manke/group/ferrari/PhD_project/reference_datasets/autoRELACS_dataset/output_snakePipes_Drosophila_merged"

bamCoverage -b $SNAKEPIPE_OUTDIR/filtered_bam/Biomek_H3K27me3.filtered.bam \
-o ../intermediate_files/Biomek_H3K27me3.filtered.RLEscaled.bw \
--scaleFactor 0.010266 -p 20 -e 

bamCoverage -b $SNAKEPIPE_OUTDIR/filtered_bam/manual_H3K27me3.filtered.bam \
-o ../intermediate_files/manual_H3K27me3.filtered.RLEscaled.bw \
--scaleFactor 0.009741 -p 20 -e 
%%bash

module load deeptools

computeMatrix scale-regions -R ../intermediate_files/Merged_Peaks_H3K27me3.bed \
-S ../intermediate_files/Biomek_H3K27me3.filtered.RLEscaled.bw \
../intermediate_files/manual_H3K27me3.filtered.RLEscaled.bw \
-o ../intermediate_files/H3K27me3_merged_plotProfile.mat.gz \
--missingDataAsZero -p 20 -b 400 -a 400
%%bash

module load deeptools

plotProfile -m ../intermediate_files/H3K27me3_merged_plotProfile.mat.gz \
-o ../figures/H3K27me3_merged_plotProfile.pdf --perGroup --averageType median \
--plotHeight 7 --plotWidth 8 --colors black dimgray<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10834507593549397, 'fraction_numerical': 0.04305620792234226, 'mean_word_length': 4.0338849487785655, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '549608', 'license': 'MIT', 'n_tokens_mistral': 2732, 'n_tokens_neox': 2469, 'n_words': 331, 'path': 'figure_2/notebooks/2b_metaprofiles.ipynb', 'repo_name': 'FrancescoFerrari88/AutoRELACS', 'size': '6387'}