Skip to content

Preprocessing

Jump to bottom
Vivekanandan Ramalingam edited this page Mar 19, 2024 · 12 revisions

1. Experimental dataset

For this tutorial we'll use experimental CHIP-seq data, for the transcription factor CTCF in the K562 cell line, which is available on the ENCODE data portal. There are 5 such experiments that we find in ENCODE, you can see them listed here CHIP-seq CTCF K562 . We'll restrict ourselves to one experiment ENCSR000EGM

Download the .bam files for the two replicates shown below in the image.

replicate bams

The two replicates are isogenic replicates (biological). A more detailed explanation of the various types of replicates can be found here.

Links to the replicate bam files provided below.

ENCFF198CVB

ENCFF488CXC

wget https://www.encodeproject.org/files/ENCFF198CVB/@@download/ENCFF198CVB.bam -O rep1.bam
wget https://www.encodeproject.org/files/ENCFF488CXC/@@download/ENCFF488CXC.bam -O rep2.bam

Now download the control for the experiment, which is available here ENCSR000EHI

Download the bam file shown in the image below.

control bam

Link provided below

ENCFF023NGN

wget https://www.encodeproject.org/files/ENCFF023NGN/@@download/ENCFF023NGN.bam -O control.bam

Finally, download the reference files. In the example below, some preprocessing is required to filter out unwanted chromosomes from the hg38.chrom.sizes file. Additionally, the blacklist file shown is specific to hg38, and should be replaced with a genome-specific blacklist if alternative genomes are used.

Available Blacklists:

For those interested in using the blacklists, a current version for dm3, dm6, ce10, ce11, mm10, hg19, and hg38 are available in the lists/ folder at https://github.com/Boyle-Lab/Blacklist/

Please cite:

Amemiya, H.M., Kundaje, A. & Boyle, A.P. The ENCODE Blacklist: Identification of Problematic Regions of the Genome. Sci Rep 9, 9354 (2019). https://doi.org/10.1038/s41598-019-45839-z

# download genome refrence
wget https://www.encodeproject.org/files/GRCh38_no_alt_analysis_set_GCA_000001405.15/@@download/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta.gz \
-O hg38.genome.fa.gz | gunzip

# index genome reference
samtools faidx hg38.genome.fa

# download chrom sizes
wget https://www.encodeproject.org/files/GRCh38_EBV.chrom.sizes/@@download/GRCh38_EBV.chrom.sizes.tsv

# exclude alt contigs and chrEBV
grep -v -e '_' -e 'chrEBV' GRCh38_EBV.chrom.sizes.tsv > hg38.chrom.sizes
rm GRCh38_EBV.chrom.sizes.tsv

# make file with chromosomes only
awk '{print $1}' hg38.chrom.sizes > chroms.txt

# download blacklist
wget https://www.encodeproject.org/files/ENCFF356LFX/@@download/ENCFF356LFX.bed.gz -O blacklist.bed.gz
gunzip blacklist.bed.gz

1.1 Preprocessing steps to generate bigwig counts tracks

For the following steps you will need samtools bamtools and bedGraphToBigWig, which are not installed as part of this repository.

The tools can be installed via the links below or using conda.

samtools

bamtools

bedGraphToBigWig (Linux 64-bit)

bedGraphToBigWig (Mac OSX 10.14.6)

conda install -c conda-forge -c bioconda samtools=1.19.2 bamtools=2.5.2 ucsc-bedgraphtobigwig=445
1.1.1 Merge the two replicates and create and index
samtools merge -f merged.bam rep1.bam rep2.bam
samtools index merged.bam

We will also index the control BAM file

samtools index control.bam
1.1.2 Create bigwig files using bedtools via intermediate bedGraph files

In addition to creating the bigwig files, at this step we will filter the bam files to keep only the chromosomes that we want to use in the model. In the example shown below, we do this using samtools view and the hg38.chrom.sizes reference file.

Experiment

# get coverage of 5’ positions of the plus strand
samtools view -b merged.bam $(cut -f 1 hg38.chrom.sizes) | \
	bedtools genomecov -5 -bg -strand + -ibam stdin | \
	sort -k1,1 -k2,2n > plus.bedGraph

# get coverage of 5’ positions of the minus strand
samtools view -b merged.bam $(cut -f 1 hg38.chrom.sizes) | \
        bedtools genomecov -5 -bg -strand - -ibam stdin | \
        sort -k1,1 -k2,2n > minus.bedGraph

# Convert bedGraph files to bigWig files
bedGraphToBigWig plus.bedGraph hg38.chrom.sizes plus.bw
bedGraphToBigWig minus.bedGraph hg38.chrom.sizes minus.bw

Control

# get coverage of 5’ positions of the control plus strand
samtools view -b control.bam $(cut -f 1 hg38.chrom.sizes) | \
        bedtools genomecov -5 -bg -strand + -ibam stdin | \
        sort -k1,1 -k2,2n > control_plus.bedGraph

# get coverage of 5' positions of the control minus strand
samtools view -b control.bam $(cut -f 1 hg38.chrom.sizes) | \
        bedtools genomecov -5 -bg -strand - -ibam stdin | \
        sort -k1,1 -k2,2n > control_minus.bedGraph

# Convert bedGraph files to bigWig files
bedGraphToBigWig control_plus.bedGraph hg38.chrom.sizes control_plus.bw
bedGraphToBigWig control_minus.bedGraph hg38.chrom.sizes control_minus.bw

1.2 Identify peaks

For the purposes of this tutorial we will use the optimal IDR thresholded peaks that are already available in the ENCODE data portal. We will use the the narrowPeak files that are in BED6+4 format. Explanation of what each of the 10 fields means can be found here. Currently, only this format is supported but in the future support for more formats will be added.

See image below that shows the file listed in the ENCODE data portal

Download the file:

wget https://www.encodeproject.org/files/ENCFF396BZQ/@@download/ENCFF396BZQ.bed.gz -O peaks.bed.gz
gunzip peaks.bed.gz

1.3 Organize you data

At this point, we suggest creating a directory structure to store the data, models, predictions, metrics, importance scores, discovered motifs, plots & visualizations etc. that will make it easier for you to organize and maintain your work. Let's start by creating a parent directory for the experiment and moving the bigwig files and peaks file from section 1.1 & 1.2 to a data directory

mkdir ENCSR000EGM
mkdir ENCSR000EGM/data
mv *.bw ENCSR000EGM/data
mv peaks.bed ENCSR000EGM/data

Once this is done, your directory hierarchy should resemble this

Note that the relative paths in subsequent pipeline steps assume that the current working directory is the one immediately above the project directory. For example, if you are in the ENCSR000EGM folder, navigate up one level with cd ...

Make a reference directory at the same level as the ENCSR000EGM experiment directory. In the reference directory we will place 4 files: the genome reference hg38.genome.fa, the indexed reference hg38.genome.fai, the chromosome sizes file hg38.chrom.sizes and one text file that contains a list of chromosomes we care about chroms.txt. blacklist.bed contains all the blacklist regions with artifact peaks.

mkdir ENCSR000EGM/reference 
mv hg38.genome.fa* ENCSR000EGM/reference
mv hg38.chrom.sizes ENCSR000EGM/refence
mv chroms.txt ENCSR000EGM/reference
mv blacklist.bed ENCSR000EGM/reference

The directory structure should look like this (TODO: update this image):

1.4 Outlier removal

Filter the peaks file for outliers. Peaks that meet either of two criteria are removed:

(1) Overlap with regions identified in blacklist.bed. (2) Number of reads in the peak is in the --quantile quantile.

First prepare a file input_outliers.json as shown below:

{
    "0": {
        "signal": {
            "source": ["ENCSR000EGM/data/plus.bw",
                       "ENCSR000EGM/data/minus.bw"]
        },
        "loci": {
            "source": ["ENCSR000EGM/data/peaks.bed"]
        },
        "bias": {
            "source": ["ENCSR000EGM/data/control_plus.bw",
                       "ENCSR000EGM/data/control_minus.bw"],
            "smoothing": [null, null]
        }
    }
}

Next, run the following command:

bpnet-outliers \
    --input-data input_outliers.json  \
    --quantile 0.99 \
    --quantile-value-scale-factor 1.2 \
    --task 0 \
    --chrom-sizes ENCSR000EGM/reference/hg38.chrom.sizes \
    --chroms $(paste -s -d ' ' ENCSR000EGM/reference/chroms.txt) \
    --sequence-len 1000 \
    --blacklist ENCSR000EGM/reference/blacklist.bed \
    --global-sample-weight 1.0 \
    --output-bed ENCSR000EGM/data/peaks_inliers.bed

Clone this wiki locally