DNase analysis

Contact: oursu@stanford.edu

chromovar3d.stanford.edu/DNase/ - contains all analysis, with the following subdirectories:

alignments
- alignments/liftOver/ - reads from Degner et al. lifted over to hg19 (tagAlign format)
- alignments/subsampled_reads/ - reads subsampled from merged replicates, at 35 million reads per individual

signal
signal/bigwig/ - bigwig files of signal for DNase
signal/align2rawsignal - signal for DNase in .mat format (for use with extractSignal)

peaks
peaks/peaks_per_individual/ - peaks called on each individual
peaks/mergedPeaks_TrimFromSummit - peaks in each individual were trimmed to be 100bp long centered on the summit, then these peaks were merged across individuals, which is contained in this directory
final list of peaks (before selecting the top 200 000 peaks): peaks/mergedPeaks_TrimFromSummit/DNase_TrimFromSummit.mergeBedremoveBlacklist_Log10PvalueThreshold_5.gz

data matrix
data_matrix/ - Mean signal extracted in the peaks from peaks/mergedPeaks_TrimFromSummit/
(For easy access, this data matrix can also be accessed as: chromovar3d.stanford.edu/QTLs/uncorrectedSignal/DNase_removeBlacklist_Log10PvalueThreshold_5_DATA_MATRIX.gz)

Methods

We downloaded DNaseI data for 70 individuals from the (http://www.ncbi.nlm.nih.gov/pubmed/22307276), as aligned reads (hg18 assembly) from: http://eqtl.uchicago.edu/dsQTL_data/MAPPED_READS/. We converted read coordinates to hg19 using LiftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver). For further analysis, we only considered the reads mapping to chromosomes 1-22, X and Y, since for chrM we expected no histone mark signal.

To call DNase peaks, we first shifted reads by 75 base pairs (75bp to the left for reads on plus strand and 75bp to the right for reads on the minus strand), leading to expected peak sizes of 150bp, and then called peaks using MACS2 on pooled replicates (subsampled to 35M reads to ensure comparable sequencing depths across samples), with the following parameters: --nomodel option, --shiftsize=75 --pvalue 1e-2. We trimmed peaks in each individual to be 100bp long, centered on the peak summit provided by MACS (Zhang et al. 2008) (column 10), and removed peaks overlapping blacklisted regions of the genome. We merged the trimmed peaks across individuals, and filtered peaks whose minimum p-value across all individuals was higher than 1e-5. This filters out peaks that have bad quality in all individuals.

To generate the signal used for QTL analysis, for each peak in the final set of peaks, we extract the mean signal in each individual, as computed by align2rawsignal (parameters: kernel (k)=epanechnikov, fragment length (l)=150, smoothing window (w)=150, normFlag(n)=5, mapFilter (f)=0) followed by extractSignal.

Since previous studies have reported around 200 000 DHS peaks in each cell type, we decided to use the top 250000 peaks (with regard to mean signal across individuals) for our QTL study.