Last updated: 2022-03-16

Checks: 7 0

Knit directory: chipseq-cross-species/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20220209)

The command set.seed(20220209) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: c447893

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version c447893. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .RData
    Ignored:    .Rhistory
    Ignored:    analysis/.RData
    Ignored:    analysis/.Rhistory
    Ignored:    data/genomic_data/
    Ignored:    output/annotations/
    Ignored:    output/bam_files/
    Ignored:    output/filtered_peaks/
    Ignored:    output/logs/
    Ignored:    output/peaks/
    Ignored:    output/plots/
    Ignored:    output/qc/A-1_input.SeqDepthNorm.bw
    Ignored:    output/qc/A-1_input.SeqDepthNorm_dunnart_downSampled.bw
    Ignored:    output/qc/A-1_input.ccurve.txt
    Ignored:    output/qc/A-1_input.extrap.txt
    Ignored:    output/qc/A-1_input_R1_trimmed.fastq.gz
    Ignored:    output/qc/A-1_input_est_lib_complex_metrics.txt
    Ignored:    output/qc/A-2_H3K4me3.SeqDepthNorm.bw
    Ignored:    output/qc/A-2_H3K4me3.SeqDepthNorm_dunnart_downSampled.bw
    Ignored:    output/qc/A-2_H3K4me3.ccurve.txt
    Ignored:    output/qc/A-2_H3K4me3.extrap.txt
    Ignored:    output/qc/A-2_H3K4me3_R1_trimmed.fastq.gz
    Ignored:    output/qc/A-2_H3K4me3_est_lib_complex_metrics.txt
    Ignored:    output/qc/A-2_H3K4me3_vs_A-1_input.frip_default.txt
    Ignored:    output/qc/A-2_H3K4me3_vs_A-1_input.frip_default_dunnart_downSampled.txt
    Ignored:    output/qc/A-2_H3K4me3_vs_B-1_input.frip_default.txt
    Ignored:    output/qc/A-2_H3K4me3_vs_B-1_input.frip_default_dunnart_downSampled.txt
    Ignored:    output/qc/A-3_H3K27ac.SeqDepthNorm.bw
    Ignored:    output/qc/A-3_H3K27ac.SeqDepthNorm_dunnart_downSampled.bw
    Ignored:    output/qc/A-3_H3K27ac.ccurve.txt
    Ignored:    output/qc/A-3_H3K27ac.extrap.txt
    Ignored:    output/qc/A-3_H3K27ac_R1_trimmed.fastq.gz
    Ignored:    output/qc/A-3_H3K27ac_est_lib_complex_metrics.txt
    Ignored:    output/qc/A-3_H3K27ac_vs_A-1_input.frip_default.txt
    Ignored:    output/qc/A-3_H3K27ac_vs_A-1_input.frip_default_dunnart_downSampled.txt
    Ignored:    output/qc/A-3_H3K27ac_vs_B-1_input.frip_default.txt
    Ignored:    output/qc/A-3_H3K27ac_vs_B-1_input.frip_default_dunnart_downSampled.txt
    Ignored:    output/qc/B-1_input.SeqDepthNorm.bw
    Ignored:    output/qc/B-1_input.SeqDepthNorm_dunnart_downSampled.bw
    Ignored:    output/qc/B-1_input.ccurve.txt
    Ignored:    output/qc/B-1_input.extrap.txt
    Ignored:    output/qc/B-1_input_R1_trimmed.fastq.gz
    Ignored:    output/qc/B-1_input_est_lib_complex_metrics.txt
    Ignored:    output/qc/B-2_H3K4me3.SeqDepthNorm.bw
    Ignored:    output/qc/B-2_H3K4me3.SeqDepthNorm_dunnart_downSampled.bw
    Ignored:    output/qc/B-2_H3K4me3.ccurve.txt
    Ignored:    output/qc/B-2_H3K4me3.extrap.txt
    Ignored:    output/qc/B-2_H3K4me3_R1_trimmed.fastq.gz
    Ignored:    output/qc/B-2_H3K4me3_est_lib_complex_metrics.txt
    Ignored:    output/qc/B-2_H3K4me3_vs_A-1_input.frip_default.txt
    Ignored:    output/qc/B-2_H3K4me3_vs_A-1_input.frip_default_dunnart_downSampled.txt
    Ignored:    output/qc/B-2_H3K4me3_vs_B-1_input.frip_default.txt
    Ignored:    output/qc/B-2_H3K4me3_vs_B-1_input.frip_default_dunnart_downSampled.txt
    Ignored:    output/qc/B-3_H3K27ac.SeqDepthNorm.bw
    Ignored:    output/qc/B-3_H3K27ac.SeqDepthNorm_dunnart_downSampled.bw
    Ignored:    output/qc/B-3_H3K27ac.ccurve.txt
    Ignored:    output/qc/B-3_H3K27ac.extrap.txt
    Ignored:    output/qc/B-3_H3K27ac_R1_trimmed.fastq.gz
    Ignored:    output/qc/B-3_H3K27ac_est_lib_complex_metrics.txt
    Ignored:    output/qc/B-3_H3K27ac_vs_A-1_input.frip_default.txt
    Ignored:    output/qc/B-3_H3K27ac_vs_A-1_input.frip_default_dunnart_downSampled.txt
    Ignored:    output/qc/B-3_H3K27ac_vs_B-1_input.frip_default.txt
    Ignored:    output/qc/B-3_H3K27ac_vs_B-1_input.frip_default_dunnart_downSampled.txt
    Ignored:    output/qc/H3K27ac_multiBAM_fingerprint_metrics.txt
    Ignored:    output/qc/H3K27ac_multiBAM_fingerprint_rawcounts.txt
    Ignored:    output/qc/H3K27ac_multibamsum.npz
    Ignored:    output/qc/H3K27ac_multibamsum.tab
    Ignored:    output/qc/H3K27ac_pearsoncor_multibamsum_matrix.txt
    Ignored:    output/qc/H3K27ac_plot_coverage_rawcounts.tab
    Ignored:    output/qc/bamPEFragmentSize_rawcounts.tab
    Ignored:    output/qc/bamPEFragmentSize_rawcounts_dunnart_downSampled.tab
    Ignored:    output/qc/multiBAM_fingerprint_metrics.txt
    Ignored:    output/qc/multiBAM_fingerprint_metrics_dunnart_downSampled.txt
    Ignored:    output/qc/multiBAM_fingerprint_rawcounts.txt
    Ignored:    output/qc/multiBAM_fingerprint_rawcounts_dunnart_downSampled.txt
    Ignored:    output/qc/multibamsum.npz
    Ignored:    output/qc/multibamsum.tab
    Ignored:    output/qc/multibamsum_dunnart_downSampled.npz
    Ignored:    output/qc/multibamsum_dunnart_downSampled.tab
    Ignored:    output/qc/pearsoncor_multibamsum_matrix.txt
    Ignored:    output/qc/pearsoncor_multibamsum_matrix_dunnart_downSampled.txt
    Ignored:    output/wga/
    Ignored:    results/

Untracked files:
    Untracked:  .snakemake/
    Untracked:  Rplots.pdf
    Untracked:  data/raw_reads/
    Untracked:  output/qc/H3K27ac_multiBAM_fingerprint.pdf
    Untracked:  output/qc/H3K27ac_overlap_default_dunnart_downSampled.frip
    Untracked:  output/qc/H3K27ac_pearsoncor_multibamsum.pdf
    Untracked:  output/qc/H3K27ac_plot_coverage.pdf
    Untracked:  output/qc/H3K4me3_overlap_default_dunnart_downSampled.frip
    Untracked:  output/qc/ucsc_alignment/
    Untracked:  output/rnaseq/

Unstaged changes:
    Deleted:    output/qc/H3K4me3_overlap_default.frip

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/mouse_dunnart_data_processing_for_comparison.Rmd) and HTML (docs/mouse_dunnart_data_processing_for_comparison.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	c447893	lecook	2022-03-16	wflow_publish(“analysis/mouse_dunnart_data_processing_for_comparison.Rmd”)
html	a90d5a0	lecook	2022-03-16	Build site.
Rmd	87c4146	lecook	2022-03-16	wflow_publish(“analysis/mouse_dunnart_data_processing_for_comparison.Rmd”)
html	93879c2	lecook	2022-03-16	Build site.
html	a552f17	lecook	2022-03-16	Build site.
Rmd	5a17b44	lecook	2022-03-16	wflow_publish(“analysis/mouse_dunnart_data_processing_for_comparison.Rmd”)
html	11fbc81	lecook	2022-03-16	Build site.
Rmd	dffebdd	lecook	2022-03-16	wflow_publish(“analysis/mouse_dunnart_data_processing_for_comparison.Rmd”)
html	4b71068	lecook	2022-03-16	Build site.
Rmd	a287d1b	lecook	2022-03-16	wflow_publish(“analysis/mouse_dunnart_data_processing_for_comparison.Rmd”)
Rmd	eaeea84	lecook	2022-03-16	updated
Rmd	81d5677	lecook	2022-03-16	fix downsampled
html	c5d4b95	lecook	2022-03-16	Build site.
html	f7966c8	lecook	2022-03-16	Build site.
Rmd	040a6a7	lecook	2022-03-16	wflow_publish(“analysis/mouse_dunnart_data_processing_for_comparison.Rmd”)
html	f593bed	lecook	2022-03-16	Build site.
Rmd	039df2a	lecook	2022-03-16	wflow_publish(“analysis/mouse_dunnart_data_processing_for_comparison.Rmd”)
html	9d125bb	lecook	2022-03-16	Build site.
html	9bcce4a	lecook	2022-03-16	Build site.
Rmd	c0a1798	lecook	2022-03-16	wflow_publish(“analysis/mouse_dunnart_data_processing_for_comparison.Rmd”)
html	841e0f4	lecook	2022-03-16	Build site.
Rmd	275510f	lecook	2022-03-16	wflow_publish(“analysis/mouse_dunnart_data_processing_for_comparison.Rmd”)
Rmd	faabd95	lecook	2022-03-03	updated
Rmd	2103ae7	lecook	2022-03-03	updates
Rmd	323d8b6	lecook	2022-03-03	copied scripts over
Rmd	ae523af	lecook	2022-03-02	fix code formatting for bam coverage rule
Rmd	a97d3c5	lecook	2022-03-01	first commit
html	275d229	lecook	2022-03-01	Build site.
Rmd	784c4ec	lecook	2022-03-01	update
Rmd	bf8fac9	lecook	2022-03-01	updated
Rmd	5fb3661	lecook	2022-03-01	clean up
Rmd	c2eeed7	lecook	2022-02-28	debugged

Pipeline for mouse and dunnart peak calling

Download mouse unfiltered alignments from ENCODE

See mouse_data_ENCODE for details on these samples.

xargs -L 1 curl -O -J -L < ENCODE_files.txt
xargs -L 1 curl -O -J -L < control_files.txt

Filter mouse alignments

Uses samtools and picards MarkDuplicates to remove low quality reads. This is based on and in accordance with the ENCODE Guidelines and pipeline (https://github.com/ENCODE-DCC/chip-seq-pipeline2)

Library Complexity ChIP-seq Standards:

PBC1	PBC2	Bottlenecking level	NRF	Complexity	Flag colors
< 0.5	< 1	Severe	< 0.5	Concerning	Orange
0.5 ≤ PBC1 < 0.8	1 ≤ PBC2 < 3	Moderate	0.5 ≤ NRF < 0.8	Acceptable	Yellow
0.8 ≤ PBC1 < 0.9	3 ≤ PBC2 < 10	Mild	0.8 ≤ NRF < 0.9	Compliant	None
≥ 0.9	≥ 10	None	> 0.9	Ideal	None

See full snakemake script for details: code/mouse_peak_calling/mouse_H3K4me3 and code/mouse_peak_calling/mouse_H3K27ac

Table. Read alignment quality metrics before and aftering filtering steps

				read counts
replicate	antibody	stage	accession ID	bwa alignment	filtered readŝ	remove duplicates	NRF*	PBC1*	PBC2*
1	H3K27ac	E10.5	ENCFF213EBC	86919825	74296920	61061524	0.83	0.83	5.53
2	H3K27ac	E10.5	ENCFF548BRR	49113729	43572898	38069875	0.88	0.88	8.03
1	H3K27ac	E11.5	ENCFF512SFE	1714339	12827968	11879316	0.93	0.93	13.76
2	H3K27ac	E11.5	ENCFF515PKL	18419756	13753883	11961286	0.88	0.88	7.7
2	H3K27ac	E12.5	ENCFF011NFM	41510259	32754035	28442356	0.88	0.88	7.55
1	H3K27ac	E13.5	ENCFF194ORC	47508791	39115139	23235076	0.64	0.62	2.49
2	H3K27ac	E13.5	ENCFF290ZNF	46476048	37416425	26624266	0.74	0.73	3.52
1	H3K27ac	E14.5	ENCFF327VAO	37480089	30104464	19228162	0.66	0.64	2.53
2	H3K27ac	E14.5	ENCFF902HAR	27840473	22156787	15203030	0.69	0.67	2.68
1	H3K27ac	E15.5	ENCFF584JFB	51675307	40125566	35551034	0.9	0.9	9.61
2	H3K27ac	E15.5	ENCFF707WKL	38499223	30334974	27522173	0.92	0.92	11.76
1	H3K4me3	E10.5	ENCFF124UYX	74283697	65480729	47429840	0.73	0.76	4.51
2	H3K4me3	E10.5	ENCFF045IPK	39305102	34283317	27327155	0.8	0.82	5.7
1	H3K4me3	E11.5	ENCFF760QYZ	24048599	17547752	16515484	0.95	0.95	20.21
2	H3K4me3	E11.5	ENCFF717QDV	19425608	14136366	13504654	0.96	0.97	26.71
1	H3K4me3	E12.5	ENCFF182ZPF	37674883	30140587	26489691	0.89	0.89	9.29
2	H3K4me3	E12.5	ENCFF941QJZ	38463605	29985012	26459778	0.89	0.9	9.44
1	H3K4me3	E13.5	ENCFF485UDC	47478873	38287309	33094082	0.89	0.9	10.81
2	H3K4me3	E13.5	ENCFF124TAB	47598002	38139564	33020017	0.89	0.91	11.12
1	H3K4me3	E14.5	ENCFF724DMU	25402316	18964886	17278836	0.92	0.93	13.9
2	H3K4me3	E14.5	ENCFF665QBJ	28852674	22484958	19804325	0.9	0.9	10.56
1	H3K4me3	E15.5	ENCFF258KCR	55588923	44138244	38278946	0.89	0.91	12.05
2	H3K4me3	E15.5	ENCFF401BKM	42617927	34005723	29687861	0.9	0.91	12.13
1	input control	E10.5	ENCFF157KEH	75962416	57382820	55915261	0.98	0.98	45.2
2	input control	E10.5	ENCFF825AVI	70250734	53358758	51796823	0.98	0.98	38.64
1	input control	E11.5	ENCFF184CUE	14581637	10190363	10031334	0.99	0.99	69.67
2	input control	E11.5	ENCFF376FGM	16001082	11239845	11033934	0.99	0.99	59.63
1	input control	E12.5	ENCFF203JQV	74367126	54966177	52895081	0.97	0.97	29.47
2	input control	E12.5	ENCFF058AUT	68220307	50400902	48480815	0.97	0.97	29.1
1	input control	E13.5	ENCFF117QRC	77084648	56275163	53582563	0.97	0.97	30.93
2	input control	E13.5	ENCFF248PGK	66492732	48343406	46564874	0.98	0.98	38.51
1	input control	E14.5	ENCFF784ORI	35480589	25431770	24573359	0.97	0.97	32.34
2	input control	E14.5	ENCFF002HZV	33021720	23658223	23199516	0.99	0.99	57.52
1	input control	E15.5	ENCFF727QTS	72318210	51828288	50444903	0.98	0.99	52.73
2	input control	E15.5	ENCFF182XFG	86322714	63335464	60305840	0.97	0.97	28.35

calculated from a subsample of ~10-15M aligned reads prior to removing MarkDuplicates

̂ excluded low quality (MAPQ30) and orphaned reads

Normalise dunnart and mouse BAM files to 10M reads

This is so we can compare to mouse peaks.

Based on https://davemcg.github.io/post/easy-bam-downsampling/ script for doing this.

Run script in the directory with the files you want to subsample.

cd output/bam_files/
bash subsample.sh

Read QC with deepTools

Plot correlation between BAM files

plotCorrelation computes the overall similarity between two or more files based on read coverage (or other scores) within genomic regions.

This helps to determine whether the different sample types can be separated, i.e., samples of different conditions are expected to be more dissimilar to each other than replicates within the same condition.

Fingerprint plots

This tool samples indexed BAM files and plots a profile of cumulative read coverages for each. All reads overlapping a window (bin) of the specified length are counted; these counts are sorted and the cumulative sum is finally plotted.

It determines how well the signal in the ChIP-seq sample can be differentiated from the background distribution of reads in the control sample.

An ideal input with perfect uniform distribution of reads along the genome (i.e. without enrichments in open chromatin etc.) and infinite sequencing coverage should generate a straight diagonal line. A very specific and strong ChIP enrichment will be indicated by a prominent and steep rise of the cumulative sum towards the highest rank. This means that a big chunk of reads from the ChIP sample is located in few bins which corresponds to high, narrow enrichments typically seen for transcription factors.

Cross-correlation analysis

Information from: https://docs.google.com/document/d/1lG_Rd7fnYgRpSIqrIfuVlAz2dW1VaSQThzk836Db99c/edit

This set of programs operate on mapped Illumina single-end read datasets in tagAlign or BAM format. Because my data is paired-end I need to only use the forward read.

A high-quality ChIP-seq experiment will produce significant clustering of enriched DNA sequence tags/reads at locations bound by the protein of interest; the expectation is that we can observe a bimodal enrichment of reads (sequence tags) on both the forward and the reverse strands.

Cross-correlation analysis is done on a filtered (but not-deduped) and subsampled BAM. There is a special fastq trimming for cross-correlation analysis. Read1 fastq is trimmed to 50bp first using trimfastq.py (last modified 2017/11/08, https://github.com/ENCODE-DCC/chip-seq-pipeline2/blob/master/src/trimfastq.py). And then it is separately mapped as SE. Reads are filtered but duplicates are not removed. Then 15 million reads are randomly sampled and used for cross-correlation analysis.

ENCODE standards:

col.	abbreviation	description
1	Filename	tagAlign/BAM filename
2	numReads	effective sequencing depth i.e. total number of mapped reads in input file
3	estFragLen	comma separated strand cross-correlation peak(s) in decreasing order of correlation.
4	corr_estFragLen	comma separated strand cross-correlation value(s) in decreasing order (COL2 follows the same order)
5	phantomPeak	Read length/phantom peak strand shift
6	corr_phantomPeak	Correlation value at phantom peak
7	argmin_corr	strand shift at which cross-correlation is lowest
8	min_corr	minimum value of cross-correlation
9	NSC	Normalized strand cross-correlation coefficient (NSC) = COL4 / COL8
10	RSC	Relative strand cross-correlation coefficient (RSC) = (COL4 - COL8) / (COL6 - COL8)
11	QualityTag	Quality tag based on thresholded RSC (codes= -2:veryLow, -1:Low, 0:Medium, 1:High, 2:veryHigh)

NSC; NSC>1.1 (higher values indicate more enrichment; 1 = no enrichment)

RSC; RSC>0.8 (0 = no signal; <1 low quality ChIP; >1 high enrichment

Quality tag based on thresholded RSC (codes: -2:veryLow,-1:Low,0:Medium,1:High; 2:veryHigh)

Table. Peak calling quality metrics for mouse ChIP-seq peaks after peak calling with MACS2

replicate	antibody	stage	accession ID	number of peaks	FRiP ̂	NSC*	RSC*	Quality Tag*	% replicate reproducibility	number of consensus peaks	FRiP^ for consensus peaks
1	H3K27ac	E10.5	ENCFF213EBC	61169	0.12	1.11	1.15	1	45.27	27688	0.14
2	H3K27ac	E10.5	ENCFF548BRR	80654	0.19	1.19	1.23	1	34.33
1	H3K27ac	E11.5	ENCFF512SFE	25352	0.07	1.08	1.61	2	47.93	12149	0.06
2	H3K27ac	E11.5	ENCFF515PKL	23732	0.06	1.09	1.97	2	51.2
1	H3K27ac	E12.5	ENCFF394TZN	27242	0.07	1.08	1.59	2	58.17	15845	0.08
2	H3K27ac	E12.5	ENCFF011NFM	62319	0.11	1.09	1.49	1	25.43
1	H3K27ac	E13.5	ENCFF194ORC	85874	0.18	1.19	2.14	2	35.72	30671	0.13
2	H3K27ac	E13.5	ENCFF290ZNF	71246	0.13	1.13	2.05	2	43.05
1	H3K27ac	E14.5	ENCFF327VAO	45785	0.15	1.19	1.65	2	48.56	22232	0.14
2	H3K27ac	E14.5	ENCFF902HAR	44042	0.14	1.16	1.72	2	50.48
1	H3K27ac	E15.5	ENCFF584JFB	35724	0.09	1.08	1.49	1	51.08	18246	0.09
2	H3K27ac	E15.5	ENCFF707WKL	37600	0.1	1.1	1.54	2	48.53
1	H3K4me3	E10.5	ENCFF124UYX	29478	0.44	2.28	1.12	1	70.77	20860	0.47
2	H3K4me3	E10.5	ENCFF045IPK	26820	0.49	2.47	1.15	1	77.78
1	H3K4me3	E11.5	ENCFF760QYZ	28388	0.34	1.86	1.25	1	66.43	18857	0.33
2	H3K4me3	E11.5	ENCFF717QDV	28368	0.33	1.81	1.24	1	66.48
1	H3K4me3	E12.5	ENCFF182ZPF	29466	0.33	1.93	1.21	1	66.35	19549	0.32
2	H3K4me3	E12.5	ENCFF941QJZ	30868	0.31	1.67	1.21	1	63.34
1	H3K4me3	E13.5	ENCFF485UDC	29715	0.4	2.16	1.24	1	69.78	20734	0.41
2	H3K4me3	E13.5	ENCFF124TAB	29547	0.41	2.19	1.24	1	70.18
1	H3K4me3	E14.5	ENCFF724DMU	24250	0.5	2.34	1.26	1	81.08	19660	0.53
2	H3K4me3	E14.5	ENCFF665QBJ	25084	0.55	2.54	1.26	1	78.38
1	H3K4me3	E15.5	ENCFF258KCR	23365	0.36	2.25	1.29	1	80.11	18717	0.39
2	H3K4me3	E15.5	ENCFF401BKM	23920	0.42	2.47	1.33	1	78.25

* computed from 15M reads subsampled ̂ fraction of reads in peaks

Peak QC

Number of peaks

Dunnart

echo "Number of peaks A3 H3K27ac dunnart"
less output/peaks/A-3_H3K27ac_vs_A-1_input_macs2_default_dunnart_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks A2 H3K4me3 dunnart"
less output/peaks/A-2_H3K4me3_vs_A-1_input_macs2_default_dunnart_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks B3 H3K27ac dunnart"
less output/peaks/B-3_H3K27ac_vs_B-1_input_macs2_default_dunnart_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks B2 H3K4me3 dunnart"
less output/peaks/B-2_H3K4me3_vs_B-1_input_macs2_default_dunnart_downSampled_peaks.narrowPeak | wc -l

Mouse

H3K4me3

echo "Number of peaks rep 1 E10.5 H3K4me3"
less output/peaks/ENCFF124UYX_vs_ENCFF157KEH_E10.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E10.5 H3K4me3"
less output/peaks/ENCFF045IPK_vs_ENCFF825AVI_E10.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E11.5 H3K4me3"
less output/peaks/ENCFF760QYZ_vs_ENCFF184CUE_E11.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E11.5 H3K4me3"
less output/peaks/ENCFF717QDV_vs_ENCFF376FGM_E11.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E12.5 H3K4me3"
less output/peaks/ENCFF941QJZ_vs_ENCFF058AUT_E12.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E12.5 H3K4me3"
less output/peaks/ENCFF182ZPF_vs_ENCFF203JQV_E12.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E13.5 H3K4me3"
less output/peaks/ENCFF485UDC_vs_ENCFF117QRC_E13.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E13.5 H3K4me3"
less output/peaks/ENCFF124TAB_vs_ENCFF248PGK_E13.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E14.5 H3K4me3"
less output/peaks/ENCFF665QBJ_vs_ENCFF002HZV_E14.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E14.5 H3K4me3"
less output/peaks/ENCFF724DMU_vs_ENCFF784ORI_E14.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E15.5 H3K4me3"
less output/peaks/ENCFF401BKM_vs_ENCFF182XFG_E15.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E15.5 H3K4me3"
less output/peaks/ENCFF258KCR_vs_ENCFF727QTS_E15.5_H3K4me3_macs2_downSampled_peaks.narrowPeak | wc -l

H3K27ac

echo "Number of peaks rep 1 E10.5 H3K27ac"
less output/peaks/ENCFF213EBC_vs_ENCFF157KEH_E10.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E10.5 H3K27ac"
less output/peaks/ENCFF548BRR_vs_ENCFF825AVI_E10.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E11.5 H3K27ac"
less output/peaks/ENCFF512SFE_vs_ENCFF184CUE_E11.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E11.5 H3K27ac"
less output/peaks/ENCFF515PKL_vs_ENCFF376FGM_E11.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E12.5 H3K27ac"
less output/peaks/ENCFF011NFM_vs_ENCFF058AUT_E12.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E12.5 H3K27ac"
less output/peaks/ENCFF394TZN_vs_ENCFF203JQV_E12.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E13.5 H3K27ac"
less output/peaks/ENCFF194ORC_vs_ENCFF117QRC_E13.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 2 E13.5 H3K27ac"
less output/peaks/ENCFF290ZNF_vs_ENCFF248PGK_E13.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E14.5 H3K27ac"
less output/peaks/ENCFF902HAR_vs_ENCFF002HZV_E14.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l 
echo "Number of peaks rep 2 E14.5 H3K27ac"
less output/peaks/ENCFF327VAO_vs_ENCFF784ORI_E14.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l
echo "Number of peaks rep 1 E15.5 H3K27ac"
less output/peaks/ENCFF707WKL_vs_ENCFF182XFG_E15.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l 
echo "Number of peaks rep 2 E15.5 H3K27ac"
less output/peaks/ENCFF584JFB_vs_ENCFF727QTS_E15.5_H3K27ac_macs2_downSampled_peaks.narrowPeak | wc -l

Consensus peaksets for biological replicates

Number of consensus peaks

Dunnart

echo "Number of consensus H3K4me3 dunnart peaks"
less output/peaks/H3K4me3_overlap_default_dunnart_downSampled_peaks.narrowPeak | wc -l
echo "Number of consensus H3K27ac dunnart peaks"
less output/peaks/H3K27ac_overlap_default_dunnart_downSampled_peaks.narrowPeak | wc -l

Mouse

echo "Number of consensus H3K4me3 mouse E10.5 peaks"
less output/peaks/H3K4me3_E10.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K4me3 mouse E11.5 peaks"
less output/peaks/H3K4me3_E11.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K4me3 mouse E12.5 peaks"
less output/peaks/H3K4me3_E12.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K4me3 mouse E13.5 peaks"
less output/peaks/H3K4me3_E13.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K4me3 mouse E14.5 peaks"
less output/peaks/H3K4me3_E14.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K4me3 mouse E15.5 peaks"
less output/peaks/H3K4me3_E15.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K27ac mouse E10.5 peaks"
less output/peaks/H3K27ac_E10.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K27ac mouse E11.5 peaks"
less output/peaks/H3K27ac_E11.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K27ac mouse E12.5 peaks"
less output/peaks/H3K27ac_E12.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K27ac mouse E13.5 peaks"
less output/peaks/H3K27ac_E13.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K27ac mouse E14.5 peaks"
less output/peaks/H3K27ac_E14.5_overlap.narrowPeak | wc -l
echo "Number of consensus H3K27ac mouse E15.5 peaks"
less output/peaks/H3K27ac_E15.5_overlap.narrowPeak | wc -l

Consensus peaks and grouping into putative enhancers and promoters

Find uniquely H3K4me3 sites (i.e. peaks that don’t overlap with H3K27ac):

conda activate chip
bedtools intersect -v -a output/peaks/H3K4me3_E10.5_overlap.narrowPeak -b output/peaks/H3K27ac_E10.5_overlap.narrowPeak > output/peaks/H3K4me3_E10.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K4me3_E11.5_overlap.narrowPeak -b output/peaks/H3K27ac_E11.5_overlap.narrowPeak > output/peaks/H3K4me3_E11.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K4me3_E12.5_overlap.narrowPeak -b output/peaks/H3K27ac_E12.5_overlap.narrowPeak > output/peaks/H3K4me3_E12.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K4me3_E13.5_overlap.narrowPeak -b output/peaks/H3K27ac_E13.5_overlap.narrowPeak > output/peaks/H3K4me3_E13.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K4me3_E14.5_overlap.narrowPeak -b output/peaks/H3K27ac_E14.5_overlap.narrowPeak > output/peaks/H3K4me3_E14.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K4me3_E15.5_overlap.narrowPeak -b output/peaks/H3K27ac_E15.5_overlap.narrowPeak > output/peaks/H3K4me3_E15.5_only.narrowPeak


bedtools intersect -v -a output/peaks/H3K4me3_overlap_default_dunnart_downSampled_peaks.narrowPeak -b output/peaks/H3K27ac_overlap_default_dunnart_downSampled_peaks.narrowPeak > output/peaks/H3K4me3_only_dunnart_downSampled.narrowPeak

Find unique H3K27ac sites (i.e. peaks that don’t overlap with H3K4me3):

bedtools intersect -v -a output/peaks/H3K27ac_E10.5_overlap.narrowPeak -b output/peaks/H3K4me3_E10.5_overlap.narrowPeak > output/peaks/H3K27ac_E10.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K27ac_E11.5_overlap.narrowPeak -b output/peaks/H3K4me3_E11.5_overlap.narrowPeak > output/peaks/H3K27ac_E11.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K27ac_E12.5_overlap.narrowPeak -b output/peaks/H3K4me3_E12.5_overlap.narrowPeak > output/peaks/H3K27ac_E12.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K27ac_E13.5_overlap.narrowPeak -b output/peaks/H3K4me3_E13.5_overlap.narrowPeak > output/peaks/H3K27ac_E13.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K27ac_E14.5_overlap.narrowPeak -b output/peaks/H3K4me3_E14.5_overlap.narrowPeak > output/peaks/H3K27ac_E14.5_only.narrowPeak
bedtools intersect -v -a output/peaks/H3K27ac_E15.5_overlap.narrowPeak -b output/peaks/H3K4me3_E15.5_overlap.narrowPeak > output/peaks/H3K27ac_E15.5_only.narrowPeak


bedtools intersect -v -a output/peaks/H3K27ac_overlap_default_dunnart_downSampled_peaks.narrowPeak -b output/peaks/H3K4me3_overlap_default_dunnart_downSampled_peaks.narrowPeak > output/peaks/H3K27ac_only_dunnart_downSampled.narrowPeak

Find peaks common between H3K27ac & H3K4me3 with a reciprocal overlap of at least 50%.

bedtools intersect -f 0.5 -r -a output/peaks/H3K4me3_E10.5_overlap.narrowPeak -b output/peaks/H3K27ac_E10.5_overlap.narrowPeak > output/peaks/H3K4me3_and_H3K27ac_E10.5.narrowPeak
bedtools intersect -f 0.5 -r -a output/peaks/H3K4me3_E11.5_overlap.narrowPeak -b output/peaks/H3K27ac_E11.5_overlap.narrowPeak > output/peaks/H3K4me3_and_H3K27ac_E11.5.narrowPeak
bedtools intersect -f 0.5 -r -a output/peaks/H3K4me3_E12.5_overlap.narrowPeak -b output/peaks/H3K27ac_E12.5_overlap.narrowPeak > output/peaks/H3K4me3_and_H3K27ac_E12.5.narrowPeak
bedtools intersect -f 0.5 -r -a output/peaks/H3K4me3_E13.5_overlap.narrowPeak -b output/peaks/H3K27ac_E13.5_overlap.narrowPeak > output/peaks/H3K4me3_and_H3K27ac_E13.5.narrowPeak
bedtools intersect -f 0.5 -r -a output/peaks/H3K4me3_E14.5_overlap.narrowPeak -b output/peaks/H3K27ac_E14.5_overlap.narrowPeak > output/peaks/H3K4me3_and_H3K27ac_E14.5.narrowPeak
bedtools intersect -f 0.5 -r -a output/peaks/H3K4me3_E15.5_overlap.narrowPeak -b output/peaks/H3K27ac_E15.5_overlap.narrowPeak > output/peaks/H3K4me3_and_H3K27ac_E15.5.narrowPeak

bedtools intersect -f 0.5 -r -a output/peaks/H3K4me3_overlap_default_dunnart_downSampled_peaks.narrowPeak -b output/peaks/H3K27ac_overlap_default_dunnart_downSampled_peaks.narrowPeak > output/peaks/H3K4me3_H3K27ac_dunnart_downSampled.narrowPeak

Combine peak files

cat output/peaks/H3K4me3_E10.5_only.narrowPeak output/peaks/H3K4me3_and_H3K27ac_E10.5.narrowPeak | sort -k 4,4 > output/peaks/E10.5_promoter_peaks.narrowPeak
cat output/peaks/H3K4me3_E11.5_only.narrowPeak output/peaks/H3K4me3_and_H3K27ac_E11.5.narrowPeak | sort -k 4,4 > output/peaks/E11.5_promoter_peaks.narrowPeak
cat output/peaks/H3K4me3_E12.5_only.narrowPeak output/peaks/H3K4me3_and_H3K27ac_E12.5.narrowPeak | sort -k 4,4 > output/peaks/E12.5_promoter_peaks.narrowPeak
cat output/peaks/H3K4me3_E13.5_only.narrowPeak output/peaks/H3K4me3_and_H3K27ac_E13.5.narrowPeak | sort -k 4,4 > output/peaks/E13.5_promoter_peaks.narrowPeak
cat output/peaks/H3K4me3_E14.5_only.narrowPeak output/peaks/H3K4me3_and_H3K27ac_E14.5.narrowPeak | sort -k 4,4 > output/peaks/E14.5_promoter_peaks.narrowPeak
cat output/peaks/H3K4me3_E15.5_only.narrowPeak output/peaks/H3K4me3_and_H3K27ac_E15.5.narrowPeak | sort -k 4,4 > output/peaks/E15.5_promoter_peaks.narrowPeak

cat output/peaks/H3K4me3_only_dunnart_downSampled.narrowPeak output/peaks/H3K4me3_H3K27ac_dunnart_downSampled.narrowPeak| sort -k 4,4 > output/peaks/dunnart_downSampled_promoter_peaks.narrowPeak


cat output/peaks/H3K27ac_E10.5_only.narrowPeak | sort -k 4,4 > output/peaks/E10.5_enhancer_peaks.narrowPeak
cat output/peaks/H3K27ac_E11.5_only.narrowPeak | sort -k 4,4 > output/peaks/E11.5_enhancer_peaks.narrowPeak
cat output/peaks/H3K27ac_E12.5_only.narrowPeak | sort -k 4,4 > output/peaks/E12.5_enhancer_peaks.narrowPeak
cat output/peaks/H3K27ac_E13.5_only.narrowPeak | sort -k 4,4 > output/peaks/E13.5_enhancer_peaks.narrowPeak
cat output/peaks/H3K27ac_E14.5_only.narrowPeak | sort -k 4,4 > output/peaks/E14.5_enhancer_peaks.narrowPeak
cat output/peaks/H3K27ac_E15.5_only.narrowPeak | sort -k 4,4 > output/peaks/E15.5_enhancer_peaks.narrowPeak
cat output/peaks/H3K27ac_only_dunnart_downSampled.narrowPeak | sort -k 4,4 > output/peaks/dunnart_downSampled_enhancer_peaks.narrowPeak

sessionInfo()

R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux

Matrix products: default
BLAS/LAPACK: /usr/local/easybuild-2019/easybuild/software/compiler/gcc/10.2.0/openblas/0.3.12/lib/libopenblas_haswellp-r0.3.12.so

locale:
 [1] LC_CTYPE=en_AU.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_AU.UTF-8        LC_COLLATE=en_AU.UTF-8    
 [5] LC_MONETARY=en_AU.UTF-8    LC_MESSAGES=en_AU.UTF-8   
 [7] LC_PAPER=en_AU.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] workflowr_1.7.0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8        jquerylib_0.1.4   pillar_1.6.1      compiler_4.1.0   
 [5] bslib_0.2.5.1     later_1.2.0       git2r_0.28.0      tools_4.1.0      
 [9] getPass_0.2-2     digest_0.6.27     jsonlite_1.7.2    evaluate_0.14    
[13] lifecycle_1.0.1   tibble_3.1.2      pkgconfig_2.0.3   rlang_1.0.2      
[17] cli_2.5.0         rstudioapi_0.13   yaml_2.2.1        xfun_0.23        
[21] httr_1.4.2        stringr_1.4.0     knitr_1.33        sass_0.4.0       
[25] fs_1.5.0          vctrs_0.3.8       rprojroot_2.0.2   glue_1.4.2       
[29] R6_2.5.0          processx_3.5.2    fansi_0.5.0       rmarkdown_2.8    
[33] callr_3.7.0       magrittr_2.0.1    whisker_0.4       ps_1.6.0         
[37] promises_1.2.0.1  ellipsis_0.3.2    htmltools_0.5.1.1 httpuv_1.6.1     
[41] utf8_1.2.1        stringi_1.6.2     crayon_1.4.1

Processing mouse and dunnart data for comparison

lecook

2022-02-23

Pipeline for mouse and dunnart peak calling

Download mouse unfiltered alignments from ENCODE

Filter mouse alignments

Normalise dunnart and mouse BAM files to 10M reads

Read QC with deepTools

Plot correlation between BAM files

Fingerprint plots

Cross-correlation analysis

Peak QC

Number of peaks

Dunnart

Mouse

Consensus peaksets for biological replicates

Number of consensus peaks

Dunnart

Mouse

Consensus peaks and grouping into putative enhancers and promoters