ChIP-seq analysis part1

sofware installation and data source

install MACS2: 06/11/2015 sudo pip -H install MACS2
macs2 version:2.1.0.20150420

install NCIS to estimate the scaling factor:
Download the package here:
wget http://www.biomedcentral.com/content/supplementary/1471-2105-13-199-s2.gz

It is a .gz file, you can open your Rstudio and install it:
install.packages("~/NCIS/1471-2105-13-199-s2.gz", repos = NULL, type="source")
For help:
library(NCIS); ?NCIS

see a post here using NCIS before MACS peak calling:
Adding a custom normalization to MACS.
when use NCIS, one needs to deduplicate the reads first. see a post on the MACS goolge group. Then the estimated scaling factor between ChIP and input control will be feeded into MACS2 with the --ratio flag:
macs2 callpeak -t ChIP.bam -c Control.bam --broad -g hs --broad-cutoff 0.1 --ration 1.4

Data source

downloaded the data from the ENNCODE project on 06/11/2015 here for MCF7 and panc1 cells. H3k27ac ChIP-seq, duplicates for each cell line. I choose these data sets because the SYDH Histone tracks contain input DNAs. The Broad Histione tracks do not have input DNA sequenced.

MCF7 cells wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep1.bam
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep1.bam.bai
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep2.bam.bai
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep2.bam
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7H3k27acUcdPk.narrowPeak.gz

input DNA wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7InputUcdAlnRep1.bam
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7InputUcdAlnRep1.bam.bai
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7InputUcdAlnRep2.bam
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistoneMcf7InputUcdAlnRep2.bam.bai

Panc1 cells
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam.bai
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep2.bam wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep2.bam.bai
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistonePanc1H3k27acUcdPk.narrowPeak.gz

input DNA
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeSydhHistone/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam.bai

project layout

It is always good to organize your project. For me, the layout is like this:

diff_ChIP_test is the project main folder and it contains four sub-folders: data, doc , results and scripts. Data are downloaded into the data folder, scripts are in the scripts folder and the output from the scripts are in the resultsfolder.
Further Reading: A Quick Guide to Organizing Computational Biology Projects

Is macs2 OK for broad peaks such as H3K27ac?

a discussion here
Yet another note from a user argues that MACS2 can be used for broad peak calling, but by choosing correctarguments:

While I can't speak to much of the first paragraph of your message, I wanted to let you know that when it comes to histone modification analysis, using the --call-supeaks option in the command line has proven to be a phenomenally successful way to use MACS to handle modifications (right now, I'm working on H3K4me1) which have both sharp defined peaks but also very broad "peaks" which are more like domains of H3K4me1 signals clustered together. This obviously created many situations, in part because MACS will automatically combine peaks which are separated by a distance of 10bp or less into one called peak, where artefactually broad regions were called (like 30kb peaks... yeah right!).

So many people have argued, including in a recent review in Nat Immunology about ChIP-seq and peak-calling, that modifications with broad distribution (even H3K27Ac can be broad in our hands) should not be analyzed with MACS; using various biological end-points as testing I have found that this is not true and that MACS outperforms ZINBA (specifically designed for broad histone modifications) WHEN the subpeaks function is called for detection of histone peaks with broad, narrow, or both types of signal distribution.

It seems to me that MACS2 has evloved a lot to deal with the broad peaks compared with the widely used MACS14. Although other tools such as SICER are designed sepcifically for histone modifications, I am still going to use MACS2 for H3K27ac ChIP-seq peak calling.

Further Reading:
Practical Guidelines for the Comprehensive Analysis of ChIP-seq Data
ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia

starting pilot analysis

I am going to use one single ChIP bam file and one input file to do some initial testing with different parameters of MACS2.

--broad

When this flag is on, MACS will try to composite broad regions in BED12 ( a gene-model-like format ) by putting nearby highly enriched regions into a broad region with loose cutoff. The broad region is controlled by another cutoff through --broad-cutoff. The maximum length of broad region length is 4 times of d from MACS. DEFAULT: False

call peaks with macs2 using --broad, building model: macs2 callpeak -t ../data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam -c ../data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam --broad -g hs --broad-cutoff 0.1 -n panc1H3k27acRep1 --outdir panc1H3k27acRep1_with_model_broad

call peaks with macs2 using --broad, bypass the model: macs2 callpeak -t ../data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam -c ../data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam --broad -g hs --broad-cutoff 0.1 -n panc1H3k27acRep1 --outdir panc1H3k27acRep1_without_model_broad --nomodel --extsize 146

-nomodel and --extsize 146 tell MACS2 use 146bp as fragment size to pileup sequencing reads.

For model building, a NAME_model.r is an R script in the output which you can use to produce a PDF image about the model based on your data. Load it to R by:

$ Rscript NAME_model.r

The model will estimate the d

A note from Tao Liu:

If the d is not small ~ < 2*tag size (for those tag size < 50bp), and the model image in PDF shows clean bimodal shape, d may be good. And several bp differences on d shouldn't affect the peak detection on general transcription factor ChIP-seq much.

However, for Pol2 or histone marks, things may be different. Pol2 is moving so it's not appropriate to say there is a fixed fragment size. I don't know the correct answer. For histone mark ChIP-seq, since they would have a underlying characteristic 147bp resolution for a nucleosome size, you can simply skip model building and use "--shiftsize 74 --nomodel" instead. Also if you want, you can try other software like SICER and NPS.

--extsize EXTSIZE The arbitrary extension size in bp. When nomodel is true, MACS will use this value as fragment size to extend each read towards 3' end, then pile them up. It's exactly twice the number of obsolete SHIFTSIZE. In previous language, each read is moved 5'->3' direction to middle of fragment by 1/2 d, then extended to both direction with 1/2 d. This is equivalent to say each read is extended towards 5'->3' into a d size fragment. DEFAULT: 200. EXTSIZE and SHIFT can be combined when necessary. Check SHIFT option.

call peaks with macs2 not using --broad, bypass the model, produce bedgraph file -B, and --call-summit macs2 callpeak -t ../data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam -c ../data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam -g hs -q 0.01 -n panc1H3k27acRep1_regular --call-summit -B --outdir panc1H3k27acRep1_without_model_regular --nomodel --extsize 146

short note fom me
The bedgraph file generated by macs2 is very huge (1Gb for this particular case), because it contains decimals(?). If you want to visualize it in IGV, you need to get a TDF file first. Remember to change the suffix .bdg to .bedgraph, otherwise IGV will not recoginize the file format.

I personally like to covert the bam files directly to bigwig files using deeptools. Using 10bp as a bin size, I get a bigwig file of 205Mb and you can directly load it into IGV.
bamCoverage -b ../data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam --normalizeTo1x 2451960000 --missingDataAsZero yes --binSize 10 -o panc1_H3k27acRep1_deeptool_normalized.bw

Results from pilot analysis

1. peak numbers

First, let's look at how many peaks are produced by three different settings of macs2 arguments.
1. use --broad, build model
it produced 54223 broadpeaks.
2. use --broad, bypass model
it produced 66288 broadpeaks.
3. regular peak calling will produce a summit.bed file and narrow peaks
it produced 182569 narrowpeaks.

using --broad definetly improve the identification of peaks(or more appropriately:enriched regions). You can find the narrow peaks in the gappedPeak file:

NAME_peaks.gappedPeak is in BED12+3 format which contains both the broad region and narrow peaks. The 5th column is 10*-log10qvalue, to be more compatible to show grey levels on UCSC browser. Tht 7th is the start of the first narrow peak in the region, and the 8th column is the end. The 9th column should be RGB color key, however, we keep 0 here to use the default color, so change it if you want. The 10th column tells how many blocks including the starting 1bp and ending 1bp of broad regions. The 11th column shows the length of each blocks, and 12th for the starts of each blocks. 13th: fold-change, 14th: -log10pvalue, 15th: -log10qvalue. The file can be loaded directly to UCSC genome browser.

The narrowPeak file downloaded from the ENCODE website contains 73953 narrow peaks. zcat ../data/wgEncodeSydhHistonePanc1H3k27acUcdPk.narrowPeak.gz| wc -l 

Take home messages for now:
1. Using different tools to call peaks will produce different number of peaks.
2. Using the same tool with different settings will produce different number of peaks.
There is no consensus to use which arguments. It depends on your data type and purpose. The bottom line is that you need to make sense of the data and find biological siginifcances from there.

As long as you document how you did the analysis (so that other people can reproduce your work) and can convince people you are doing it reasonably, you are fine.

I will stick to use these argements for all my subsequent analysis:
macs2 callpeak -t ../data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam -c ../data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam --broad -g hs --broad-cutoff 0.1 -n panc1H3k27acRep1 --outdir panc1H3k27acRep1_without_model_broad --nomodel --extsize 146

2. peak quality in terms of qvalues (FDR) and overlapping.

using bedtools for overlapping:

bedtools intersect -a panc1H3k27acRep1_with_model_broad/panc1H3k27acRep1_peaks.broadPeak -b panc1H3k27acRep1_without_model_broad/panc1H3k27acRep1_peaks.broadPeak -wa | cut -f1-3 | sort | uniq | wc -l 45228

bedtools intersect -a panc1H3k27acRep1_with_model_broad/panc1H3k27acRep1_peaks.broadPeak -b panc1H3k27acRep1_without_model_broad/panc1H3k27acRep1_peaks.broadPeak -wb | cut -f1-3 | sort | uniq | wc -l 61507
largely they are overlapped with each other.

compare with the narrowpeaks downloaded from ENCODE:
no-model:
bedtools intersect -a panc1H3k27acRep1_without_model_broad/panc1H3k27acRep1_peaks.broadPeak -b ../data/wgEncodeSydhHistonePanc1H3k27acUcdPk.narrowPeak -wa | cut -f1-3 | sort | uniq | wc -l 35722

build model
bedtools intersect -a panc1H3k27acRep1_with_model_broad/panc1H3k27acRep1_peaks.broadPeak -b ../data/wgEncodeSydhHistonePanc1H3k27acUcdPk.narrowPeak -wa | cut -f1-3 | sort | uniq | wc -l 28943

it seems that without model building is giving more concordant peaks compared with the ENCODE narrowpeak.

check FDRs

columns of the broadPeak file:
7th: fold-change, 8th: -log10pvalue, 9th: -log10qvalue
How many peaks have a FDR of 0.01 and fold-change of 2:
peaks with no-model building
cat panc1H3k27acRep1_without_model_broad/panc1H3k27acRep1_peaks.broadPeak | awk '$9 >2 && $7 >2' | wc -l 41136

peaks after building model
cat panc1H3k27acRep1_with_model_broad/panc1H3k27acRep1_peaks.broadPeak | awk '$9 >2 && $7 >2' | wc -l 35487

it seems that without model building gives more confident peaks

I did some exploratory analysis using R and published at Rpub

mannual visualization of called peaks.

call peaks for all the samples.

calculate scaling factor for each ChIP bam file using NCIS library

By default "Larger dataset will be scaled towards smaller dataset" for macs2, one can call peaks with macs2 using --ratio flag. NCIS library needs the ChIP bam and input control bam file as two arguments. I need to write a R script and excute on the command line by using Rscript:
how to get commands line arguments for R

## calculate the scaling factor using NCIS for ChIP-seq data
# see links https://groups.google.com/forum/#!searchin/macs-announcement/NCIS/macs-announcement/0EF4cQF09FI/2-zlu2rqfOkJ
# http://searchvoidstar.tumblr.com/post/52594053877/adding-a-custom-normalization-to-macs
# Ming Tang 06/15/2015

library(NCIS)

library(ShortRead)

options(echo=TRUE) # set to FALSE if you not  want see commands in output 
args <- commandArgs(trailingOnly = TRUE)
print(args)
# trailingOnly=TRUE means that only your arguments are returned, check:
# print(commandsArgs(trailingOnly=FALSE))

ChIP_bam<- args[1]
input_control_bam<- args[2]

# NCIS usese the Aligned Reads object from the shortRead package, however, it is recommended
# to use GenomicAignments package to read in the bam files
# ga_ChIP<- readGAlignments(ChIP_bam)
# ga_input<-readGAlignments(input_control_bam)
# However, the resulting GenomicAlignment object is not recognized by NCIS.
# I have to use the legacy readAligned function from the ShortRead package.
# it takes around 15mins to finish

ga_ChIP<- readAligned(ChIP_bam, type="BAM")
ga_input<-readAligned(input_control_bam, type="BAM")

res<- NCIS(ga_ChIP, ga_input, data.type="AlignedRead")
res
res$est
res$r.seq.depth

To use it, save the R script as NCIS_scaling.r and on terminal:
Rscript ChIP.bam control.bam
It will output the scaling factor of ChIP/control. Note that the MACS2 flag --ratio is also for ChIP/control.

 --ratio RATIO         When set, use a custom scaling ratio of ChIP/control
                        (e.g. calculated using NCIS) for linear scaling.
                        DEFAULT: ingore

scaling factor for all four samples

Rscript NCIS_scaling.r ../data/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep1.bam ../data/wgEncodeSydhHistoneMcf7InputUcdAlnRep1.bam
Rscript NCIS_scaling.r ../data/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep2.bam ../data/wgEncodeSydhHistoneMcf7InputUcdAlnRep1.bam
Rscript NCIS_scaling.r ../data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam ../data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam
Rscript NCIS_scaling.r ../data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep2.bam ../data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam

The Scaling factor and the sequencing depth ratio are all for ChIP/control

file_name	scaling factor	seq depth ratio
Mcf7Rep1	1.761996	2.091474
Mcf7Rep2	1.709105	1.917619
panc1Rep1	0.5861799	1.110217
panc1Rep2	0.5160342	0.8233171

we can get a rough idea of the size of each bam file:

ls -sh data/*
4.0K data/README.md 1.3G data/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep1.bam 6.2M data/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep1.bam.bai 1.2G data/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep2.bam 6.1M data/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep2.bam.bai 468K data/wgEncodeSydhHistoneMcf7H3k27acUcdPk.narrowPeak.gz773M data/wgEncodeSydhHistoneMcf7InputUcdAlnRep1.bam 5.9M data/wgEncodeSydhHistoneMcf7InputUcdAlnRep1.bam.bai718M data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam 6.1M data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam.bai 522M data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep2.bam5.9M data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep2.bam.bai 4.4M data/wgEncodeSydhHistonePanc1H3k27acUcdPk.narrowPeak 678M data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam 6.1M data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam.bai

MACS2 peak calling with --ratio
macs2 callpeak -t data/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep1.bam -c data/wgEncodeSydhHistoneMcf7InputUcdAlnRep1.bam --broad -g hs --broad-cutoff 0.1 -n Mcf7H3k27acUcdAlnRep1_ratio --outdir results/Mcf7H3k27acUcdAlnRep1_ratio --nomodel --extsize 146 --ratio 1.761996
macs2 callpeak -t data/wgEncodeSydhHistoneMcf7H3k27acUcdAlnRep2.bam -c data/wgEncodeSydhHistoneMcf7InputUcdAlnRep1.bam --broad -g hs --broad-cutoff 0.1 -n Mcf7H3k27acUcdAlnRep2_ratio --outdir results/Mcf7H3k27acUcdAlnRep2_ratio --nomodel --extsize 146 --ratio 1.709105
macs2 callpeak -t data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep1.bam -c data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam --broad -g hs --broad-cutoff 0.1 -n Panc17H3k27acUcdAlnRep1_ratio --outdir results/Panc1H3k27acUcdAlnRep1_ratio --nomodel --extsize 146 --ratio 0.5861799
macs2 callpeak -t data/wgEncodeSydhHistonePanc1H3k27acUcdAlnRep2.bam -c data/wgEncodeSydhHistonePanc1InputUcdAlnRep1.bam --broad -g hs --broad-cutoff 0.1 -n Panc17H3k27acUcdAlnRep2_ratio --outdir results/Panc1H3k27acUcdAlnRep2_ratio --nomodel --extsize 146 --ratio 0.5861799

MACS2 peak calling without --ratio

To bulk process the bam files downloaded from ENCODE,one can write a bash script.

#! /bin/bash

set -e
set -u
set -o pipefail -o errexit -o nounset

# we loop for the ChIP bam files
for bam in ../data/*H3k27ac*bam
do 
    # strip out only the meaningful filename to be used for output
    file_name=$(echo "$bam" | sed -E "s/..\/data\/wg.+Histone(.+)(H3k27ac.+).bam/\1\2/")

    # need to retain the ../data/ path. it could be simply: sed -E "s/H3k27ac/Input/" if 
    # every bam file has a input control
    input_control=$(echo "$bam" | sed -E "s/(wg.+Histone)(.+)(H3k27ac.+).bam/\1\2InputUcdAlnRep1.bam/")

    echo "processing ${file_name} bam file"
    echo "the input control file is ${input_control}"
    echo "calling peaks with macs2"
    macs2 callpeak -t "$bam" -c "${input_control}"  --broad -g hs --broad-cutoff 0.1 -n "${file_name}" --outdir ../results/"${file_name}"  --nomodel --extsize 146

done

The sed regular expression caused me some headache. for + operator to function in sed, one needs to turn on the -Eflag for extended regular expression on macOS, or -r on GNU sed.
To caputure part of the pattern, use \1, and the parentheses do not need to be escaped:
If you wanted to keep the first word of a line, and delete the rest of the line: sed 's/\([a-z]*\).*/\1'
or turn on the -E flag:
sed -E 's/([a-z]*).*/\1'
see a very good tutorial on sed.

After execute the bash script, 4 folders are created in the resultsfolder:Mcf7H3k27acUcdAlnRep1,Mcf7H3k27acUcdAlnRep2,Panc1H3k27acUcdAlnRep1 and Panc1H3k27acUcdAlnRep2

Peak number with and without --ratio

file_name	with --ratio	without --ratio	overlapping
Mcf7Rep1	25744	18572	18572
Mcf7Rep2	17564	12257	12256
panc1Rep1	348968	66288	66282
panc1Rep2	71033	54312	54229

It looks like that including the --ratio will generally increase the peak number. Especially for Pan1Rep1, it has 5 times more peaks with --ratio than without --ratio.
It will be interesting to check the peak quality of the increasing number of peaks after adding --ratio. I mannually checked several peaks and found that the newly found peaks are mostly very weak ones. I decided to do my subsequent analysis with the peaks got from the MACS2 without --ratio.

filter peaks from the blacklists

blacklists

Functional genomics experiments based on next-gen sequencing (e.g. ChIP-seq, MNase-seq, DNase-seq, FAIRE-seq) that measure biochemical activity of various elements in the genome often produce artifact signal in certain regions of the genome. It is important to keep track of and filter artifact regions that tend to show artificially high signal (excessive unstructured anomalous reads mapping). Below is a list of comprehensive empirical blacklists identified by the ENCODE and modENCODE consortia. Note that these blacklists were empirically derived from large compendia of data using a combination of automated heuristics and manual curation. These blacklists are applicable to functional genomic data based on short-read sequencing (20-100bp reads). These are not directly applicable to RNA-seq or any other transcriptome data types. The blacklisted regions typically appear u niquely mappable so simple mappability filters do not remove them. These regions are often found at specific types of repeats such as centromeres, telomeres and satellite repeats. It is especially important to remove these regions that computing measures of similarity such as Pearson correlation between genome-wide tracks that are especially affected by outliers.

download the hg19 blacklist into the data folder by wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/wgEncodeDacMapabilityConsensusExcludable.bed.gz

creat a new folder broad_peaks_all inside the results folder, and copy all the broad peaks to this folder.
ls ./broad_peaks_all
Mcf7H3k27acUcdAlnRep1_peaks.broadPeak Panc1H3k27acUcdAlnRep1_peaks.broadPeakMcf7H3k27acUcdAlnRep2_peaks.broadPeak Panc1H3k27acUcdAlnRep2_peaks.broadPeak

#! /bin/bash

## filter the broad peaks from MACS output against the blacklist regions
## bedtools intersect -v

set -e
set -u
set -o pipefail -o errexit -o nounset


for peak in ../results/broad_peaks_all/*broadPeak
do
    file_name=$(basename $peak .broadPeak)
    bedtools intersect -a "$peak" -b ../data/wgEncodeDacMapabilityConsensusExcludable.bed -v \
    > ../results/broad_peaks_all/"${file_name}".filtered.bed
done

save the bash script as filter_blacklist.sh and run it at terminal in the script folder ./fiter_blacklist.sh.

chmod u+x filter_blacklist.sh

After filtering:
wc -l /broad_peaks_all/*
18572 Mcf7H3k27acUcdAlnRep1_peaks.broadPeak 18548 Mcf7H3k27acUcdAlnRep1_peaks.filtered.bed 12257 Mcf7H3k27acUcdAlnRep2_peaks.broadPeak 12239 Mcf7H3k27acUcdAlnRep2_peaks.filtered.bed 66288 Panc1H3k27acUcdAlnRep1_peaks.broadPeak 66248 Panc1H3k27acUcdAlnRep1_peaks.filtered.bed 54312 Panc1H3k27acUcdAlnRep2_peaks.broadPeak 54235 Panc1H3k27acUcdAlnRep2_peaks.filtered.bed
Then remove all the broadpeaks:
rm *broadPeak

merge peaks

After get the filtered peaks, I need to merge all the peaks into a superset that contains all the peaks. I will first select out the peaks that overlap with each other between biological replicates, and then merge all four of them using bedtools.
bedtools intersect -a Mcf7H3k27acUcdAlnRep1_peaks.filtered.bed -b Mcf7H3k27acUcdAlnRep2_peaks.filtered.bed -wa | cut -f1-3 | sort | uniq > Mcf7Rep1_peaks.bed
bedtools intersect -a Mcf7H3k27acUcdAlnRep1_peaks.filtered.bed -b Mcf7H3k27acUcdAlnRep2_peaks.filtered.bed -wb | cut -f1-3 | sort | uniq > Mcf7Rep2_peaks.bed
bedtools intersect -a Panc1H3k27acUcdAlnRep1_peaks.filtered.bed -b Panc1H3k27acUcdAlnRep2_peaks.filtered.bed -wa | cut -f1-3 | sort | uniq > Panc1Rep1_peaks.bed
bedtools intersect -a Panc1H3k27acUcdAlnRep1_peaks.filtered.bed -b Panc1H3k27acUcdAlnRep2_peaks.filtered.bed -wb | cut -f1-3 | sort | uniq > Panc1Rep2_peaks.bed

wc -l *
18548 Mcf7H3k27acUcdAlnRep1_peaks.filtered.bed 12239 Mcf7H3k27acUcdAlnRep2_peaks.filtered.bed 9679 Mcf7Rep1_peaks.bed
11319 Mcf7Rep2_peaks.bed 66248 Panc1H3k27acUcdAlnRep1_peaks.filtered.bed 54235 Panc1H3k27acUcdAlnRep2_peaks.filtered.bed 35218 Panc1Rep1_peaks.bed
44358 Panc1Rep2_peaks.bed

rm *filtered*
cat *bed | sort -k1,1 -k2,2n | bedtools merge | tee merge.bed | wc -l
39046
we got a final superset containing 39046 peaks.

change the mac terminal color scheme and use oh my zsh!

I just got my lab mac pro. I wanted to change the color scheme of the default terminal, which is just black and white.

I googled around, and found solarized .

you can download it and unzip it.
Then, open your terminal (you can search "terminal" at the upper right corner)

go to the preferences:

and import the color scheme for terminals in the folder osx-terminal.app-colors-solarized, which is in the unzipped the folder solarized. I preferred the Solarized Dark ansi.terminal , and set it as Default.

Now,  download oh my zsh! by:

$ curl -L https://raw.github.com/robbyrussell/oh-my-zsh/master/tools/install.sh | sh

You will have to install git first for the upper command to work.
Now, when you fire up your terminal, it looks much more prettier! (there are many other schemes for oh my zsh, I found the default is good. You can change it by modifying the .zshrc file in your home directory.)

Using Awk to format the gene expression file to gct format for Gene Set Enrichment Analysis (GSEA)

After one runs microarray or RNA-seq analysis, usually he would do a Gene Set Enrichment Analysis (GSEA) analysis.  There are many tools to use. One of the most commonly used one is GSEA developed in Broad Institute.

It requires four data files to be loaded:
1. Expression dataset in res, gct, pcl or txt format
2. Phenotype labels in cls format
3. Gene sets in gmx or gmt formt
4. Chip annotations

The first impression of mine is that: Oh my, why there are so many different formats? Yes, after merging into the computational biology field for a while, I find that most of the time I spend is on data formatting. That's in consistence with many others' experiences.

Well, for this post, I will specifically show you how to format gene expression data file output from affy (for microarray) to gct format using awk. For RNA-seq data, you can do it similarly for DESeq2 and EdgR outputs (using normalized counts).

Let's look at the expression file output by affy:

# R code
library(affy)
## read in the data
Data<- ReadAffy()
## RMA normalization and get the eset (expressionSet) object

eset<- rma(Data)

e<- exprs(eset)
write.table( e, "raw_expression.txt", row.names=F, quote=F, sep="\t")

The file we have:

The required file format:

we see that the first column is the probe name and the other columns are expression values for different samples. The first problem is that the first line is one grid off; the first column should have a name "Name". In addition, we need to add two lines, and we need to add a dummy column in the second column. We will fix it step by step:


Now, we have the desired format:

You can certainly open the file in excel and edit it very easily. The file is only several MB big. However, when you have a file that is several hundred MB or several GB big, you can not open it with excel. Avoiding using excel is my ultimate goal, although it comes very handy for small data sets. Using excel for bioinformatics can cause problems:
It can change gene names to dates
https://www.techdirt.com/articles/20140727/03133828025/using-spreadsheets-bioinformatics-can-corrupt-data-changing-gene-names-into-dates.shtml

https://nsaunders.wordpress.com/2012/10/22/gene-name-errors-and-excel-lessons-not-learned/

https://www.techdirt.com/articles/20140727/03133828025/using-spreadsheets-bioinformatics-can-corrupt-data-changing-gene-names-into-dates.shtml

Again, learning how to use awk is invaluable!

Using Awk to join two files based on several columns

I was reading a thread on stackoverflow and found that this post was very interesting. I will go through the problem and the awk solution. Again, Awk is awesome!

Note:

It looks a bit messy when I copied the ipython notebook directly into Blogger. Ideally, one can write the blog using ipython notebook directly, but it does not work well with Blogger. see here http://blog.fperez.org/2012/09/blogging-with-ipython-notebook.html

I do not want to use other blogging platforms for now. So, just bear it and you can go to here to see the clean ipython notebook. Github now randers ipython notebook.

I created some dummy files.
file_a is a tab-delimited bed file with 6 colums:

In [1]:

cat file_a.bed

chr1 123 aa b c d
chr1 234 a b c d
chr1 345 aa b c d
chr1 456 a b c d

file_b is the file that contain additional infomation, which we want to add to file_a:

In [2]:

cat file_b.bed

xxxx abcd chr1 123 aa c d e
yyyy defg chr1 345 aa e f g

we want to annotate file_a based on the fact that columns 3,4,5 in file_b are the same as columns 1,2,3 in file_a.
To do this, we are going to use Awk associated array. see a link

Let me execute the awk one-liner first and then explain what's going on here:

In [7]:

awk 'NR==FNR{a[$3,$4,$5]=$1OFS$2;next}{$6=a[$1,$2,$3];print}' OFS='\t' \
file_b.bed file_a.bed

chr1 123 aa b c xxxx abcd
chr1 234 a b c 
chr1 345 aa b c yyyy defg
chr1 456 a b c 

we annotated file_a using file_b. Aka, we added first two columns from file_b to file_a.
There are several things happening here:
we see built-in variables in awk: NR and FNR. NR is the line number of the
current processing line.
when awk read in multiple files, awk NR variable will give the total number of records
relative to all the input file. Awk FNR will give you number of records for each inpu
file. see a link
here
for all the built-in variables in awk.
Let's deomonstrate the difference between NR and FNR:

In [5]:

awk '{print FILENAME, NR}' file_a.bed file_b.bed

file_a.bed 1
file_a.bed 2
file_a.bed 3
file_a.bed 4
file_b.bed 5
file_b.bed 6

FILENAME is another built-in variable for the input file name of awk.
There are 4 lines in file_a and 2 lines in file_b, and NR increments for the total lines.
compare with FNR:

In [6]:

awk '{print FILENAME, FNR}' file_a.bed file_b.bed

file_a.bed 1
file_a.bed 2
file_a.bed 3
file_a.bed 4
file_b.bed 1
file_b.bed 2

Now, awk prints out the line numbers in respect to each file.

From the awk code, we are reading file_b first. NR==FNR means when NR equals to FNR
(this is true only for file_b) do the following: {a[$3,$4,$5]=$1OFS$2;next}.
We created an associated array named a using columns 3,4,5 in file_b as keys and
the columns 1 and 2: $1"\t"$2 as values. we set OFS="\t" in the end of the command.
next means to proceed for the next line, rather than execute the following { } code block.

when awk reads in the second file (file_a.bed), NR==FNR is not true, awk
exectues the second { } code block: {$6=a[$1,$2,$3];print}
we look up the associated array a we created from file_b.bed using
the first three columns in file_a.bed as keys, and assign column 6 to
the looked-up values and print it out the whole line.

Conclusion: Awk is very powerful in text wrangling. Once get used to the
syntax, you can do fairly complicated formatting in an awk one-liner.
I strongly recommand you to learn it.

IBash Notebook for reproducible research

I use command line a lot. It is awesome for data processing, text data formatting and even exploratory data analysis.

Last week, one of my colleagues complained that she forgot how she got the data and processed the data. With a future "ME" in mind, one needs to do extensive documentations of where, when and how you download and process the data, and document the versions of the tools used in the analysis.

Although there are many ways to make tasks on command line reproducible such as using  Drake and GNU make, it is still not as straightforward as using Ipython Notebook for python and R markdown files for R, respectively.

Luckily, I got to know from Jeroen Janssens, who wrote the "Data science at command line" Book, that there is a bash_kernal for Ipython notebook, and I gave it a try.
see a screenshot of the notebook:

see the whole notebook here: http://nbviewer.ipython.org/gist/crazyhottommy/71e0dcb6d678c137733c#

Essentially, I copied the .ipynb file (it is a JSON file) and pasted it to a gist, and insert the gist link to the nbviewer website http://nbviewer.ipython.org/

With IBash Notebook, one can document the linux commands in a real-time manner and make his research more reproducible!

get all the promoter sequences of human hg19 genome

One of my former friends (a biologist who does not know much about computer) wants me to help her to get all the promoter sequences from the human genome. It is a very simple task and I think all the biologists should know how to do it. There are many ways to do it, but I will show you how to do it using bioconductor.

Using less I see all the sequences are there, to make sure the sequences are right, one can manually inspect the UCSC genome browser for several sequences.

I really do not want to dig in (google is your friend) to find a way to write the name using the SYMBOL rather than the ENTREZID....
You can convert the gene ids by many ways too.
I have two posts for that http://crazyhottommy.blogspot.com/2014/09/converting-gene-ids-using-bioconductor.html
and http://crazyhottommy.blogspot.com/2014/09/mapping-gene-ids-with-mygene.html
In addition, I prefer to prepare a bed file for all the promoters using bedtools slop (RefSeq table from UCSC, or from a GENCODE GTF file). Then, use bedtools to extract DNA sequences using bedtools getfasta. To me, it is more flexible on the command lines.
see my previous post here http://crazyhottommy.blogspot.com/2015/02/fetch-genomic-sequences-from-coordinates.html

simulation of distribution of means draw from exponential distribution

"central limit theorem (CLT) states that, given certain conditions, the arithmetic mean of a sufficiently large number of iterates of independent random variables, each with a well-defined expected value and well-defined variance, will be approximately normally distributed, regardless of the underlying distribution." [1]

I am going to draw 40 numbers from exponential distribution [2] (you can do it for any distribution) for 1000 times and examine the distribution of the means. In R, you can do it by rexp(40, lambda),  where lambda is the rate parameter. The mean of exponential distribution is 1/lambda and the standard deviation is also 1/lambda. Set lambda = 0.2 for all of the simulations for this specific test.

see a gist below.

[1] http://en.wikipedia.org/wiki/Central_limit_theorem

[2] http://en.wikipedia.org/wiki/Exponential_distribution

My first Software Carpentry workshop as an instructor

I just came back from the software-carpentry workshop held on April 2 and 3 at the University of Miami. The workshop link is here http://xuf12.github.io/2015-04-02-umiami/
and the Etherpad:https://etherpad.mozilla.org/2015-04-02-umiami

I taught Shell and R on the morning sessions and I enjoyed it very much. Most of the participants (biologists) are beginners for Shell and R, so I taught the very basics and keep my pace slow. It looks like they also enjoyed learning with me. It was my first time teaching to such a group of (around 15) people and I got myself refreshed on the basics.

The organizer Sawsan is a very nice lady and we enjoyed talking about teaching and (data) science. I met the other three instructors: Xu Fei, Matthew and Ashiwin. They are all wonderful instructors and I learned a lot from them. I also leaned a lot from the Git and SQL sessions. I think I should at least finish the online instructor training by Greg Wilson. I was in a job transition and missed the last session.  I wish I could attend more such workshops to refine my teaching skills.

Despite it was my first time teaching, I got encouraged by the students comments:

I still have a lot to improve though. I will need to speak better english (As a non-native speaker, it sometimes can be quite challenging). I am happy that I made most of the students understood :)

benchmarking for shuf vs fast_sample

Sometimes we want to randomly select a proportion of lines from a txt file. The easiest way is to use the Unix command shuf. On a mac machine, you can install it by home brew.

brew install coreutils

But you need to invoke it as gshuf. https://www.topbug.net/blog/2013/04/14/install-and-use-gnu-command-line-tools-in-mac-os-x/
I also came across a tool called fast_sample that can do the same thing
https://github.com/earino/fast_sample
I did some benchmarking for them.

The take home message for the benchmarking is that Unix tools sometimes are better than tools you write in terms of speed and memory efficiency.

understanding p value, multiple comparisons,FDR and q value

UPDATE 09/23/1016.
please also read this Nature Biotech primer

How does multiple testing correction work?

I am writing this post for my own later references. Deep understanding of p-value, FDR and q-value is not trivial, and many biologists are misusing and/or misinterpretering them.

For biologists' sake, I will use an example of gene expression. Suppose we have two groups of cells: control and treatment (can be anything like chemical treatment, radiation treatment etc..). We are looking if Gene A is differentially expressed or not under treatment. Each group we have 12 replicates.

What we usually do is take the average of 12 replicates of each group and do a t-test to compare if the difference is significant or not (assume normal distribution). We then get a p-value, say p = 0.035. We know it is smaller than 0.05 (a threshold we set), and we conclude that after treatment, expression of Gene A is significantly changed. However, what does it mean by saying a p value of 0.035?

Everything starts with a null hypothesis:
H0 : There are no difference of gene expression for Gene A after treatment.

and an alternative hypothesis:
H1:  After treatment, expression of Gene A changes.

The definition of every P value begins by assuming a null hypothesis is True[1]. With a p-value of 0.035, it means that under the Null, the probability that we see the difference of gene expression after treatment is 0.035, which is very low. If we choose a significant level of alpha=0.05, we then reject the Null hypothesis and accept the alternative hypothesis. So, if you can not state what the null hypothesis is, you can not understand the P value[1].

For a typical genomic study, there are thousands of genes we want to compare. How do we report the gene list containing the genes that are differentially expressed? We can perform a-test for each single gene and if the p-value is smaller than 0.05, we report it. However, it will give us a lot of false positives because we did not consider multiple tests. see the gist below.

Even if we randomly generated the data, you still see some pvalues are smaller than 0.05 (lower figure)!! We randomly generated data, there should be no genes that differentially expressed. However, we see a flat line across different p values.
p values are random variables. Mathmetically, one can demonstrate that under the null hypothesis (and some assumptions are met, in this case,  the test statistic T follows standard normal distritubtion), p-values follow a uniform (0,1) distribution, which means that P(p < p1) = p1. This means that the probability see a p value smaller than p1 is equal to p1. That being said, with a 100 t-tests, under the null (no difference between control and treatment), we will see 1 test with a p value smaller than 0.01. And we will see 2 tests with a p value smaller than 0.02 etc...
This explains why we see some p-values are smaller than 0.05 in our randomly generated numbers.

How do we control the false positives for multiple comparisons?
One way is to use the Bonferroni correction to correct the familywise error rate (FWER):
define a particular comparison as statistically significant only when the P value is less than alpha(often 0.05) divided by the number of comparisons (p < alpha/m) [2].  Say we computed 100 t-tests, and got 100 p values, we only consider the genes with a p value smaller than 0.05/100 as significant. This approach is very conservative and is used in Genome-wide association studies (GWAS). Since we often compare millions of genetic variations between (tens of thousands) cases and controls, this threshold will be very small! [3]

Alternatively, we can use False Discovery Rate (FDR) to report the gene list.
FDR = #false positives/# called significant.
This approach does not use the term statistically significant but instead use the term discovery.
Let's control FDR for a gene list with FDR = 0.05.
It means that of all the discoveries, 5% of them is expected to be false positives.

Benjamini & Hochberg (BH method) in 1995 proposed a way to control FDR:
Let k be the largest i such that p(i) <= (i/m) * alpha, (m is the number of comparisons)
then reject H(i) for i =1, 2, ...k
This process controls the FDR at level alpha. The method sets a different threshold p value for each comparison.  Say we computed 100 t-tests, and got 100 p values, and we want to control the FDR =0.05. We then rank the p values from small to big.
if p(1) <= 1/100 * 0.05, we then reject null hypothesis and accept the alternative.
if p(2) < = 2/100 * 0.05, we then reject the null and accept the alternative..
.....

Let's zoom in to look at the first 15 smallest p values

we can see that the 14th p value is bigger than its own threshold ,which is computed by (0.05/m) * 14 = 7.960878e-05
we will use p.adjust function and the method "fdr" or "BH" to correct the p value, what the p.adjust function does is to recalculate the p-values.
p(i)<= (i/m) * alpha
p(i) * m/i <= alpha
we can then only accept the returned the p values if p.adjust(pvals) <= alpha
sum( p.adjust(pvals, method="fdr") < 0.05 )
# 13, the same as we saw from the figure

Another method by Storey in 2002 is the direct approach to FDR:
Let K be the largest i such that pi_0 * p(i) < (i/m) * alpha
then reject H(i) for i =1,2,...k
pi_0 is the estimate of the proportion of null hypothesis in the gene list is true, range from 0 to 1.
so when pi_0 is 1, then we have the Benjamini & Hochberg correction.
This method is less conservative than the BH method.
Use the qvalue function in the bioconductor package "qvalue"
sum( qvalue(pvals)$qvalues < 0.05)
# 17, less conservative than the BH method

Note that FDR is a property of a list of genes.
q value is defined for a specific gene:

"But if you do want to assign a number to each gene, a simple thing you can do, is you can go gene by gene, and decide what would be the smallest FDR I would consider, that would include this gene in the list. And once you do that, then you have defined a q-value. And this is something that is very often reported in the list of genes"[4]

"To define the q-value we order features we tested by p-value then compute the FDRs for a list with the most significant, the two most significant, the three most significant, etc... The FDR of the list with the, say, m most significant tests is defined as the q-value of the m-th most significant feature. In other words, the q-value of a feature, is the FDR of the biggest list that includes that gene" [5]

updated on 03/16/2015
See a paper in PlosOne

The Extent and Consequences of P-Hacking in Science

http://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.1002106

and a post here: http://acceleratingscience.com/proteomics/whats-true-whats-false-proteostats-and-the-fdr/?utm_source=twitterfeed&utm_medium=twitter

References:

[1] Intuitive biostatistics by Harvey Motulsky Third Edition page 127.
[2] Intuitive biostatistics by Harvey Motulsky Third Edition page 187.
[3] Intuitive biostatistics by Harvey Motulsky Third Edition page 188.
[4] HarvardX: PH525.3x Advanced Statistics for the Life Sciences, week1, video lecture for FDR.
[5] HarvardX: PH525.3x Advanced Statistics for the Life Sciences, week1, quiz for FDR.
https://courses.edx.org/courses/HarvardX/PH525.3x/1T2015/courseware/92dbe89dc80c4fc084cad3f00b0381f7/