Epigenomics, Session 2
Bioinformatics Resource Center - Rockefeller University
http://rockefelleruniversity.github.io/RU_course_template/
brc@rockefeller.edu
CUT&RUN/ATAC (part 2) - Alignment and Peak Calling
Set the Working directory
Before running any of the code in the practicals or slides we need to set the working directory to the folder we unarchived.
You may navigate to the unarchived RU_Course_help folder in the Rstudio menu.
Session -> Set Working Directory -> Choose Directory
or in the console.
What we will cover
Earlier we learned how to asses the quality of a FASTQ file and filter the reads based on quality scores.
We will take these filtered FASTQ files and continue with the analysis:
- Align the FASTQ reads to the genome to generate BAM alignment files.
- Identify areas of significant read build up by calling peaks.
- These steps are very similar for CUT&RUN and ATACseq, but there are some differences we will highlight.
Aligning CUT&RUN and ATACseq reads
Aligning CUT&RUN and ATACseq reads
Following assessment of read quality and any read filtering we applied, we will want to align our reads to the genome so as to identify any genomic locations showing enrichment for aligned reads above background.
Since CUT&RUN and ATACseq reads will align continuously against our reference genome we can use genomic aligners, as opposed to aligners that account for splice junctions (e.g. for RNAseq).
The resulting BAM file will contain aligned sequence reads for use in further analysis.
Creating a reference genome
First we need to retrieve the sequence information for the genome of interest in FASTA format. We can use the BSgenome libraries to retrieve the full sequence information.
For the mouse mm10 genome we load the package BSgenome.Mmusculus.UCSC.mm10.
## | BSgenome object for Mouse
## | - organism: Mus musculus
## | - provider: UCSC
## | - genome: mm10
## | - release date: Sep 2017
## | - 239 sequence(s):
## | chr1 chr2 chr3
## | chr4 chr5 chr6
## | chr7 chr8 chr9
## | chr10 chr11 chr12
## | chr13 chr14 chr15
## | ... ... ...
## | chrX_KZ289094_fix chrX_KZ289095_fix chrY_JH792832_fix
## | chrY_JH792833_fix chrY_JH792834_fix chr1_KK082441_alt
## | chr11_KZ289073_alt chr11_KZ289074_alt chr11_KZ289075_alt
## | chr11_KZ289077_alt chr11_KZ289078_alt chr11_KZ289079_alt
## | chr11_KZ289080_alt chr11_KZ289081_alt
## |
## | Tips: call 'seqnames()' on the object to get all the sequence names, call
## | 'seqinfo()' to get the full sequence info, use the '$' or '[[' operator to
## | access a given sequence, see '?BSgenome' for more information.Creating a reference genome
We will only use the major chromosomes for our analysis so we may exclude random and unplaced contigs. Here we cycle through the major chromosomes and create a DNAStringSet object from the retrieved sequences.
mainChromosomes <- paste0("chr", c(1:19, "X", "Y", "M"))
mainChrSeq <- lapply(mainChromosomes, function(x) BSgenome.Mmusculus.UCSC.mm10[[x]])
names(mainChrSeq) <- mainChromosomes
mainChrSeqSet <- DNAStringSet(mainChrSeq)
mainChrSeqSet## DNAStringSet object of length 22:
## width seq names
## [1] 195471971 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr1
## [2] 182113224 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr2
## [3] 160039680 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr3
## [4] 156508116 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr4
## [5] 151834684 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr5
## ... ... ...
## [18] 90702639 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr18
## [19] 61431566 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chr19
## [20] 171031299 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chrX
## [21] 91744698 NNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNN chrY
## [22] 16299 GTTAATGTAGCTTAATAACAA...TACGCAATAAACATTAACAA chrMCreating a reference genome
Now we have a DNAStringSet object we can use the writeXStringSet to create our FASTA file of sequences to align to.
Creating an Rsubread index
We will be aligning using the seed-and-vote algorithm from the folks behind subread. We can therefore use the Rsubread package. Before we attempt to align our FASTQ files, we will need to first build an index from our reference genome using the buildindex() function.
NOTE: Building an index is memory intensive and by default is set to 8GB. This may be too large for your laptop or Desktop.
Aligning CUT&RUN Reads
In the following slides we will focus on aligning the CUT&RUN data, but alignment of ATACseq data will be the same. In the post processing and peak calling section there are some considerations for each method that we will cover.
CUT&RUN and ATACseq data are typically paired-end sequencing, and most often comes as two files typically with *_1* and *_2* or *_R1* and *_R2* in the file name to denote which number in pair a file is.
The two matched paired-end read files will (most often) contain the same number of reads and the reads will be ordered the same in both files.
Aligning Sequence Reads
Below we read in the first 1000 reads in replicate 1 of the Week 6 CUT&RUN FASTQ file.
The read names will match across files for paired reads with the exception of a 1 or 2 in name to signify if the read was first or second in a pair.
require(ShortRead)
# first 1000 reads in each file
read1 <- readFastq("data/SOX9CNR_W6_rep1_1K_R1.fastq.gz")
read2 <- readFastq("data/SOX9CNR_W6_rep1_1K_R2.fastq.gz")
id(read1)[1:2]## BStringSet object of length 2:
## width seq
## [1] 17 SRR20110409.9 9/1
## [2] 19 SRR20110409.10 10/1## BStringSet object of length 2:
## width seq
## [1] 17 SRR20110409.9 9/2
## [2] 19 SRR20110409.10 10/2Aligning Sequence Reads
Now that we have an index and our FASTQ files, we can perform the alignment.
Our two FASTQ files that have been filtered with rfastp for QC metrics, and are available on dropbox if you have built the index and want to try the alignment on your own.
Rsubread alignment of CUT&RUN
We can align our raw sequence data in FASTQ format to the FASTA file of our mm10 genome sequence using the Rsubread package. Specifically we will be using the align function as this utilizes the subread genomic alignment algorithm.
The align() function accepts arguments for the index to align to, the FASTQ(s) to align, the name of output BAM, the mode of alignment (rna or dna), and fragment length filters.
We will go through alignment of the CUT&RUN.
Sort and index CUT&RUN reads
As before, we sort and index our files using the Rsamtools packages sortBam() and indexBam() functions respectively.
The resulting sorted and indexed BAM file is now ready for use in external programs such as IGV as well as for further downstream analysis in R.
Create a bigWig
We can also create a bigWig from our sorted, indexed BAM file to allow us to quickly review our data in IGV.
First we use the coverage() function to create an RLElist object containing our coverage scores.
RleList of length 22
$chr1
integer-Rle of length 195471971 with 1828459 runs
Lengths: 3004130 40 121 40 1365 34 6 34 113 ... 13 9 8 6 4 30 6 100490
Values : 0 1 0 1 0 1 2 1 0 ... 6 3 2 3 4 2 1 0
$chr2
integer-Rle of length 182113224 with 1855901 runs
Lengths: 3050000 31 4 71 9 30 1 3 3 ... 119 39 1 61 6 34 6 100214
Values : 0 3 2 0 1 2 1 2 1 ... 0 2 1 0 1 2 1 0
Create a bigWig
We can now export our RLElist object as a bigWig using the rtracklayer package’s export.bw() function.
Create a bigWig
We may wish to normalize our coverage to allow us to compare enrichment across samples.
First we we can retrieve and plot the number of mapped reads using the idxstatsBam() function.
mappedReads <- idxstatsBam("SOX9CNR_W6_rep1_sorted.bam")
TotalMapped <- sum(mappedReads[, "mapped"])
TotalMapped[1] 24848310Create a bigWig
To make the normalized bigwig we can use the weight parameter in the coverage() to scale our reads to the number of mapped reads multiplied by a million (reads per million).
This bigwig is available here on dropbox for download if you want to check it out.
BAM and bigWig
Post alignment processing - CUT&RUN
Peak calling software - CUT&RUN
To identify regions of SOX9 transcription factor binding we can use a peak caller.
Although many peak callers are available in R and beyond, one the most popular and widely used peak caller remains MACS.
Another peak caller, SEACR, was developed specificlly for CUT&RUN. We will use both.
Consider an appropriate input
While we don’t have one in this test data set, inputs are sometimes used in peak calling. They are theoretically less critical in CUT&RUN compared to ChIPseq due to lower background signal, but this is not always the case.
Input samples in CUT&RUN are typically made from IgG pulldowns
Allows for control of artefact regions which occur across samples.
For example: When using tumour samples for CUT&RUN, it might be important to have matched input samples. Differing conditions of same tissue may share common input.
Install software with conda and Herper
We will install a version of MACS that is not available in R. It is possible to install it on Mac and Linux using the Anaconda package repository (unfortunately there is not a Windows implementation).
Anaconda is a huge collection of version controlled packages that can be installed through the conda package management system. With conda it is easy to create and manage environments that have a variety of different packages in them.
We can interact with the Anaconda package system using the R package Herper. This was created by us here at the BRC and is part of Bioconductor.
Install software with conda and Herper
First, we will use Herper to install MACS with the
install_CondaTools function. This will automatically
install miniconda, a lightweight version of Anaconda.
We just need to tell install_CondaTools the
following:
- what tool/s you want
- the name of the environment you want to build.
- path where the miniconda currently is or should be installed
- whether you are updating a current environment (optional)
We will install MACS (version 3), in addition to samtools and bedtools, which will be used to prepare our BAM files for use by MACS and SEACR.
Install software with conda and Herper
Behind the scenes, Herper will install the most minimal version of conda (called miniconda), and then will create a new environment into which MACS, bedtools, and samtools will be installed.
When you run the function it returns and prints out where the software is installed.
$pathToConda
[1] "miniconda/bin/conda"
$environment
[1] "CnR_analysis"
$pathToEnvBin
[1] "miniconda/envs/CnR_analysis/bin"
Install software with conda and Herper
After the software are installed, we can then use them using the
with_CondaEnv function from Herper paired with the
system function that allows us to run terminal commands
from R.
samtools: invalid option -- h
Usage: samtools sort [options...] [in.bam]
Options:
-l INT Set compression level, from 0 (uncompressed) to 9 (best)
-u Output uncompressed data (equivalent to -l 0)
-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]
-M Use minimiser for clustering unaligned/unplaced reads
-R Do not use reverse strand (only compatible with -M)
-K INT Kmer size to use for minimiser [20]
-I FILE Order minimisers by their position in FILE FASTA
-w INT Window size for minimiser indexing via -I ref.fa [100]
-H Squash homopolymers when computing minimiser
-n Sort by read name (natural): cannot be used with samtools index
-N Sort by read name (ASCII): cannot be used with samtools index
-t TAG Sort by value of TAG. Uses position as secondary index (or read name if -n is set)
-o FILE Write final output to FILE rather than standard output
-T PREFIX Write temporary files to PREFIX.nnnn.bam
--no-PG
Do not add a PG line
--template-coordinate
Sort by template-coordinate
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
-O, --output-fmt FORMAT[,OPT[=VAL]]...
Specify output format (SAM, BAM, CRAM)
--output-fmt-option OPT[=VAL]
Specify a single output file format option in the form
of OPTION or OPTION=VALUE
--reference FILE
Reference sequence FASTA FILE [null]
-@, --threads INT
Number of additional threads to use [0]
--write-index
Automatically index the output files [off]
--verbosity INT
Set level of verbosityProper pairs - CUT&RUN
Now that we have our aligned paired-end data contained in BAM files, and the software installed to eventually call peaks, we can assess some initial QC metrics like alignment quality and fragment length.
We will read our newly aligned data using the GenomicAlignments package.
Here we only want reads which are properly paired so we will use the ScanBamParam() and scanBamFlag() functions to control what will be read into R.
We set the scanBamFlag() function parameters isProperPair to TRUE so as to only read in reads paired in alignment.
Proper pairs - CUT&RUN
We can now use these flags with the ScanBamParam() function to read in only properly paired reads.
We additionally specify information to be read into R using the what parameter. Importantly we specify the insert size information - isize. To reduce memory footprint we read only information from chromosome 18 by specifying a GRanges object to which parameter.
## class: ScanBamParam
## bamFlag (NA unless specified): isProperPair=TRUE
## bamSimpleCigar: FALSE
## bamReverseComplement: FALSE
## bamTag:
## bamTagFilter:
## bamWhich: 0 ranges
## bamWhat: qname, mapq, isize, flag
## bamMapqFilter: NAGAlignmentPairs - CUT&RUN
Now we have set up the ScanBamParam object, we can use the readGAlignmentPairs() function to read in our paired-end CUT&RUN data using the readGAlignments() function.
In the data folder is a BAM file contianing the reads from only chromosome 18 that we will read into R.
The resulting object is a GAlignmentPairs object that contains information on our paired reads
sortedBAM <- "data/SOX9CNR_W6_rep1_chr18_sorted.bam"
cnrReads <- readGAlignmentPairs(sortedBAM, param = myParam)
cnrReads[1:2, ]## GAlignmentPairs object with 2 pairs, strandMode=1, and 0 metadata columns:
## seqnames strand : ranges -- ranges
## <Rle> <Rle> : <IRanges> -- <IRanges>
## [1] chr18 + : 3001178-3001217 -- 3001186-3001225
## [2] chr18 + : 3002517-3002556 -- 3002544-3002583
## -------
## seqinfo: 22 sequences from an unspecified genomeGAlignmentPairs - CUT&RUN
We access the GAlignments objects using the first() and second() accessor functions to gain information on the first or second read respectively.
read1 <- GenomicAlignments::first(cnrReads)
read2 <- GenomicAlignments::second(cnrReads)
read2[1:2, ]## GAlignments object with 2 alignments and 4 metadata columns:
## seqnames strand cigar qwidth start end width
## <Rle> <Rle> <character> <integer> <integer> <integer> <integer>
## [1] chr18 - 40M 40 3001186 3001225 40
## [2] chr18 - 40M 40 3002544 3002583 40
## njunc | qname mapq isize flag
## <integer> | <character> <integer> <integer> <integer>
## [1] 0 | SRR20110409.2378427 0 -48 147
## [2] 0 | SRR20110409.10295909 10 -67 147
## -------
## seqinfo: 22 sequences from an unspecified genomeRetrieve MapQ scores - CUT&RUN
One of the first we can do is obtain the distribution of MapQ scores for our read1 and read2. We can access this for each read using mcols() function to access the mapq slot of the GAlignments object for each read.
For Rsubread, a MAPQ of 0 means that a read maps equally well to multiple locations, we want to assess if this is common and potentially remove these reads.
## [1] 0 40 40 10 10MapQ score frequencies - CUT&RUN
We can then use the table() function to summarize the frequency of scores for each read in pair.
## read1MapQ
## 0 6 8 10 13 20 40
## 3123 32 307 1848 9611 31077 357317## read2MapQ
## 0 6 8 10 13 20 40
## 3123 50 396 2232 11486 31529 354499Plot MapQ scores - CUT&RUN
Finally we can plot the distributions of MapQ for each read in pair using ggplot2.
library(ggplot2)
toPlot <- data.frame(MapQ = c(names(read1MapQFreqs), names(read2MapQFreqs)), Frequency = c(read1MapQFreqs,
read2MapQFreqs), Read = c(rep("Read1", length(read1MapQFreqs)), rep("Read2",
length(read2MapQFreqs))))
toPlot$MapQ <- factor(toPlot$MapQ, levels = unique(sort(as.numeric(toPlot$MapQ))))
ggplot(toPlot, aes(x = MapQ, y = Frequency, fill = MapQ)) + geom_bar(stat = "identity") +
facet_grid(~Read)
Retrieving insert sizes - CUT&RUN
Now we have reads in the paired aligned data into R, we can retrieve the insert sizes from the mcols() attached to GAlignments objects of each read pair.
CUT&RUN should represent a mix of fragment lengths depending on the target. A TF will generally only be in nucleosome-free regions and will have small fragments, while a histone mark might have longer fragments.
Since properly paired reads will have the same insert size length we extract insert sizes from read1.
## [1] 48 67 67 67 67 67Plotting insert sizes - CUT&RUN
We can use the table() function to retrieve a vector of the occurrence of each fragment length.
## insertSizes
## 30 31 32 33 34
## 1 574 746 1067 1386Plotting insert sizes - CUT&RUN
Now we can use the newly acquired insert lengths for chromosome 18 to plot the distribution of all fragment lengths.
library(ggplot2)
toPlot <- data.frame(InsertSize = as.numeric(names(fragLenSizes)), Count = as.numeric(fragLenSizes))
fragLenPlot <- ggplot(toPlot, aes(x = InsertSize, y = Count)) + geom_line()
fragLenPlot + theme_bw()
Plotting insert sizes - CUT&RUN
We can now annotate our nucleosome free (< 120bp) regions and see a clear bi-modal distribution of reads separated by this distinction.
Filter reads - CUT&RUN
We will filter out fragments where either read has a MAPQ of 0 and reads shorter than 1000bp. MAPQ of zero means a read mapped equally well to at least two places in the genome.
Note: For TF CUT&RUN where we expect binding in nucleosome free regions, sometimes it is common to only take reads shorter than 120bp. In the exercises you have a chance to compare both strategies.
Remove blacklist - CUT&RUN
ChIPseq and CUT&RUN will often show the presence of common artifacts, such as ultra-high signal regions. Such regions can confound peak calling, fragment length estimation and QC metrics. Anshul Kundaje created the DAC blacklist as a reference to help deal with these regions.
We will use this reference to remove these reads from the BAM file.
Remove blacklist regions - CUT&RUN
You can find a Blacklist for most genomes at the Boyle lab Github page. The most recent version (V2) for mm10 is available in the data folder of this course.
We read it into R as a GRanges object using the import
function from rtracklayer.
library(rtracklayer)
blacklist <- "data/mm10-blacklist.v2.bed"
bl_regions <- rtracklayer::import(blacklist)
bl_regions## GRanges object with 3435 ranges and 1 metadata column:
## seqnames ranges strand | name
## <Rle> <IRanges> <Rle> | <character>
## [1] chr10 1-3135400 * | High Signal Region
## [2] chr10 3218901-3276600 * | Low Mappability
## [3] chr10 3576901-3627700 * | Low Mappability
## [4] chr10 4191101-4197600 * | Low Mappability
## [5] chr10 4613501-4615400 * | High Signal Region
## ... ... ... ... . ...
## [3431] chrY 6530201-6663200 * | High Signal Region
## [3432] chrY 6760201-6835800 * | High Signal Region
## [3433] chrY 6984101-8985400 * | High Signal Region
## [3434] chrY 10638501-41003800 * | High Signal Region
## [3435] chrY 41159201-91744600 * | High Signal Region
## -------
## seqinfo: 21 sequences from an unspecified genome; no seqlengthsRemove blacklist regions - CUT&RUN
The granges function finds the range spanning the whole
fragment (5’ read 1 to 3’ read 2) when given a GAlignmentPairs
object.
We then identify the fragments from this GRanges object that overlap the blacklist regions.
fragment_spans <- granges(cnrReads_filter)
bl_remove <- overlapsAny(fragment_spans, bl_regions)
table(bl_remove)## bl_remove
## FALSE TRUE
## 394289 5893Write out filtered BAM file - CUT&RUN
This logical vector is then used to filter out these fragments that overlap blacklist.
It is helpful to include the names for each read in the BAM file for the peak callers, so we unlist the GAlignmentPairs object so each read now has its own row, and then the names are set from the metadata of the object.
Lastly, the BAM file is written to file.
Prepare BAM for peak calling - CUT&RUN
The samtools fixmate function is helpful to clean up
paired end BAM files. It will repair some of the metadata in the BAM
file, which can be especially helpful when doing filtering and
reading/writing of paired-end data in R.
Below we use Herper to call samtools as the fixmate function is not
available with Rsamtools. You’ll notice a few extra tricks in the
with_CondaEnv function call.
We can use brackets {} to run multiple lines of code, and we also
utilize the tempfile function to create a temporary path.
These files will be deleted with the R session ends.
forPeak_bam <- gsub("_filter.bam", "_forPeak.bam", filter_bam)
Herper::with_CondaEnv("CnR_analysis", pathToMiniConda = "miniconda", {
tempBam <- paste0(tempfile(), ".bam")
system(paste("samtools sort", "-n", "-o", tempBam, filter_bam, sep = " "))
system(paste("samtools fixmate", "-m", tempBam, forPeak_bam, sep = " "))
})Prepare BAM for peak calling - CUT&RUN
The peak caller MACS requires the file to be sorted by coordinate.
This can be done in R with the sortBam function.
Note: This could also be done using samtools sort on the
command line using Herper, as shown in the previous slide when we sorted
by name before using samtools fixmate.
Peak calling - CUT&RUN
CUT&RUN peak calling with MACS
MACS3 is installed and ready to use in our conda environment. The
macs3 callpeak function is used to call peaks.
To run MACS to we simply need to supply:
- A BAM file to find enriched regions in. (specified after -t)
- File format (BAMPE for paired end, specified after -f)
- A Name for peak calls (specified after –name).
- An output folder to write peaks into (specified after –outdir).
- Optionally, and sometimes recommended, we can identify a control file (we don’t have one, but would be specified after –c).
Example command line function call:
macs3 callpeak -t pairedEnd.bam -f BAMPE
--outdir path/to/output/
--n pairedEndPeakNameCUT&RUN peak calling with MACS
If we have sequenced paired-end data then we do know the fragment lengths and can provide BAM files to MACS which have been prefiltered to properly paired
We have to simply tell MACS that the data is paired using the format argument. This is important as MACS will guess it is single-end BAM by default.
MACS output files
Following MACS peak calling we would get 3 files.
- {name}.narrowPeak – Narrow peak format suitable for IGV and further analysis
- {name}.peaks.xls – Peak table suitable for review in excel.(not actually a xls but a tsv)
- {name}.summits.bed – Summit positions for peaks useful for finding motifs and plotting
CUT&RUN peak calling with SEACR
SEACR is a peak calling algorithm that was written by same lab that developed CUT&RUN (Henikoff). It is meant to deal with the lower read depths and sparser background signal from CUT&RUN.
SEACR requires a bedgraph file, which like a BAM file contains information about genomic CUT&RUN signal. In the following slides we will make a paired end bed file, filter it and then write out as a bedgraph.
SEACR is built for paired end data, and first we make a bedpe file
using the function bamtobed from bedtools. We can use our
conda environment we made previously to do this.
CUT&RUN peak calling with SEACR
Lastly we make a bedgraph file with the genomecov
function from bedtools.
A key input to this function is a file with chromosome lengths, and we can extract this from the FASTA file.
NOTE: The chrom.lengths.txt is already available in the ‘data’ folder of the course.
library(dplyr)
library(tibble)
indexFa("BSgenome.Mmusculus.UCSC.mm10.mainChrs.fa")
seqlengths(Rsamtools::FaFile("BSgenome.Mmusculus.UCSC.mm10.mainChrs.fa")) %>%
as.data.frame() %>%
rownames_to_column() %>%
write.table(file = "chrom.lengths.txt", row.names = FALSE, col.names = FALSE,
quote = FALSE, sep = "\t")CUT&RUN peak calling with SEACR
Then we make the bedgraph file with genomecov using the
file containing chromosome lengths in the ‘-g’ argument and the bedpe
file in the ‘-i’ argument.
The result is written to a file that is our bedgraph.
CUT&RUN peak calling with SEACR
Now we have the file required to run SEACR, but we need the SEACR executable script to call peaks.
We can download this from the Github page for SEACR. After it’s downloaded you might see the following error when running the script, indicating the permissions need to be changed.
download.file("https://github.com/FredHutch/SEACR/archive/refs/tags/v1.4-beta.2.zip",
destfile = "~/Downloads/SEACR_v1.4.zip")
unzip("~/Downloads/SEACR_v1.4.zip", exdir = "~/Downloads/SEACR_v1.4")
seacr_path <- "~/Downloads/SEACR_v1.4/SEACR-1.4-beta.2/SEACR_1.4.sh"
system(paste(seacr_path, "-h"))trying URL 'https://github.com/FredHutch/SEACR/archive/refs/tags/v1.4-beta.2.zip'
downloaded 3.0 MB
sh: /Users/douglasbarrows/Downloads/SEACR_v1.4/SEACR-1.4-beta.2/SEACR_1.4.sh: Permission deniedCUT&RUN peak calling with SEACR
The permissions can be changed by running the chmod
function. Here we make it accessible to anyone, but this can be
fine-tuned depending on your system.
SEACR: Sparse Enrichment Analysis for CUT&RUN
Usage: bash SEACR_1.4.sh -b <experimental bedgraph>.bg -c [<control bedgraph>.bg | <FDR threshold>] -n [norm | non] -m [relaxed | stringent] -o output_prefix -e [double] -r [yes | no]
Required input fields:
-b|--bedgraph Experimental bedgraph file
-c|--control Control bedgraph file
-n|--normalize Internal normalization (norm|non)
-m|--mode Stringency mode (stringent|relaxed)
-o|--output Output prefix
Optional input fields:
-e|--extension Peak extension constant for peak merging (Default=0.1)
-r|--remove Remove peaks overlapping IgG peaks (yes|no, Default=yes)
-h|--help Print help screen
Output file:
<output prefix>.[stringent | relaxed].bed (Bed file of enriched regions)
Output data structure:
<chr> <start> <end> <AUC> <max signal> <max signal region>
Description of output fields:
Field 1: Chromosome
Field 2: Start coordinate
Field 3: End coordinate
Field 4: Total signal contained within denoted coordinates
Field 5: Maximum bedgraph signal attained at any base pair within denoted coordinates
Field 6: Region representing the farthest upstream and farthest downstream bases within the denoted coordinates that are represented by the maximum bedgraph signal
Examples:
bash SEACR_1.4.sh -b target.bedgraph -c IgG.bedgraph -n norm -m stringent -o output
Calls enriched regions in target data using normalized IgG control track with stringent threshold and default 10% peak extension for merging
bash SEACR_1.4.sh -b target.bedgraph -c IgG.bedgraph -n non -m relaxed -o output -e 0.5
Calls enriched regions in target data using non-normalized IgG control track with relaxed threshold and 50% peak extension for merging
bash SEACR_1.4.sh -b target.bedgraph -c 0.01 -n non -m stringent -o output
Calls enriched regions in target data by selecting the top 1% of regions by total signal
CUT&RUN peak calling with SEACR
Now we can call peaks with the bedgraph file we generated and the SEACR.sh script.
Here we set a few key arguments for peaks calling: * -b: the path to the bedgraph file * -c: either a path to a control file (eg IgG) or a number between 0 and 1. If its a number, then this fraction will control the top peaks returned based on the signal in the peak + below we will only return the top 1% of peaks * -n: ‘norm’ if using control, ‘non’ if not * -m: ‘relaxed’ or ‘stringent’ mode
CUT&RUN peak calling with SEACR
Following peak calling we get one file.
- {prefix}.stringent.bed – a bed file with columns showing signal information that is suitable for IGV and further analysis
We can look at this peak file in IGV and compare to the MACS peak calls:
Post alignment processing and peak calling - ATACseq
Proper pairs - ATACseq
ATACseq processing and peak calling is very similar to CUT&RUN, with minor differences. We will quickly go through the same workflow for ATACseq.
flags = scanBamFlag(isProperPair = TRUE)
myParam = ScanBamParam(flag = flags, what = c("qname", "mapq", "isize", "flag"))
sortedBAM_atac <- "data/W6_ATAC_rep1_chr18_sorted.bam"
atacReads <- readGAlignmentPairs(sortedBAM_atac, param = myParam)
head(atacReads, 2)## GAlignmentPairs object with 2 pairs, strandMode=1, and 0 metadata columns:
## seqnames strand : ranges -- ranges
## <Rle> <Rle> : <IRanges> -- <IRanges>
## [1] chr18 - : 3000146-3000189 -- 3000146-3000189
## [2] chr18 + : 3000257-3000300 -- 3000280-3000325
## -------
## seqinfo: 44 sequences from an unspecified genomePlotting insert sizes - ATACseq
ATACseq should represent a mix of fragment lengths corresponding to nucleosome free, mononucleosome and polynucleosome fractions. This is different from TF CUT&RUN where we expect almost exclusively short fragements.
We can use the table() function to retrieve a vector of the occurrence of each fragment length.
atacReads_read1 <- GenomicAlignments::first(atacReads)
atacReads_read2 <- GenomicAlignments::second(atacReads)
insertSizes_atac <- abs(elementMetadata(atacReads_read1)$isize)
fragLenSizes_atac <- table(insertSizes_atac)
toPlot_atac <- data.frame(InsertSize = as.numeric(names(fragLenSizes_atac)), Count = as.numeric(fragLenSizes_atac))
fragLenPlot_atac <- ggplot(toPlot_atac, aes(x = InsertSize, y = Count)) + geom_line() +
theme_bw()Plotting insert sizes - ATACseq
The nucleosome pattern becomes much more clear when we apply a log2 transformation to the counts.
A good quality ATACseq experiment typically shows this nucleosome patterning.
Plotting insert sizes - ATACseq
We can now annotate our nucleosome free (< 100bp), mononucleosome (180bp-247bp) and dinucleosome (315-437) as in the Greenleaf study.
fragLenPlot_atac + scale_y_continuous(trans = "log2") + geom_vline(xintercept = c(180,
247), colour = "red") + geom_vline(xintercept = c(315, 437), colour = "darkblue") +
geom_vline(xintercept = c(100), colour = "darkgreen") + theme_bw()
Filter reads by MAPQ - ATACseq
We will filter out fragments where either read has a MAPQ of 0. For CUT&RUN we also filtered out longer fragments, we will not do that here to get a full picture of the chromatin structure.
Remove blacklist regions - ATACseq
The granges function finds the range spanning the whole
fragment (5’ read 1 to 3’ read 2) when given a GAlignmentPairs
object.
We then identify the fragments from this GRanges object that overlap the blacklist regions.
Write out filtered BAM file - ATACseq
This logical vector is then used to filter out these fragments that overlap blacklist.
It is helpful to include the names for each read in the BAM file for the peak callers, so we unlist the GAlignmentPairs object so each read now has its own row, and then the names are set from the metadata of the object.
Lastly, the BAM file is written to file.
atacReads_filter_noBL <- atacReads_filter[!bl_remove_atac]
atacReads_unlist <- unlist(atacReads_filter_noBL)
names(atacReads_unlist) <- mcols(atacReads_unlist)$qname
filter_bam_atac <- gsub("sorted.bam", "filter.bam", basename(sortedBAM_atac))
rtracklayer::export(atacReads_unlist, filter_bam_atac, format = "bam")Prepare BAM for peak calling - ATACseq
Similar to CUT&RUN, the BAM file needs to be cleaned up and prepared for MACS3. Here we do this using Herper and the command line version of samtools to sort by name, fix the metadata with fixmate, and lastly sort by coordinate.
forPeak_bam_atac <- gsub("_filter.bam", "_macs.bam", filter_bam_atac)
Herper::with_CondaEnv("CnR_analysis", pathToMiniConda = "miniconda", {
tempBam <- paste0(tempfile(), ".bam")
tempBam2 <- paste0(tempfile(), ".bam")
system(paste("samtools sort", "-n", "-o", tempBam, filter_bam_atac, sep = " "))
system(paste("samtools fixmate", "-m", tempBam, tempBam2, sep = " "))
system(paste("samtools sort -o", forPeak_bam_atac, tempBam2, sep = " "))
})ATACseq peak calling with MACS
Then we call peaks with MACS3 using the BAMPE setting.
While its not common to perform ATAC-seq using single end sequencing, there are some specific strategies to call peaks for single end data. In this case, see our extended ATACseq course.
ATACseq peak calling with MACS
BAM file for Nucleosome free regions
It is also common to only use nucleosome free regions (e.g. < 100bp) to focus on TF binding sites. This can be helpful for some down stream analyses, such as motif analysis.
In this case, we filter the GAlignmentPairs object based on isnert size.
read1_noBL <- GenomicAlignments::first(atacReads_filter_noBL)
insertSizes_atac_noBL <- abs(elementMetadata(read1_noBL)$isize)
atacReads_NF <- atacReads_filter_noBL[insertSizes_atac_noBL < 100]
atacReads_NF_unlist <- unlist(atacReads_NF)
names(atacReads_NF_unlist) <- mcols(atacReads_NF_unlist)$qname
NF_bam_atac <- gsub("sorted.bam", "filterNF.bam", basename(sortedBAM_atac))
rtracklayer::export(atacReads_NF_unlist, NF_bam_atac, format = "bam")NF BAM for peak calling - ATACseq
As before, the BAM file needs to be cleaned up and prepared for MACS3.
forPeak_NFbam_atac <- gsub("_filterNF.bam", "_macsNF.bam", NF_bam_atac)
Herper::with_CondaEnv("CnR_analysis", pathToMiniConda = "miniconda", {
tempBam <- paste0(tempfile(), ".bam")
tempBam2 <- paste0(tempfile(), ".bam")
system(paste("samtools sort", "-n", "-o", tempBam, NF_bam_atac, sep = " "))
system(paste("samtools fixmate", "-m", tempBam, tempBam2, sep = " "))
system(paste("samtools sort -o", forPeak_NFbam_atac, tempBam2, sep = " "))
})Call NF peaks with MACS
Then we call peaks for the nucleosome free regions with MACS3 using the BAMPE setting.
peaksNF_name_atac <- gsub("_macsNF.bam", "_macsNF", basename(forPeak_NFbam_atac))
with_CondaEnv("CnR_analysis", pathToMiniConda = "miniconda", system2(command = "macs3",
args = c("callpeak", "-t", forPeak_NFbam_atac, "--format", "BAMPE", "--outdir",
".", "--name", peaksNF_name_atac)))
Time for an exercise!
Exercise on CUT&RUN data can be found here
Answers to exercise
Answers can be found here