Question: ATAC-seq DE analysis
1
gravatar for collmerr
9 days ago by
collmerr10
collmerr10 wrote:

Hello all. I'm trying to do differential expression analysis using ATAC-seq data sets, but I'm having trouble getting HTseq-count or featureCounts to actually count the data. Any good protocols to accomplish this?

Thanks!

ADD COMMENTlink modified 8 days ago by trausch1.0k • written 9 days ago by collmerr10
1

HT-seq or featureCounts are not really suitable (tedious work to prepare an input format with bed regions) for ChIP/ATAC data. Use bamCoverage or multiBamSummary from deeptools. All you need is BAMs and a bed file.

ADD REPLYlink modified 8 days ago • written 8 days ago by venu5.3k

Disagreed. See my comment below. bedtools does not really support paired-end data, while featureCounts perfectly does. Preparing a featureCounts input file from BED is as easy as:

awk 'OFS="\t" {print $1"."$2+1"."$3, $1, $2+1, $3, "."}' peaks.narrowPeak > peaks.saf
ADD REPLYlink written 8 days ago by ATpoint4.4k
1
gravatar for b.nota
8 days ago by
b.nota4.0k
Netherlands
b.nota4.0k wrote:

The way I approach this, is to first use MACS2 for peak calling using all BAM files together. This will generate a bed file with coordinates of all peaks. Then I use bedtools, something like this:

bedtools coverage -b ${bamFile} -a All_peaks.bed > CoverageFiles/${bamFile%.bam}.txt

And then merge all columns with the read counts from the files for edgeR.

ADD COMMENTlink modified 8 days ago • written 8 days ago by b.nota4.0k

The issue with bedtools is one can not input multiple bam files at once. With multiBamSummary from deeptools, you can input as many BAM files as possible, it will generate one counts table for all BAMs.

ADD REPLYlink written 8 days ago by venu5.3k
1

Bedtools can count alignments from multiple BAM files (see http://bedtools.readthedocs.io/en/latest/content/tools/multicov.html)

ADD REPLYlink written 8 days ago by James Ashmore2.4k

+1. I completely forgot about multicov function. Thanks for adding.

ADD REPLYlink written 8 days ago by venu5.3k

As ATAC-seq is typically paired-end, I recomennd featureCounts from the subread package, which can handle this (bedtools cannot as far as I know), and is multithreaded.

ADD REPLYlink written 8 days ago by ATpoint4.4k

It would be helpful then to show how you do it with featureCount, see OP's remark in the question: "but I'm having trouble getting HTseq-count or featureCounts to actually count the data".

ADD REPLYlink written 8 days ago by b.nota4.0k
1

Sorry, did not notice that sentence. There we go:

featureCounts cannot take bed or narrowPeak files directly, but wants a custom format called SAF:

## Make SAF file (+1 because SAF is 1-based, BED/narrowPeak is 0-based)
awk 'OFS="\t" {print $1"."$2+1"."$3, $1, $2+1, $3, "."}' peaks.narrowPeak > peaks.saf

## Then run fc:
featureCounts -a peaks.saf -F SAF -T <int> -p -o peaks_countMatrix.txt *.bam

-T = number of threads; -p indicates paired-end data

ADD REPLYlink written 8 days ago by ATpoint4.4k
1
gravatar for Devon Ryan
8 days ago by
Devon Ryan80k
Freiburg, Germany
Devon Ryan80k wrote:

Differential accessibility, not differential expression.

We tend to use CSAW for this. The full pipeline is here, which is essentially:

  1. Fragments are filtered (alignmentSieve) so only those <=150 bases (i.e., likely from open chromatin) are present.
  2. MACS2 is run on that.
  3. CSAW is run in R and given the peaks from 2 and the unfiltered BAM files that were input into step 1.

Alternatively, you could use diffBind in the last step. As Venu pointed out, you can use multiBamSummary with a BED file to get your counts, but CSAW at least has functions that also do that.

ADD COMMENTlink written 8 days ago by Devon Ryan80k
0
gravatar for trausch
8 days ago by
trausch1.0k
Germany
trausch1.0k wrote:

Another ATAC-Seq pipeline is available here. As part of the pipeline a count matrix is generated and an Rscript template for differential peak calling is part of the pipeline (script is here).

ADD COMMENTlink written 8 days ago by trausch1.0k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1319 users visited in the last hour