Misconcept of Epigenetics

There is a great buzz on twitter discussing what is epigenetics and what is not.

Epigenetics as a term is meaningless unless it refers to heritability. Otherwise it’s all transcriptional regulation @chrisw_thinkdox @nih

— Einstein Epigenomics (@EpgntxEinstein) April 8, 2016

I found this PNAS paper is worth to read: Epigenetis: Core misconcept

Some quotes from the paper:
Development is a process involves no changes of underlying DNA. It is driven primary by regulatory proteins called transcription factors (TFs) which binds to specific regulatory DNA sequence to activate or repress gene expression. This process starts from fertilized eggs where TFs are mainly contributed by mother. 

In the first cell division these transcription factors are distributed asymmetrically to the daughter cells. Then, and in subsequent cell generations, new patterns of gene expression arise in two ways: as a matter of course (some transcription factors activating expression of genes encoding other regulatory proteins, etc.), and in response to signals sent by other cells. Signaling used in development affects expression of genes encoding yet more transcription factors—it does so by changing activities of transcription factors, already present, which target those genes. Retinoic acid (a small molecule) and growth hormone (a protein) are examples of such signaling molecules.

The Core Concepts Misconception

Curiously, the picture I have just sketched is absent from the Core Concepts article. Rather, it is said, chemical modifications to DNA (e.g., methylation) and to histones—the components of nucleosomes around which DNA is wrapped in higher organisms—drive gene regulation. This obviously cannot be true because the enzymes that impose such modifications lack the essential specificity: All nucleosomes, for example, “look alike,” and so these enzymes would have no way, on their own, of specifying which genes to regulate under any given set of conditions.

Histone modifications are called “epigenetic” in the Core Concepts article, a word that for years has implied memory (see Epigenetic). This is odd: It is true that some of these modifications are involved in the process of transcription per se—facilitating removal and replacement of nucleosomes as the gene is transcribed, for example (9, 10). And some are needed for certain forms of repression (11). But all attempts to show that such modifications are “copied along with the DNA,” as the article states, have, to my knowledge, failed.

Epigenetics involves two aspects: specificity and memory (heritability). In a word, it is not histone modification drives gene expression and histone modifications != epigenetics.

For more readings:

Another paper on Genetics: what do you mean: "epigenetic?"

Paper: Epigenetics, cellular memory and gene regulation

Two consecutive posts on Epigenetics: what is it and what isn't it?

A post from professor Michael Eisen from UC berkely: The imprinter of all maladies

The New Yorker screws up big time with science: researchers criticize the Mukerjee piece on epigenetics.

From the same author above:  TED gets epigenetics wrong, but the juggernaut rolls on

A post from the Center for Epigenomics ant the Albert Einstein College of Medicine: The bright future of epigenetics

Growing pains for field of epigenetics as some call for overhaul.

How to deal with paired samples for snakemake

I am currently experimenting snakemake for pipelines. We will have paired samples when processing data. e.g. ChIP-seq (IP vs Input), DNA-seq (Normal vs Tumor). How to specify them in snakemake?
The trick is to use a python dictionary to pair samples and a function to define Input files.

I will show an toy example with ChIP-seq processing.

My first ever minimal working ChIP-seq pipeline using snakemake

Why using snakemake

Snakemake is a python3 based pipeline building tool (a python variant of GNU make) specialized for bioinformatics. I put my notes managing different versions of python here. You can write any python codes inside the Snakefile. Using snakemake is to simplify the tedious pre-processing work for large genomic data sets and for the sake of reproducibility. There are many other tools you can find here for this purpose.

Key features of snakemake

• Snakemake automatically creates missing directories.
• wildcards and Input function

To access wildcards in a shell command: {wildcards.sample}

{wildcards} is greedy (.+): {sample}.fastq could be matching sampleA.fastq if there is no sub-folder anymore, but evenwhateverfolder/sampleA.fastq can be matched as well.

One needs to think snakemake in a bottom-up way: snakemake will first look for the output files, and substitue the {wildcards} with the file names, and look for which rule can be used to creat the output, and then look for input files that are defined by the {wildcards}.

Read the following

flexible bioinformatics pipelines with snakemake
Build bioinformatics pipelines with Snakemake
snakemake ChIP-seq pipeline example
submit all the jobs immediately
snakemake-parallel-bwa
RNA-seq snakemake example
functions as inputs and derived parameters
snakemake FAQ
snakemake tutorial from the developer

examples

https://github.com/slowkow/snakefiles/blob/master/bsub.py
https://github.com/broadinstitute/viral-ngs/tree/master/pipes

A working snakemake pipeline for ChIP-seq

The folder structure is like this:

├── README.md
├── Snakemake
├── config.yaml
└── rawfastqs
    ├── sampleA
    │   ├── sampleA_L001.fastq.gz
    │   ├── sampleA_L002.fastq.gz
    │   └── sampleA_L003.fastq.gz
    ├── sampleB
    │   ├── sampleB_L001.fastq.gz
    │   ├── sampleB_L002.fastq.gz
    │   └── sampleB_L003.fastq.gz
    ├── sampleG1
    │   ├── sampleG1_L001.fastq.gz
    │   ├── sampleG1_L002.fastq.gz
    │   └── sampleG1_L003.fastq.gz
    └── sampleG2
        ├── sampleG2_L001.fastq.gz
        ├── sampleG2_L002.fastq.gz
        └── sampleG2_L003.fastq.gz

There is a folder named rawfastqs containing all the raw fastqs. each sample subfolder contains multiple fastq files from different lanes.

In this example, I have two control (Input) samples and two corresponding case(IP) samples.

CONTROLS = ["sampleG1","sampleG2"]
CASES = ["sampleA", "sampleB"]

putting them in a list inside the Snakefile. If there are many more samples, need to generate it with pythonprogrammatically.

## dry run
snakemake -np

## work flow diagram
snakemake --forceall --dag | dot -Tpng | display

To Do:

• Make the pipeline more flexiable. e.g. specify the folder name containing raw fastqs, now it is hard coded.
• write a wrapper script for submitting jobs in moab. Figuring out dependencies and --immediate-submit

Downsampling bam files to a certain number of reads

When compare different ChIP-seq data sets or analyze a set of ChIP-seq data sets together (e.g. ChromHMM analysis), it is desirable to subsample the deeply sequenced ones to a certain number of reads (say 15million or 30 million).

In paper Integrative analysis of 111 reference human epigenomes:

To avoid artificial differences in signal strength due to differences in sequencing depth, all consolidated histone mark data sets (except the additional histone marks the seven deeply profiled epigenomes, Fig. 2j) were uniformly subsampled to a maximum depth of 30 million reads (the median read depth over all consolidated samples). For the seven deeply profiled reference epigenomes (Fig. 2j), histone mark data sets were subsampled to a maximum of 45 million reads (median depth). The consolidated DNase-seq data sets were subsampled to a maximum depth of 50 million reads (median depth). These uniformly subsampled data sets were then used for all further processing steps (peak calling, signal coverage tracks, chromatin states).

After reading several posts here and here. It seems samtools and sambamba are the tools to use, but they both output a proportion number of reads.

time samtools view -s 3.6 -b my.bam -o subsample.bam
real    6m9.141s
user    5m59.842s
sys 0m8.912s

time sambamba view -f bam -t 10 --subsampling-seed=3 -s 0.6 my.bam -o subsample.bam
real    1m34.937s
user    11m55.222s
sys 0m29.872s

-s 3.6 set seed of 3 and 60% of the reads by samtools. Using multiple cpu with sambamba is much faster and an index file is generated on the fly.

If one wants to get say 15 million reads, one needs to do samtools flag stat or samtools idxstats to get the total number of reads, and then calculate the proportion by: 15 million/total = proportion.

samtools idxstats is much faster when the bam is sorted and indexed:

Retrieve and print stats in the index file. The output is TAB delimited with each line consisting of reference sequence name, sequence length, # mapped reads and # unmapped reads.

Total number of reads: samtools idxstats example.bam | cut -f3 | awk 'BEGIN {total=0} {total += $1} END {print total}'

Finally, feed the proportion to -s flag. One might want to remove the unmapped the reads and the duplicated reads in the bam file before downsampling. One might also need to sort the subsampled bam file again and index it.

extract all kinds of archives on the fly using extract command

If you are tired of remembering (googling) how to extract different kinds of archives like me, use this gem I found.



The extract function can deal with all kinds of archives.  I found it is very handy.
see it from the original author on github.

library size normalization for ChIP-seq

I have discussed how to use DESeq2 to do differential binding for ChIP-seq at here.

I am experimenting DiffBind to do the same thing, which internally uses EdgR, DESeq and DESeq2. The author Rory Stark is very responsive on the bioconductor support site and has answered several of my questions.

Today, I am going to keep a note here for normalizing the ChIP-seq data. If one compares ChIP-seq versus RNA-seq data, they are in the end all counts data. For RNA-seq, we usually get a read count table for the counts in the exons (union of them is for a gene); for ChIP-seq, we get a read count table for counts within the peaks. The peaks have to be identified by other tools such as MACS first. The counts data follow a (negative) binomial distribution. That's why tools such as DESeq2, which was developed for RNAseq is used for ChIP-seq.

After we get a count table, it comes to the normalization problem. If you are interested, read this paper Beyond library size: a field guide to NGS normalization. In the DiffBind package, the counts table is obtained by a function ?dba.count.

There are several ways to specify how the counts are normalized for the binding affinity matrix:

score   
which score to use in the binding affinity matrix. Note that all raw read counts are maintained for use by dba.analyze, regardless of how this is set. One of:
DBA_SCORE_READS raw read count for interval using only reads from ChIP
DBA_SCORE_READS_FOLD    raw read count for interval from ChIP divided by read count for interval from control
DBA_SCORE_READS_MINUS   raw read count for interval from ChIP minus read count for interval from control
DBA_SCORE_RPKM  RPKM for interval using only reads from ChIP
DBA_SCORE_RPKM_FOLD RPKM for interval from ChIP divided by RPKM for interval from control
DBA_SCORE_TMM_READS_FULL    TMM normalized (using edgeR), using ChIP read counts and Full Library size
DBA_SCORE_TMM_READS_EFFECTIVE   TMM normalized (using edgeR), using ChIP read counts and Effective Library size
DBA_SCORE_TMM_MINUS_FULL    TMM normalized (using edgeR), using ChIP read counts minus Control read counts and Full Library size
DBA_SCORE_TMM_MINUS_EFFECTIVE   TMM normalized (using edgeR), using ChIP read counts minus Control read counts and Effective Library size
DBA_SCORE_TMM_READS_FULL_CPM    same as DBA_SCORE_TMM_READS_FULL, but reporrted in counts-per-million.
DBA_SCORE_TMM_READS_EFFECTIVE_CPM   same as DBA_SCORE_TMM_READS_EFFECTIVE, but reporrted in counts-per-million.
DBA_SCORE_TMM_MINUS_FULL_CPM    same as DBA_SCORE_TMM_MINUS_FULL, but reporrted in counts-per-million.
DBA_SCORE_TMM_MINUS_EFFECTIVE_CPM   Tsame as DBA_SCORE_TMM_MINUS_EFFECTIVE, but reporrted in counts-per-million.

DBA_SCORE_TMM_READS_FULL vs DBA_SCORE_TMM_READS_EFFECTIVE:

Diffbind let's you to choose use full library size or effective library size for trimmed mean of M values(TMM) normalization which was proposed by Mark D Robinson for RNAseq.

Full library size is the number of reads in the bam files.
Effective library size is the number of reads mapped in the exons or within the peaks. It is the column sums for the matrix.

Note that effective library size (bFullLibrarySize =FALSE) may be more appropriate for situations when the overall signal (binding rate) is expected to be directly comparable between the samples.

If one wants to subtract the input reads, one can use DBA_SCORE_TMM_MINUS_FULL and DBA_SCORE_TMM_MINUS_EFFECTIVE

No matter what score you choose, for differential binding analysis in Diffbind, it is always the raw counts is used for the binding matrix. Diffbind (by default) subtract the input raw reads for subsequent analysis. Whether or not this is good was discussed here.

For example, if one uses DESeq2, the details are as follows:

For each contrast, a separate analysis is performed. First, a matrix of counts is constructed for the contrast, with columns for all the samples in the first group, followed by columns for all the samples in the second group. The raw read count is used for this matrix; if the bSubControl parameter is set to TRUE (as it is by default), the raw number of reads in the control sample (if available) will be subtracted. Next the library size is computed for each sample for use in subsequent normalization. By default, this is the total number of reads in peaks (the sum of each column). Alternatively, if the bFullLibrarySize parameter is set to TRUE, the total number of reads in the library (calculated from the source BAM/BED file) is used. The first step concludes with a call to DESeq2’s DESeqDataSetFromMatrix function, which returns a DESeqDataSet object. If bFullLibrarySize is set to TRUE, then sizeFactors is called with the number of reads in the BAM/BED files for each ChIP sample, divided by the minimum of these; otherwise, estimateSizeFactors is invoked. estimateDispersions is then called with the DESeqDataSet object and fitType set to local. Next the model is fitted and tested using nbinomWaldTest

estimateSizeFactors in DESeq2:

Given a matrix or data frame of count data, this function estimates the size factors as follows: Each column is divided by the geometric means of the rows. The median (or, ir requested, another location estimator) of these ratios (skipping the genes with a geometric mean of zero) is used as the size factor for this column.

The most powerful unix commands I have learned so far: find + parallel

I have been using unix since 2013, yet I learn unix tricks almost everyday. The most powerful commands I have learned so far is the find, xargs and parallel commands.

please check parallel GNU page for documentations.

parallel has changed my way to do repetitive works. Now, I use fewer and fewer for loops.

Use case 1: I have 100 folders with names starting with H3K4me3, inside each folder, I have 5 .gz files that I want to cat together. The usual way to do it:

# !/bin/bash

for dir in H3K4me3*/
do
    cd $dir && cat *H3K4me3.bed.gz > ${dir}_merged.gz 

    cd ..

done

Note that cat works well with *gz files.

The parallel way:

 ls -d H3K4me3* | parallel 'find {} -name "*H3K4me3*bed.gz" | xargs cat > {}_H3K4me3.bed.gz'

Using parallel, I can take full advantage of the multi-core nodes on the computing cluster, so it is much faster.

Use case 2: I have 100 folders (50 folder names start with H3K4me3, 50 start with H3K4me), each folder has multiple levels of sub-folders. I want to delete some bam files in 50 of them with name starting with H3K4me3, but I do not know which sub-folder the bam files may exist.

I do not really know a way to do it without using find. My solution would be:

ls -d H3K4me3* | parallel 'find {} -name "*bam"' | parallel rm {}

piping to two parallel is the magic of this solution. Unix commands are elegant and efficient!!

Edit on 04/04/2016:

With greater power comes greater responsibility. When you have too many files to process,
it is good to restrict parallel to only use certain number of CPUs with -j and not use swap-memory --noswap.

Breakpoints clustering for structural variants

Ming Tang

March 2, 2016

When structural variants are represented in bedpe format, each structural variant is represented as two linked breakpoints. Each breakpoint is represented as a genomic interval (GRanges). It is very similar to representations of Genomic interactions (data from Hi-C, ChIA-PET):

e.g.

chr1 100 200 chr2 300 400

one of the common task is to merge the breakpoints when they are overlapping.

———|———|———————–|———-|—–

—–|———|————————-|—————|—

merged to

—–|————-|———————|—————|—

It is not a trivial problem as stated by Aaron Qunlan in this post

It is the so called “breakpoint clustering” problem. I am going to use InteractionSet bioconductor package to solve this problem.

Please also check clusterPairs function in the diffHic package.

see answer from the author of InteractionSet. Note that InteractionSet is still under development.  Please refer to the tutorial for more usages. Install it using devtools:install_github("LTLA/InteractionSet").

Make some dummy GRanges.

library(InteractionSet)
all.regions <- GRanges(rep("chrA",8), 
    IRanges(c(1,6,2,9,5,2,15,20), c(3,10,4,12,7,4,18,23)))
index.1 <- c(1,3,5,7)
index.2 <- c(2,4,6,8) 

Using mode=strict when constructing the GInteraction object or using the swapAnchors method is to ensure that the first anchor index is always less than the second anchor index for each interaction. This eliminates redundant permutations of anchor regions and ensures that an interaction between regions #1 and #2 is treated the same as an interaction between regions #2 and #1. Obviously, this assumes that redundant permutations are uninteresting.

gi <- GInteractions(index.1, index.2, all.regions, mode ="strict")

gi

## StrictGInteractions object with 4 interactions and 0 metadata columns:
##       seqnames1   ranges1     seqnames2   ranges2
##           <Rle> <IRanges>         <Rle> <IRanges>
##   [1]      chrA  [ 1,  3] ---      chrA  [ 6, 10]
##   [2]      chrA  [ 2,  4] ---      chrA  [ 9, 12]
##   [3]      chrA  [ 2,  4] ---      chrA  [ 5,  7]
##   [4]      chrA  [15, 18] ---      chrA  [20, 23]
##   -------
##   regions: 8 ranges and 0 metadata columns
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

The first three pairs of interactions can be merged. Note that two anchors of each pair are overlapping. First, use findOverlaps for GInteraction object to find out which pairs with two anchors are overlapping. (2D overlapping)

out<- findOverlaps(gi)

out

## Hits object with 8 hits and 0 metadata columns:
##       queryHits subjectHits
##       <integer>   <integer>
##   [1]         1           1
##   [2]         1           2
##   [3]         1           3
##   [4]         2           1
##   [5]         2           2
##   [6]         3           1
##   [7]         3           3
##   [8]         4           4
##   -------
##   queryLength: 4 / subjectLength: 4

This will identify all pairs of interactions in gi that have two-dimensional overlaps with each other, i.e., both anchor regions in one interaction overlap with corresponding anchor regions in the other interaction. Using graph algorithm to cluster the anchors (I myself is by no means an algorithm person :)). We then do:

library(RBGL)

## need to trick the function by specifying edgemode to be directional
g <- ftM2graphNEL(as.matrix(out), W=NULL, V=NULL, edgemode="directed")

## if you want to see what's going on
library(Rgraphviz)
plot(g)

## change it back to undirected
edgemode(g) <- "undirected"

plot(g)

connections <- connectedComp(g)

To identify groups of interactions that overlap at least one other interaction in the same group (i.e., “single-linkage clusters” of interactions that have overlapping areas in the two-dimensional interaction space). The clusters can be explicitly identified using:

cluster <- integer(length(gi))
for (i in seq_along(connections)) {
    cluster[as.integer(connections[[i]])] <- i
}

We can then identify the bounding box (merge the overlapping interactions) for each cluster of interactions with:

boundingBox(gi, cluster)

## GInteractions object with 2 interactions and 0 metadata columns:
##     seqnames1   ranges1     seqnames2   ranges2
##         <Rle> <IRanges>         <Rle> <IRanges>
##   1      chrA  [ 1,  4] ---      chrA  [ 5, 12]
##   2      chrA  [15, 18] ---      chrA  [20, 23]
##   -------
##   regions: 4 ranges and 0 metadata columns
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

overlap with one anchor and get the other anchor

I have of.interest overlaps with the second anchor, and I want to return the first anchor. how can I do it?

e.g.

    first (promoters)                      second

——|————-|———————-|————|——– gi

                             |------------|        of.interest (enhancers)

of.interest <- GRanges(c("chrA","chrA"), IRanges(start=c(2, 5),end=c(8,8)))

# get overlaps between of.interest and first anchor
hits1 <- findOverlaps(gi, of.interest, use.region="first")
anchors(gi[queryHits(hits1),], type="first") # overlapping

## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrA    [1, 3]      *
##   [2]     chrA    [2, 4]      *
##   [3]     chrA    [2, 4]      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

anchors(gi[queryHits(hits1),], type="second") # "other"

## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrA   [6, 10]      *
##   [2]     chrA   [9, 12]      *
##   [3]     chrA   [5,  7]      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

of.interest[subjectHits(hits1)]

## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrA    [2, 8]      *
##   [2]     chrA    [2, 8]      *
##   [3]     chrA    [2, 8]      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

# get overlaps # between of.interest and second anchor
hits2 <- findOverlaps(gi, of.interest, use.region="second")
anchors(gi[queryHits(hits2),], type="second") # overlapping

## GRanges object with 4 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrA   [6, 10]      *
##   [2]     chrA   [6, 10]      *
##   [3]     chrA   [5,  7]      *
##   [4]     chrA   [5,  7]      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

anchors(gi[queryHits(hits2),], type="first") # "other"

## GRanges object with 4 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrA    [1, 3]      *
##   [2]     chrA    [1, 3]      *
##   [3]     chrA    [2, 4]      *
##   [4]     chrA    [2, 4]      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

of.interest[subjectHits(hits2)]

## GRanges object with 4 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chrA    [2, 8]      *
##   [2]     chrA    [5, 8]      *
##   [3]     chrA    [2, 8]      *
##   [4]     chrA    [5, 8]      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

If you’re willing to sacrifice some information, you might consider using the linearize method:

gi$index <- seq_along(gi)
linearize(gi, of.interest[1])

## GRanges object with 3 ranges and 1 metadata column:
##       seqnames    ranges strand |     index
##          <Rle> <IRanges>  <Rle> | <integer>
##   [1]     chrA   [1, 10]      * |         1
##   [2]     chrA   [9, 12]      * |         2
##   [3]     chrA   [2,  7]      * |         3
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

linearize(gi, of.interest[2])

## GRanges object with 2 ranges and 1 metadata column:
##       seqnames    ranges strand |     index
##          <Rle> <IRanges>  <Rle> | <integer>
##   [1]     chrA    [1, 3]      * |         1
##   [2]     chrA    [2, 4]      * |         3
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

This will return a GRanges containing the “other” anchor region for all interactions that overlap each entry in of.interest (you can figure out what those interactions were by looking at index in the returned object). linearize also provides options for handling “internal” interactions where both anchor regions overlap the specified entry of of.interest; by default, the union of the two anchor regions is returned for such interactions, but they can be removed entirely by setting internal=FALSE.

Alternatively, if you already have two sets of regions(promoters, enhancers) you want to know which enhancer links to which promoter

all.enhancers<- GRanges(c("chrA","chrA"), IRanges(c(2,17), c(5,19)))
all.genes<- GRanges(c("chrA","chrA"), IRanges(c(4,6), c(7,16)))

linkOverlaps(gi, all.genes, all.enhancers)

##   query subject1 subject2
## 1     1        1        1
## 2     1        2        1
## 3     2        2        1
## 4     3        1        1
## 5     3        2        1

which two genes are linked?

linkOverlaps(gi, all.genes)

##   query subject1 subject2
## 1     2        2        1
## 2     3        1        1
## 3     3        2        1

which two enhancers are linked?

linkOverlaps(gi, all.enhancers)

##   query subject1 subject2
## 1     3        1        1