<jupyter_start><jupyter_text>Read-depth as a CovariateThe output from most genomics experiments (ChIP-seq, ATAC-seq, etc) is a set of reads that can be mapped to a reference genome. Usually, before use these reads are then converted to tracks that encode the number of reads that map to each position in the genome. Potentially, these tracks are subsequently processed to produce other signal values such as fold change versus a control or -log10 p-values from a peak caller. However, all of these outputs can be confounded by a simple value: the read-depth of the experiment. In this notebook we will produce a simple visualizing demonstrating this effect.<jupyter_code>import matplotlib
matplotlib.use('pdf')
%matplotlib inline

import numpy
import pandas
import pyBigWig
import seaborn; seaborn.set_style('whitegrid')
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>We've started by collecting data from the ENCODE Portal measuring H3K27ac in smooth muscle cells that originated from H9 (ENCODE Accession: ENCSR210ZPC). We downloaded the two signal bams and the two control bams for that experiment and re-ran the ENCODE ChIP-seq pipeline using them as input with the default parameters, except that we also subsampled the read-depth of the signal bams to between 10M and 60M reads. After processing, we end up with peak calls and -log10 p-values from MACS for each read-depth. You can download the processed files from Zenodo using the code below.<jupyter_code>!wget -nc https://zenodo.org/record/5518246/files/h3k27ac.10M.bw
!wget -nc https://zenodo.org/record/5518246/files/h3k27ac.20M.bw
!wget -nc https://zenodo.org/record/5518246/files/h3k27ac.30M.bw
!wget -nc https://zenodo.org/record/5518246/files/h3k27ac.40M.bw
!wget -nc https://zenodo.org/record/5518246/files/h3k27ac.50M.bw
!wget -nc https://zenodo.org/record/5518246/files/h3k27ac.60M.bw
!wget -nc https://zenodo.org/record/5518246/files/overlap.optimal_peak.narrowPeak.gz<jupyter_output>File ‘h3k27ac.10M.bw’ already there; not retrieving.

File ‘h3k27ac.20M.bw’ already there; not retrieving.

File ‘h3k27ac.30M.bw’ already there; not retrieving.

File ‘h3k27ac.40M.bw’ already there; not retrieving.

File ‘h3k27ac.50M.bw’ already there; not retrieving.

File ‘h3k27ac.60M.bw’ already there; not retrieving.

File ‘overlap.optimal_peak.narrowPeak.gz’ already there; not retrieving.<jupyter_text>Let's start off by loading chr17 from each of the subsampled data set bigwigs.<jupyter_code>bw_filename = "h3k27ac.{}M.bw"
tracks = []

read_depths = 10, 20, 30, 40, 50, 60
for reads in read_depths:
    filename = bw_filename.format(reads)
    bw = pyBigWig.open(filename)
    track = bw.values('chr17', 0, -1)
    track = numpy.nan_to_num(track)
    tracks.append(track)<jupyter_output><empty_output><jupyter_text>Now, let's go through the peak file and, for each peak, get the maximum signal value in the peak for each of the six subsampled depths. To make things consistent, let's use the peak file from the lowest depth data set. This makes it more likely that each called peak is actually a peak.<jupyter_code>col_names = ['chrom', 'start', 'end']
peak_filename = "overlap.optimal_peak.narrowPeak.gz"
peaks = pandas.read_csv(peak_filename, sep='\t', header=None, names=col_names, usecols=(0, 1, 2))

peak_vals = []
for reads in read_depths:
    filename = bw_filename.format(reads)
    bw = pyBigWig.open(filename)
    peak_vals.append([])
    
    for i, (chrom, start, end) in peaks.iterrows():
        peak_max = numpy.max(bw.values(chrom, start, end)) 
        peak_vals[-1].append(peak_max)<jupyter_output><empty_output><jupyter_text>Finally, let's visualize the data.<jupyter_code>from matplotlib.gridspec import GridSpec

gs = GridSpec(nrows=1, ncols=5)

mid = tracks[-1].argmax()
start, end = mid - 6000, mid + 4000

plt.figure(figsize=(14, 3))
plt.subplot(gs[:3])
plt.title("H3K27ac Signal Value At Differing Depths", fontsize=14)
for i in range(6):
    plt.plot(range(start, end), tracks[i][start:end], label="{}M Reads".format(read_depths[i]))
plt.plot(range(start, end), numpy.zeros(end-start), color='0.7')
plt.legend(fontsize=10)
plt.ylabel("-log10 signal p-value")
plt.xlabel("Genomic Coordinate (chr17)")
plt.ylim(0, 550)
plt.xlim(start, end)

plt.subplot(gs[3:])
for i in range(6):
    x = peak_vals[i]
    plt.plot(sorted(x), numpy.arange(len(x)) / len(x))

plt.title("Distribution of Maximum Signal Value in Peaks", fontsize=14)
plt.xlabel("-log10 signal p-value")
plt.ylabel("CDF")
plt.ylim(0.8, 1.02)
plt.xlim(0, 350)

plt.tight_layout()
plt.savefig("nrg_fig_read-depth.pdf")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11248107289638763, 'fraction_numerical': 0.04261302184728531, 'mean_word_length': 5.020833333333333, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 20, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '572830', 'license': 'MIT', 'n_tokens_mistral': 1762, 'n_tokens_neox': 1592, 'n_words': 480, 'path': 'notebooks/Notebook B - Read-depth as a Covariate.ipynb', 'repo_name': 'shwhalen/ml-pitfalls', 'size': '4623'}