Scientific Writing from ENCODE Provenance

Generate publication-quality scientific writing from ENCODE analysis records. This skill integrates with data-provenance and cite-encode to auto-generate methods from logged pipeline runs. Every generated section follows rigorous scientific documentation standards -- complete reporting of all experimental and computational parameters with zero ambiguity.

When to Use

• User wants to write publication-ready methods sections, figure legends, or data availability statements
• User asks about "methods section", "figure legend", "scientific writing", or "manuscript preparation"
• User needs to auto-generate methods text from logged provenance/analysis steps
• User wants templates for supplementary tables, Key Resources Tables, or tool citation formatting
• Example queries: "write a methods section for my ChIP-seq analysis", "generate a figure legend for my heatmap", "format my data availability statement"

Overview

Most methods sections in genomics papers are incomplete. They omit software versions, skip reference file details, conflate technical and biological replicates, and use phrases like "default parameters" without stating what those defaults are. Reviewers catch these omissions, and readers cannot reproduce the analysis.

This skill solves the problem by generating methods text directly from the provenance chain. When every processing step has been logged (via data-provenance), the methods section writes itself. When metadata has been captured from ENCODE (via track-experiments), the experimental details are already recorded. This skill assembles these records into publication-ready prose, figure legends, supplementary tables, and data availability statements.

This standard is not aspirational -- it is the minimum bar for reproducible science.

Scientific Documentation Standards -- Required Metadata

Every methods section MUST report the following fields. Omitting any of these fields produces an incomplete methods section that reviewers will flag and readers cannot reproduce.

Field	Example	Why Required
Library preparation	TruSeq ChIP	Affects fragment size distribution and GC bias
Biological replicates	n=2 per condition	Statistical power and reproducibility
Cells/nuclei per replicate	50,000 cells	Input sufficiency for the assay
Sequencing reads	30M paired-end	Coverage depth determines sensitivity
Read length	2x150 bp	Alignment accuracy and mappability
Paired/single-end	Paired-end	Fragment size estimation, structural variants
Sequencer	NovaSeq 6000	Quality profile, error model, binning
Lab/batch	Snyder Lab, Stanford	Batch effect awareness
Reference genome	GRCh38/hg38	Coordinate system for all downstream analysis
Gene annotation	GENCODE v44	Gene definitions change between versions
ENCODE accessions	ENCSR133RZO	Exact data provenance for reproducibility
Blacklist version	ENCODE Blacklist v2	Artifact exclusion affects all peak-based analyses

How to Populate These Fields

# Track the experiment to capture metadata
encode_track_experiment(accession="ENCSR...", fetch_publications=True)

# Get full experiment details
encode_get_experiment(accession="ENCSR...")

# Get file-level metadata
encode_get_file_info(accession="ENCFF...")

# Get provenance for derived files
encode_get_provenance(file_path="/path/to/derived/file.bed")

Methods Section Templates

Each template below is a fill-in-the-blank paragraph that reads like a real methods section. Bracketed fields [like this] are populated from ENCODE metadata and provenance records. Every template follows these documentation standards.

ChIP-seq Methods

Chromatin immunoprecipitation followed by sequencing (ChIP-seq) data for
[target] in [biosample] were obtained from the ENCODE Project (ENCODE
Project Consortium 2020) under accession [ENCSR accession]. [Library
preparation method] libraries were prepared from [number] biological
replicates ([cells/nuclei] per replicate) and sequenced on an Illumina
[sequencer model] to generate [read count]M [paired-end/single-end]
reads of [read length] bp per replicate.

Raw reads were assessed with FastQC (v[version]; Andrews 2010) and
trimmed with Trim Galore (v[version]; Krueger 2015) to remove adapter
sequences and low-quality bases (Phred < 20). Trimmed reads were aligned
to the [organism] reference genome ([assembly]) using BWA-MEM (v[version];
Li 2013) with default parameters. Duplicate reads were marked and removed
using Picard MarkDuplicates (v[version]; Broad Institute). Reads with
mapping quality < 30 were excluded using samtools (v[version]; Danecek
et al. 2021). Reads mapping to ENCODE Blacklist v2 regions (Amemiya et al.
2019) were removed using bedtools intersect (v[version]; Quinlan & Hall
2010).

Peaks were called using MACS2 (v[version]; Zhang et al. 2008) with
parameters [--broad for broad marks / -q 0.05 for narrow marks]. For
narrow-peak targets, IDR analysis (Li et al. 2011) was performed on
replicate peak sets with a threshold of [0.05]. Signal tracks (fold
change over control) were generated using MACS2 bdgcmp and converted
to bigWig format using bedGraphToBigWig (Kent et al. 2010). Of [N]
called peaks, [N] ([%]) passed IDR filtering and [N] ([%]) remained
after blacklist removal.

ATAC-seq Methods

Assay for Transposase-Accessible Chromatin with sequencing (ATAC-seq) data
for [biosample] were obtained from the ENCODE Project (ENCODE Project
Consortium 2020) under accession [ENCSR accession]. [Number] biological
replicates of [cells/nuclei] [cells/nuclei] each were transposed with
Tn5 transposase ([library kit]) and sequenced on an Illumina [sequencer]
to generate [read count]M [paired-end/single-end] reads of [read length]
bp per replicate.

Raw reads were assessed with FastQC (v[version]; Andrews 2010) and
adapter-trimmed with Trim Galore (v[version]; Krueger 2015). Trimmed
reads were aligned to [assembly] using Bowtie2 (v[version]; Langmead &
Salzberg 2012) with parameters --very-sensitive -X 2000 --no-mixed
--no-discordant. Mitochondrial reads were removed. Duplicate reads were
removed using Picard MarkDuplicates (v[version]; Broad Institute). Reads
with mapping quality < 30 were excluded. Tn5 transposase offset
correction was applied (+4 bp on the positive strand, -5 bp on the
negative strand; Buenrostro et al. 2013). ENCODE Blacklist v2 regions
(Amemiya et al. 2019) were excluded.

Peaks were called using MACS2 (v[version]; Zhang et al. 2008) with
parameters --nomodel --shift -75 --extsize 150 --keep-dup all -q 0.05.
Nucleosome-free fragments (< 150 bp) were used for peak calling. Signal
tracks were generated as fold change over background. TSS enrichment
score was [value] (threshold >= 6; ENCODE data standards; Yan et al. 2020). Of [N]
called peaks, [N] ([%]) passed quality filtering.

RNA-seq Methods

RNA sequencing (RNA-seq) data for [biosample] were obtained from the
ENCODE Project (ENCODE Project Consortium 2020) under accession [ENCSR
accession]. Total RNA was extracted from [number] biological replicates
and [library preparation method] libraries were prepared. Libraries were
sequenced on an Illumina [sequencer] to generate [read count]M
[paired-end/single-end] reads of [read length] bp per replicate.

Raw reads were assessed with FastQC (v[version]; Andrews 2010) and
MultiQC (v[version]; Ewels et al. 2016). Adapter sequences were trimmed
with Trim Galore (v[version]; Krueger 2015). Reads were aligned to
[assembly] with [GENCODE annotation version] gene annotations using STAR
(v[version]; Dobin et al. 2013) in two-pass mode. Gene-level
quantification was performed using RSEM (v[version]; Li & Dewey 2011)
for expected counts and TPM values. Transcript-level quantification was
obtained with Kallisto (v[version]; Bray et al. 2016). Mapping rate was
[%] and rRNA contamination was [%] (thresholds: mapping 70-90%, rRNA <
10%; Conesa et al. 2016; ENCODE data standards). Replicate Pearson correlation was [r value]
(threshold >= 0.9).

Differential expression analysis was performed using DESeq2 (v[version];
Love et al. 2014) in R (v[version]). The design formula was [~ batch +
condition]. Genes with |log2 fold change| > [threshold] and Benjamini-
Hochberg adjusted p-value < [threshold] were considered differentially
expressed. Of [N] genes tested, [N] ([%]) were significantly
upregulated and [N] ([%]) were significantly downregulated.

WGBS Methods

Whole-genome bisulfite sequencing (WGBS) data for [biosample] were
obtained from the ENCODE Project (ENCODE Project Consortium 2020) under
accession [ENCSR accession]. Genomic DNA from [number] biological
replicates was bisulfite-converted using [conversion kit] and sequenced
on an Illumina [sequencer] to generate [read count]M [paired-end/
single-end] reads of [read length] bp per replicate.

Raw reads were assessed with FastQC (v[version]; Andrews 2010) and
trimmed with Trim Galore (v[version]; Krueger 2015) in --rrbs mode
[if applicable] or standard mode. Trimmed reads were aligned to
[assembly] using Bismark (v[version]; Krueger & Andrews 2011) with
Bowtie2 as the underlying aligner. Duplicate reads were removed using
Bismark deduplicate_bismark. Methylation calls were extracted using
MethylDackel (v[version]; https://github.com/dpryan79/MethylDackel)
with parameters [--minDepth X --mergeContext]. Bisulfite conversion
efficiency was [%] (threshold >= 98%; ENCODE data standards). CpG sites with
coverage >= [X]x were retained for analysis ([N] CpGs, [%] of all
genomic CpGs).

Differentially methylated regions (DMRs) were identified using DMRcate
(v[version]; Peters et al. 2021) with parameters [lambda, C]. CpG
sites within ENCODE Blacklist v2 regions (Amemiya et al. 2019) were
excluded prior to DMR calling.

Hi-C Methods

Hi-C data for [biosample] were obtained from the ENCODE Project (ENCODE
Project Consortium 2020) under accession [ENCSR accession]. Chromatin
was crosslinked, digested with [restriction enzyme], proximity-ligated,
and sequenced on an Illumina [sequencer] to generate [read count]M
[paired-end] reads of [read length] bp across [number] biological
replicates.

Read pairs were aligned to [assembly] using BWA-MEM (v[version]; Li
2013) and processed with pairtools (v[version]; Open2C et al. 2024) for
pair classification, deduplication, and filtering. Valid interaction
pairs (cis-reads with MAPQ >= 30) were retained. Contact matrices were
generated using cooler (v[version]; Abdennur & Mirny 2020) at
resolutions of [1kb, 5kb, 10kb, 25kb, 50kb, 100kb, 250kb, 500kb, 1Mb].
Matrices were balanced using iterative correction (ICE; Imakaev et al.
2012). Cis/trans ratio was [%] (threshold > 60%) and long-range cis
interactions (> 20kb) comprised [%] of total (threshold > 40%; Yardimci
et al. 2019).

Chromatin loops were called using HiCCUPS (Juicer v[version]; Durand
et al. 2016) at [resolutions] with default FDR thresholds.
Topologically associating domains (TADs) were identified using [method]
at [resolution]. [N] loops and [N] TADs were identified.

CUT&RUN Methods

Cleavage Under Targets and Release Using Nuclease (CUT&RUN) data for
[target] in [biosample] were obtained from the ENCODE Project (ENCODE
Project Consortium 2020) under accession [ENCSR accession]. [Number]
biological replicates of [cells/nuclei] each were incubated with
[antibody] and processed with protein A-MNase. Released DNA fragments
were sequenced on an Illumina [sequencer] to generate [read count]M
[paired-end] reads of [read length] bp per replicate.

Raw reads were assessed with FastQC (v[version]; Andrews 2010) and
trimmed with Trim Galore (v[version]; Krueger 2015). Reads were aligned
to [assembly] using Bowtie2 (v[version]; Langmead & Salzberg 2012) with
parameters --very-sensitive --no-mixed --no-discordant -X 700. Spike-in
reads (E. coli or S. cerevisiae) were aligned separately and used for
calibration normalization ([spike-in method]). Duplicate reads were
removed with Picard MarkDuplicates (v[version]; Broad Institute). ENCODE
Blacklist v2 and suspect list regions (Amemiya et al. 2019; Nordin et
al. 2023) were excluded.

Peaks were called using SEACR (v[version]; Meers et al. 2019) in
[stringent/relaxed] mode with [IgG control / numeric threshold]. Note:
CUT&RUN has different QC profiles than ChIP-seq; standard ChIP-seq
metrics (FRiP, NSC, RSC) should not be directly applied (Nordin et al.
2023). [N] peaks were identified across replicates.

Single-cell RNA-seq Methods

Single-cell RNA sequencing (scRNA-seq) data for [biosample] were obtained
from the ENCODE Project (ENCODE Project Consortium 2020) under accession
[ENCSR accession]. [Cell/nuclei isolation method] was performed on
[number] biological replicates. Libraries were prepared using [10x
Genomics Chromium / Smart-seq2 / other platform] targeting [N] cells per
replicate, and sequenced on an Illumina [sequencer] to generate [read
count]M reads per sample.

Raw reads were processed with [CellRanger v[version] (10x Genomics) /
STARsolo (STAR v[version]; Dobin et al. 2013)] to generate gene-cell
count matrices against [assembly] with [GENCODE annotation version].
Ambient RNA contamination was estimated and removed using [SoupX
(v[version]; Young & Behjati 2020) / CellBender (v[version]; Fleming
et al. 2023)]. Doublets were detected and removed using [Scrublet
(v[version]; Wolock et al. 2019) / DoubletFinder (v[version]; McGinnis
et al. 2019)]. Cells with [< N genes, > N% mitochondrial reads, < N
UMIs] were excluded.

Data were normalized using [SCTransform (v[version]; Hafemeister & Satija
2019) / scran (v[version]; Lun et al. 2016)] and integrated across
samples using [Harmony (v[version]; Korsunsky et al. 2019) / scVI
(v[version]; Lopez et al. 2018)]. Dimensionality reduction was performed
using PCA ([N] components) followed by UMAP (n_neighbors=[N],
min_dist=[value]; McInnes et al. 2018). Clustering was performed using
the Leiden algorithm (resolution=[value]; Traag et al. 2019) implemented
in [Seurat v[version] (Hao et al. 2024) / Scanpy v[version] (Wolf et
al. 2018)]. Cell types were annotated based on [marker genes / reference
mapping / automated annotation with method]. Final dataset comprised [N]
cells across [N] clusters representing [N] cell types.

Figure Legend Templates

Genome Browser Screenshot

Figure [N]. Genome browser view of [mark/signal] at the [gene name]
locus ([chr]:[start]-[end], [assembly]). Tracks shown from top to
bottom: [track 1 description, e.g., "H3K27ac ChIP-seq signal in
pancreatic islets (ENCSR...)"], [track 2], [track 3]. Signal tracks
represent fold change over input control. [Gene models from GENCODE
v[version] are shown at bottom.] [Shaded region highlights the
[promoter/enhancer/regulatory element] of interest.] Data were obtained
from the ENCODE Project (ENCODE Project Consortium 2020).

Heatmap

Figure [N]. Heatmap of [signal type, e.g., "H3K27ac ChIP-seq signal"]
across [N] [regions/genes/peaks] in [N] [samples/conditions/cell types].
Rows represent [individual peaks / genes / genomic regions], sorted by
[k-means clustering (k=[N]) / hierarchical clustering (method=[ward.D2],
distance=[euclidean]) / signal intensity]. Columns represent [samples /
cell types / conditions]. Color scale indicates [log2 fold change /
z-scored signal / CPM / RPKM] with [color scheme, e.g., "blue-white-red
diverging scale, range [-2, 2]"]. [Normalization method: quantile /
library size / spike-in.] Heatmap was generated using [deepTools
computeMatrix + plotHeatmap (v[version]; Ramirez et al. 2016) /
ComplexHeatmap (v[version]; Gu et al. 2016) / pheatmap (v[version];
Kolde 2019)].

Volcano Plot

Figure [N]. Volcano plot of differentially expressed genes between
[condition A] and [condition B] in [biosample/cell type]. X-axis shows
log2 fold change; y-axis shows -log10 adjusted p-value (Benjamini-
Hochberg correction). Significance thresholds: |log2FC| > [threshold]
(vertical dashed lines) and adjusted p-value < [threshold] (horizontal
dashed line). Red points: [N] significantly upregulated genes. Blue
points: [N] significantly downregulated genes. Gray points: [N] non-
significant genes. [Selected genes are labeled.] Differential expression
analysis was performed using [DESeq2 / edgeR / limma-voom] (v[version]).
[N] total genes were tested.

Hi-C Contact Map

Figure [N]. Hi-C contact frequency map of [region, e.g., "chromosome 7:
25-30 Mb"] in [biosample] at [resolution, e.g., "10 kb"] resolution.
Upper triangle: [observed / Knight-Ruiz balanced / ICE-normalized]
contact frequencies displayed on a [log / linear] color scale. [Lower
triangle: [comparison condition / O/E ratio / difference map].] [Loops
identified by HiCCUPS are marked with [circles/squares].] [TAD
boundaries identified by [method] are shown as [lines/triangles].]
Contact matrix was generated using [cooler (Abdennur & Mirny 2020) /
Juicer (Durand et al. 2016)] from [N]M valid read pairs. Data were
obtained from ENCODE accession [ENCSR...].

UMAP / tSNE

Figure [N]. [UMAP / tSNE] projection of [N] single cells from
[biosample] colored by [cluster identity / gene expression / sample
origin / cell cycle phase]. [UMAP parameters: n_neighbors=[N],
min_dist=[value], computed on [N] principal components. / tSNE
parameters: perplexity=[N], computed on [N] principal components.]
[N] clusters were identified using the [Leiden / Louvain] algorithm
(resolution=[value]). Cell types were annotated based on [canonical
marker gene expression / reference-based mapping using [SingleR /
Azimuth / scArches] / manual curation]. [Inset: expression of [gene
name] across clusters, showing enrichment in [cell type].]

Venn / UpSet Diagram

Figure [N]. [Venn diagram / UpSet plot] showing overlap of [peak sets /
gene lists / regulatory elements] across [N] [conditions / cell types /
datasets]. Set definitions: [Set A] = [N] [peaks/genes] from [source],
[Set B] = [N] [peaks/genes] from [source][, Set C = ...]. Overlaps were
computed using [bedtools intersect (minimum [N] bp overlap) / exact gene
ID matching]. [N] [peaks/genes] ([%]) were shared across all sets; [N]
([%]) were unique to [Set A]. [Statistical significance of overlap was
assessed using [Fisher exact test / hypergeometric test / permutation
test (N=[iterations])]; p = [value].] [UpSet plot generated with
UpSetR (v[version]; Conway et al. 2017).]

Peak Genomic Distribution Bar Chart

Figure [N]. Genomic feature distribution of [N] [peak type, e.g.,
"H3K27ac peaks"] in [biosample]. Peaks were annotated using ChIPseeker
(v[version]; Yu et al. 2015) with [GENCODE v[version] / TxDb annotation
version] gene models. Bars show percentage of peaks overlapping each
genomic feature category: promoter (<= [N] kb from TSS, [%]), 5' UTR
([%]), 3' UTR ([%]), exon ([%]), intron ([%]), downstream (<= [N] kb,
[%]), and intergenic ([%]). [Background genomic distribution is shown
for comparison (gray bars).] [A second panel shows the distribution of
peak distance to nearest TSS.]

Supplementary Table Templates

Table S1: Experiment Metadata

A master table listing all ENCODE experiments used in the study.

Column	Description	Source
ENCODE Accession	ENCSR identifier	encode_get_experiment
Assay	Assay type (ChIP-seq, ATAC-seq, etc.)	encode_get_experiment
Target	Antibody target (for ChIP/CUT&RUN)	encode_get_experiment
Biosample	Tissue or cell type	encode_get_experiment
Organism	Species	encode_get_experiment
Assembly	Genome build	encode_get_experiment
Lab	Submitting laboratory	encode_get_experiment
Biological Replicates	Number of bio reps	encode_get_experiment
Sequencing Reads	Total reads per replicate	File metadata
Read Length	bp, paired/single	File metadata
Sequencer	Instrument model	encode_get_experiment
Library	Library preparation method	encode_get_experiment
Date Released	Release date on ENCODE portal	encode_get_experiment
Audit Status	ENCODE quality audit level	encode_get_experiment

Generate with:

encode_export_data(format="csv")

Table S2: Quality Control Metrics

Assay-specific QC metrics for all experiments. Columns depend on the assay type.

ChIP-seq columns: Accession, Target, Total reads, Mapped reads, Mapping rate, Duplicate rate, NRF, PBC1, PBC2, NSC, RSC, FRiP, Peak count, IDR peaks

ATAC-seq columns: Accession, Total reads, Mapped reads, Mapping rate, Mitochondrial rate, Duplicate rate, NRF, TSS enrichment, Fragment size distribution (NFR/mono/di), FRiP, Peak count

RNA-seq columns: Accession, Total reads, Mapped reads, Mapping rate, rRNA rate, Exonic rate, Intronic rate, Intergenic rate, Genes detected (TPM > 1), Replicate correlation (Pearson r)

WGBS columns: Accession, Total reads, Mapped reads, Mapping rate, Bisulfite conversion rate, CpG coverage (mean), CpGs at 1x, CpGs at 5x, CpGs at 10x, Global methylation (%)

Hi-C columns: Accession, Total read pairs, Valid pairs, Cis pairs (%), Trans pairs (%), Long-range cis (%), Short-range cis (%), Duplicate rate, Library complexity

QC thresholds (include as table footnotes):

• ChIP-seq: FRiP >= 1%, NSC > 1.05, RSC > 0.8, NRF >= 0.8 (Landt et al. 2012)
• ATAC-seq: TSS enrichment >= 5 (GRCh38), >= 6 (hg19), >= 10 (mm10); nucleosomal ladder visible (ENCODE data standards; Yan et al. 2020)
• RNA-seq: Mapping 70-90%, rRNA < 10%, replicate r >= 0.9 (Conesa et al. 2016; ENCODE data standards)
• WGBS: Conversion >= 98%, CpG coverage >= 10x for DMRs (ENCODE data standards)
• Hi-C: Cis/trans > 60%, long-range cis > 40% (Yardimci et al. 2019)

Table S3: Differential Results

For differential expression, differential accessibility, or differential methylation analyses.

Column	Description
Feature ID	Gene ID / Peak ID / CpG ID
Gene Symbol	Gene name (for expression)
Chromosome	chr
Start	Start coordinate
End	End coordinate
log2 Fold Change	Effect size
Standard Error	SE of log2FC
Stat	Test statistic
P-value	Raw p-value
Adjusted P-value	BH-corrected FDR
Base Mean / Mean Signal	Average expression or signal

Table S4: Peak / Region Catalog

For aggregated or filtered peak sets.

Column	Description
Chromosome	chr
Start	Start coordinate (0-based)
End	End coordinate
Name	Peak identifier
Score	Scaled score (0-1000)
Signal Value	Fold enrichment or signal
P-value	-log10(p-value)
Q-value	-log10(q-value)
Summit	Distance from start to summit
Source Accession	ENCODE file accession
Biosample	Source tissue/cell type

Table S5: Software Versions

Every tool used in the analysis with version and citation.

Column	Description
Software	Tool name
Version	Exact version string
Purpose	What it was used for
Citation	Publication reference
DOI	Digital Object Identifier
URL	Download or documentation URL

Data Availability Statements

Standard ENCODE Data Availability

All sequencing data used in this study are publicly available through
the ENCODE Project portal (https://www.encodeproject.org) under the
following experiment accessions: [ENCSR list]. Processed files including
[peak calls / gene quantifications / contact matrices / methylation
calls] used in this analysis are available under file accessions [ENCFF
list]. ENCODE data are released under unrestricted use policy with no
embargo period (ENCODE Project Consortium 2020).

Derived Data Availability