Aligns and QCs methylated-RNA-immunoprecipitation (MeRIP / m6A-seq) IP and input libraries using STAR or HISAT2 splice-aware mapping, samtools sort/index, IP/input matched-pair tracking, antibody-lot metadata recording, replicate concordance via deepTools multiBamSummary + plotCorrelation, IP enrichment QC via plotFingerprint and per-transcript IP/input ratio distributions, library-complexity saturation curves via PreSeq, and the explicit do-NOT-deduplicate convention for standard non-UMI MeRIP. Use when preparing paired IP and input BAM files for exomePeak2 / MeTPeak / MACS3 peak calling, evaluating MeRIP replicate concordance and IP enrichment, deciding whether to deduplicate (standard MeRIP typically NOT), choosing genome-vs-transcriptome alignment for downstream peak vs m6Anet workflows, recording antibody clone and lot metadata for cross-batch reconciliation, detecting failed IPs via saturation curves and IP/input distribution shape, or generating IP-over-Input bigWig tracks for visualisation.
Instrucciones de origen · Vista previa de solo lectura
name
bio-epitranscriptomics-merip-preprocessing
description
Aligns and QCs methylated-RNA-immunoprecipitation (MeRIP / m6A-seq) IP and input libraries using STAR or HISAT2 splice-aware mapping, samtools sort/index, IP/input matched-pair tracking, antibody-lot metadata recording, replicate concordance via deepTools multiBamSummary + plotCorrelation, IP enrichment QC via plotFingerprint and per-transcript IP/input ratio distributions, library-complexity saturation curves via PreSeq, and the explicit do-NOT-deduplicate convention for standard non-UMI MeRIP. Use when preparing paired IP and input BAM files for exomePeak2 / MeTPeak / MACS3 peak calling, evaluating MeRIP replicate concordance and IP enrichment, deciding whether to deduplicate (standard MeRIP typically NOT), choosing genome-vs-transcriptome alignment for downstream peak vs m6Anet workflows, recording antibody clone and lot metadata for cross-batch reconciliation, detecting failed IPs via saturation curves and IP/input distribution shape, or generating IP-over-Input bigWig tracks for visualisation.
Before using code patterns, verify installed versions match. If versions differ:
CLI: <tool> --version then <tool> --help to confirm flags
Python: pip show <package> then help(module.function)
If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
STAR --outSAMtype accepts BAM SortedByCoordinate since 2.5.x; check the Log.final.out file for input/output statistics. deepTools bamCompare --operation log2 is the modern syntax (older --ratio log2 still works but is being phased out). PreSeq c_curve and lc_extrap have stable interfaces; HISAT2 reports unique vs multi-mapped in the summary log.
MeRIP-seq Preprocessing
"Get my MeRIP IP and input libraries ready for peak calling" -> Trim adapters with MeRIP-appropriate defaults (do NOT trim UMIs unless the library is UMI-MeRIP — most are not), splice-aware-align IP and input to the GENOME (not transcriptome) with STAR / HISAT2, sort and index, evaluate replicate concordance and IP enrichment with deepTools, build a saturation curve per library so peak counts can be honestly compared across libraries, record antibody clone and lot metadata so cross-batch comparison is later auditable, and produce IP-over-Input log2 bigWig tracks for downstream visualisation. Crucially, do NOT deduplicate non-UMI MeRIP — see the failure-modes section.
CLI: STAR --runMode alignReads -- splice-aware genome alignment, the field default
The Single Most Important Modern Insight -- Peak counts are library-size-dependent; saturation curves are the only honest cross-library comparison
A MeRIP library sequenced to 20 million unique reads finds substantially fewer peaks than the same biology at 60 million reads. Per-sample peak counts reported without saturation curves (PreSeq c_curve / lc_extrap; Daley & Smith 2013 Nat Methods 10:325) are uninterpretable across studies and often across replicates within a study. Subsample BAMs to a common unique-read depth before peak calling for any cross-condition peak-count comparison, OR report peaks alongside the saturation curve. A corollary: do NOT deduplicate standard MeRIP — the typical MeRIP protocol (Synaptic Systems 202-003 / Abcam ab151230 / NEB EpiMark E1610 antibody pull-down on fragmented poly(A)-selected RNA) has NO unique molecular identifiers, and picard MarkDuplicates on such libraries collapses real biological replicates of high-coverage transcripts (the opposite of what dedup achieves in DNA ChIP-seq). UMI-MeRIP is the only exception — and most MeRIP libraries in print are NOT UMI. McIntyre et al. 2020 Sci Rep 10:6590 demonstrated that replicate-to-replicate peak overlap is ~80% within a single lab but drops to a median 45% between labs using nominally identical conditions; this irreducible technical noise constrains how strongly any single MeRIP study can support biological claims, and the preprocessing pipeline is where the variance is set.
Algorithmic Taxonomy
Tool / step
Mechanism
Output
Strength
Fails when
STAR 2.7+ (Dobin 2013 Bioinformatics 29:15)
Two-pass splice-aware alignment with on-the-fly splice junction database
Sorted BAM + splice-junction TSV
Field default; multi-mapper retention configurable; STAR splice-junction DB
Memory-heavy (~30 GB human); slower than HISAT2
HISAT2 2.2.1+ (Kim 2019 Nat Biotechnol 37:907)
Hierarchical graph FM-index; splice-aware
Sorted BAM
~5x lighter memory than STAR; comparable accuracy
Less mature splice-junction handling for novel introns
BWA-MEM2 (Vasimuddin 2019 IPDPS 314)
DNA-style local alignment; NO splice awareness
Sorted BAM
Use ONLY for transcriptome-aligned MeRIP (rare)
Splits reads across exon junctions if used on genome
fastp 0.23+ (Chen 2018 Bioinformatics 34:i884)
Streaming adapter detection + quality trim
Trimmed FASTQ + JSON QC
Fast; JSON-readable QC output
UMI handling disabled by default; do NOT pass --umi for standard non-UMI MeRIP (the opposite of the failure direction in some other library types)
Trim Galore
Wrapper over cutadapt with paired-end auto-detect
Trimmed FASTQ
Conservative defaults; widely cited
Slower than fastp on large datasets
samtools sort / index
BAM coordinate sort + .bai index
Sorted BAM + index
Standard
None at default
Picard MarkDuplicates
Identifies PCR duplicates by 5' alignment start
Marked / removed BAM
Standard in DNA / ChIP
The dominant MeRIP convention is to SKIP dedup for non-UMI libraries (collapses real biology at high-coverage transcripts); a minority of pipelines dedup MeRIP — record the choice in metadata
deepTools multiBamSummary + plotCorrelation
Per-bin read counts; Spearman / Pearson matrix
Heatmap + clustering
Standard replicate-concordance plot
Bin size sensitive (use 10 kb for transcriptome-genome)
deepTools plotFingerprint (Diaz 2012 Stat Appl Genet Mol Biol 11:9)
Cumulative read-fraction vs cumulative-bin-fraction Lorenz curve
PDF + raw counts
Direct IP-vs-input enrichment QC; "good" IP has steep tail
Aligning to transcriptome (rare; specific downstream tools)
BWA-MEM2 or bowtie2; defer to read-alignment/
STAR splice-aware on transcriptome causes spurious splice calls inside transcripts
Methodology evolves; before any high-stakes preprocessing pipeline, web-search "STAR vs HISAT2 MeRIP 2024" and "MeRIP saturation curve preseq" for current consensus parameters.
Adapter Trimming for MeRIP
Goal: Remove sequencing adapters and low-quality 3' ends WITHOUT removing biological signal (UMI-MeRIP must keep UMI sequences; standard MeRIP does not have UMIs and trimming should be minimal).
Approach: Use fastp or Trim Galore with adapter auto-detection; require minimum read length 25-30 nt (shorter reads multi-map and confound exomePeak2); for standard non-UMI MeRIP, do NOT pass --umi flags; preserve random-hexamer-priming artifacts ONLY if downstream pipeline expects them (most do not).
For UMI-MeRIP (rare), insert --umi --umi_loc read1 --umi_len 8 BEFORE the alignment step. Default fastp output preserves base quality information needed by downstream variant-aware tools; do NOT pass --disable_quality_filtering for MeRIP libraries.
STAR Splice-Aware Alignment for IP and Input
Goal: Produce coordinate-sorted, indexed GENOME BAM files for each IP and input library with splice-junction-aware mapping, retaining a moderate number of multi-mappers for accurate per-window read counts at multi-isoform loci.
Approach: Build STAR genome index once with the matched GENCODE / Ensembl GTF used downstream; loop IP and input samples with identical parameters; retain up to 20 multi-mappers per read (MeRIP read counts at multi-isoform genes need this); request explicit BAM SortedByCoordinate; emit splice-junction tables for QC.
STAR \
--runMode genomeGenerate \
--genomeDir star_index \
--genomeFastaFiles genome.fa \
--sjdbGTFfile annotation.gtf \
--sjdbOverhang 100 \
--runThreadN 12
mkdir -p aligned
for sample in IP_rep1 IP_rep2 IP_rep3 Input_rep1 Input_rep2 Input_rep3; do
STAR \
--runMode alignReads \
--genomeDir star_index \
--readFilesIn trimmed/${sample}_R1.fq.gz trimmed/${sample}_R2.fq.gz \
--readFilesCommand zcat \
--outSAMtype BAM SortedByCoordinate \
--outFilterMultimapNmax 20 \
--outSAMattributes NH HI AS nM NM MD \
--outFileNamePrefix aligned/${sample}_ \
--runThreadN 12
samtools index -@ 4 aligned/${sample}_Aligned.sortedByCoord.out.bam
done
--outFilterMultimapNmax 20 is intentional: MeRIP at rRNA / snoRNA / pseudogene-rich loci needs multi-mapper retention. Reduce to 1 only if downstream analysis explicitly cannot tolerate multi-mappers. --sjdbOverhang should equal (read length - 1) but 100 is the common-enough default for 100-150 bp reads.
HISAT2 Alternative for Memory-Constrained Environments
Goal: Achieve splice-aware alignment in ~5x less memory than STAR (12-16 GB suffices for human), with comparable accuracy for MeRIP applications.
Approach: Build HISAT2 graph index; align with --dta for downstream-transcript-assembly compatibility; pipe directly to samtools sort.
HISAT2 multi-mapper handling is governed by -k; the default reports the primary alignment only. For MeRIP, pass -k 5 if multi-mapper-aware downstream counting is required.
Per-Sample QC: flagstat and idxstats
mkdir -p qc
for bam in aligned/*sortedByCoord.out.bam aligned/*sorted.bam; do
name=$(basename${bam} .bam)
samtools flagstat ${bam} > qc/${name}.flagstat
samtools idxstats ${bam} > qc/${name}.idxstats
done
Inspect flagstat for properly-paired rate (>=85% indicates good pairing); inspect idxstats for unexpected chromosome-level read piles (rRNA bleed-through, mitochondrial domination — both are MeRIP red flags).
Replicate Concordance via deepTools
Goal: Quantify how similar replicate IP libraries are to each other (and likewise input libraries) using a Spearman correlation matrix; flag a divergent replicate before it propagates into peak calling.
Approach: Compute genome-wide per-bin read counts at 10 kb resolution across all IP and input BAMs; convert to a clustered Spearman heatmap with deepTools plotCorrelation.
IP replicates within a condition should cluster (Spearman >= 0.85 typical); input replicates should cluster with each other; IP and input should NOT cluster together. A failed IP looks like input.
IP Enrichment via plotFingerprint
Goal: Confirm IP libraries are enriched (a few transcripts have many reads) and input libraries are uniform (reads spread across transcripts); fail-fast on poor IP before peak calling.
Approach: deepTools plotFingerprint builds a cumulative Lorenz-style curve; a steep tail = signal concentrated in few regions (good IP); a diagonal = uniform coverage (input or failed IP). The framework is from ChIP-seq (Diaz 2012 Stat Appl Genet Mol Biol 11:9) and transfers cleanly to MeRIP.
Good MeRIP IP: cumulative-fraction-of-reads vs cumulative-fraction-of-bins curve sits well below the diagonal in the right half (top-X% of bins capture >50% of reads). Input: near-diagonal. The --outQualityMetrics file reports JS distance and synthetic JS distance; the IP-vs-Input JS distance is a single-number IP-quality summary (higher = more concentrated signal).
Library Complexity / Saturation Curves via PreSeq
Goal: Compute per-library complexity so peak counts can be honestly compared across libraries and conditions of different sequencing depth.
Approach: PreSeq c_curve (interpolation up to observed depth) and lc_extrap (extrapolation beyond observed) on the sorted BAM. Daley & Smith 2013 Nat Methods 10:325 capture-recapture model.
mkdir -p complexity
for bam in aligned/*.bam; do
name=$(basename${bam} .bam)
preseq c_curve -B -o complexity/${name}_c_curve.txt ${bam}
preseq lc_extrap -B -o complexity/${name}_lc_extrap.txt ${bam}done
Inspect: the lc_extrap curve plots distinct molecules vs total reads; a plateau indicates saturation. For cross-condition peak-count comparison: pick a common depth (often 30M unique reads), subsample with samtools view -s 0.<frac> to that depth, THEN call peaks.
IP-over-Input bigWig for Downstream Visualisation
Goal: Produce a per-bin log2 (IP / Input) coverage track per replicate, ready for downstream metagene / browser plots.
Approach: deepTools bamCompare with --operation log2; choose a sensible pseudocount to avoid divide-by-zero at low-coverage bins.
--pseudocount 1 prevents division-by-zero at zero-coverage bins; --binSize 25 is fine-grained enough to preserve peak topology while keeping bigWig files reasonably sized.
Per-Method Failure Modes
Dedup applied to non-UMI MeRIP
Trigger:picard MarkDuplicates REMOVE_DUPLICATES=true invoked on a standard MeRIP BAM that has no UMI.
Mechanism: Standard MeRIP libraries have no unique molecular identifiers. PCR duplicates and biological re-sampling at high-coverage transcripts look identical at the alignment level. Dedup removes both, collapsing real coverage at the most-abundant transcripts to an artificially flat profile. This is the opposite of dedup's intent in DNA ChIP-seq.
Symptom: Coverage at housekeeping mRNAs (e.g., GAPDH, ACTB) drops 5-20x after dedup; downstream peak counts at highly-expressed transcripts collapse; volcano plot of differential peaks shows expression-driven false positives.
Fix: Skip dedup for standard non-UMI MeRIP. If the library is UMI-MeRIP, use umi_tools dedup (Smith 2017 Genome Res 27:491) which respects UMI rather than alignment position alone. Record dedup status in sample-sheet metadata.
Transcriptome alignment for downstream peak calling
Trigger: STAR or bowtie2 aligned to transcriptome FASTA, then BAM passed to exomePeak2 / MeTPeak / MACS3.
Mechanism: exomePeak2 and MeTPeak expect a GENOME BAM plus GTF; they project peaks back to transcript features internally. A transcriptome BAM has reads in per-transcript coordinates which the GTF cannot resolve back to genome coordinates without re-alignment.
Symptom: exomePeak2 throws errors on TxDb-genome consistency; MeTPeak returns zero peaks; MACS3 calls peaks on transcript IDs as if they were chromosomes.
Fix: Align to GENOME with STAR / HISAT2 for downstream MeRIP peak calling. Transcriptome alignment is correct only for m6anet-analysis (ONT DRS) and rare quantification-only downstream tools.
Failed IP indistinguishable from input
Trigger: A single replicate IP library has IP/input ratio distribution centred at 1.0 across all transcripts (no enrichment); fingerprint Lorenz curve sits at the diagonal.
Mechanism: Failed IP — antibody-RNA binding did not enrich m6A-containing fragments. Causes include antibody-batch defect, insufficient pulldown wash, RNA degradation during IP, or accidental mock IgG IP.
Symptom: plotFingerprint shows IP overlaying input on the Lorenz plot; per-transcript IP/input ratio histogram is centred at 1.0; downstream peak callers find few or no peaks AT THE FAILED REPLICATE while other replicates produce normal counts.
Fix: Identify the failed replicate via plotFingerprint AND IP/input ratio distribution BEFORE peak calling; exclude or re-do. Single failed IP in a 3-replicate design routinely produces "differential" peaks driven entirely by the failure.
Antibody lot mismatch across samples
Trigger: A multi-condition MeRIP study uses Synaptic Systems 202-003 antibody lot A for the control IPs and lot B for the treatment IPs (because lot A ran out mid-study).
Mechanism: Anti-m6A polyclonals (Synaptic Systems 202-003, Abcam ab151230, NEB EpiMark E1610, Cell Signaling 56593, Active Motif 61755) have batch-to-batch variability in pulldown efficiency and m6A-vs-m6Am cross-reactivity. Pooling lot-A and lot-B counts in a downstream differential model attributes lot-effect to condition.
Symptom: "Differential" peaks at high abundance transcripts; effect sizes track antibody lot rather than condition; reanalysis with lot in the design matrix removes most differential peaks.
Fix: Record antibody_clone and antibody_lot per sample in metadata; include lot as a fixed effect in downstream differential analysis. Within a single study, ideally use a single lot for ALL replicates and ALL conditions.
Peak counts compared across libraries of different depth
Trigger: "Condition A has 14,000 peaks; condition B has 22,000 peaks; condition B has more m6A."
Mechanism: Peak count is library-size-dependent. A library at 60M unique reads finds more peaks than 30M. Without rarefaction or saturation correction, peak-count comparisons across libraries are dominated by sequencing depth.
Symptom: Peak counts track total mapped reads more closely than they track biological condition; downstream "biological m6A change" claims do not survive rarefaction-to-common-depth.
Fix: Either rarefy all BAMs to common unique-read depth before peak calling, OR fit saturation curves with PreSeq lc_extrap and compare at matched depth, OR report peak count alongside the saturation curve.
Random hexamer priming over-trim
Trigger: Aggressive 5' trimming of the first 6-12 nt to remove "random hexamer priming bias" applied to MeRIP libraries.
Mechanism: Random hexamer priming bias affects the 5' nucleotide composition of reads but does NOT degrade downstream peak-calling accuracy. Over-trimming removes biological signal and shortens reads enough to inflate multi-mapper fraction.
Fix: Standard adapter trimming with --length_required 25 is sufficient; do not 5'-trim for hexamer bias unless downstream tooling explicitly requires unbiased 5' ends (most do not). The bias is a known artifact in the RNA-seq community and is robust to standard analytical pipelines.
Reconciliation: When QC Signals Disagree
Pattern
Likely cause
Action
plotFingerprint diagonal but IP/input ratio shows enrichment
Mismatched chromosome naming (chr1 vs 1) between samples
Verify `samtools view -H bam
Replicate Spearman 0.95 but plotFingerprint diverges
Replicates correlate in bulk but differ in IP enrichment depth
Check per-sample sequencing depth; reduce to common depth
Saturation curve plateaus early but peak count low
Library too shallow OR genome too small (BAM under 1M unique reads)
Use c_curve only; or sequence deeper
MultiQC misses STAR Log.final.out
STAR output naming non-standard
Re-run with --outFileNamePrefix and rerun MultiQC; check multiqc_config.yaml search patterns
picard MarkDuplicates collapses all reads to 1 per position
Tiny BAM or single read pair per fragment
Verify BAM has many properly-paired reads; do NOT dedup non-UMI MeRIP regardless
Empty fingerprint output
All BAMs have identical bin coverage
Verify BAMs are different files; check multiBamSummary --outRawCounts
bigWig file size too large
Bin size too small at deep coverage
Increase --binSize from 25 to 50; bigWig is lossy at large bin sizes
Saturation curve never plateaus
Library deeply under-sampled
Sequence deeper OR accept curve does not plateau and report accordingly
fastp --umi errors on non-UMI library
UMI flag passed but library has no UMI
Drop --umi flag for standard non-UMI MeRIP
Anticipated Reviewer Pushback
Pushback
Response
"Was deduplication applied?"
No — standard non-UMI MeRIP protocol; PCR duplicate vs biological resampling indistinguishable without UMI; dedup collapses real coverage at high-expression transcripts
"What is the IP enrichment QC?"
deepTools plotFingerprint reported per replicate; JS distance >=0.5 vs input
"Are the replicates concordant?"
Spearman correlation matrix reported via deepTools plotCorrelation on 10 kb bins; IP-IP within condition >=0.85
"Saturation curve?"
PreSeq lc_extrap per library; libraries rarefied to common depth before downstream peak calling
"What antibody clone and lot?"
Recorded per sample in metadata; same lot for all replicates within study
"Why STAR instead of HISAT2?"
STAR splice-junction-DB-based vs HISAT2 graph-based; both valid for MeRIP; choice driven by memory budget
"How many biological replicates?"
N >=3 per condition (per McIntyre 2020); N=2 is under-powered for differential downstream
"Was alignment to genome or transcriptome?"
Genome (required for exomePeak2 / MeTPeak / MACS3 downstream); transcriptome alignment is for m6anet-analysis only
References
Dobin A, Davis CA, Schlesinger F et al (2013) STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29(1):15-21. doi:10.1093/bioinformatics/bts635
Kim D, Paggi JM, Park C, Bennett C, Salzberg SL (2019) Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol 37(8):907-915. doi:10.1038/s41587-019-0201-4
Vasimuddin Md, Misra S, Li H, Aluru S (2019) Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems. IPDPS 314-324. doi:10.1109/IPDPS.2019.00041
Chen S, Zhou Y, Chen Y, Gu J (2018) fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34(17):i884-i890. doi:10.1093/bioinformatics/bty560
Ramírez F, Ryan DP, Grüning B et al (2016) deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res 44(W1):W160-W165. doi:10.1093/nar/gkw257
Diaz A, Park K, Lim DA, Song JS (2012) Normalization, bias correction, and peak calling for ChIP-seq. Stat Appl Genet Mol Biol 11(3):Article 9. doi:10.1515/1544-6115.1750
Daley T, Smith AD (2013) Predicting the molecular complexity of sequencing libraries. Nat Methods 10(4):325-327. doi:10.1038/nmeth.2375
Smith T, Heger A, Sudbery I (2017) UMI-tools: modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy. Genome Res 27(3):491-499. doi:10.1101/gr.209601.116
McIntyre ABR, Gokhale NS, Cerchietti L, Jaffrey SR, Horner SM, Mason CE (2020) Limits in the detection of m6A changes using MeRIP/m6A-seq. Sci Rep 10(1):6590. doi:10.1038/s41598-020-63355-3
Ewels P, Magnusson M, Lundin S, Käller M (2016) MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32(19):3047-3048. doi:10.1093/bioinformatics/btw354
Dominissini D, Moshitch-Moshkovitz S, Schwartz S et al (2012) Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature 485(7397):201-206. doi:10.1038/nature11112
Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR (2012) Comprehensive analysis of mRNA methylation reveals enrichment in 3' UTRs and near stop codons. Cell 149(7):1635-1646. doi:10.1016/j.cell.2012.05.003
Related Skills
m6a-peak-calling - Immediate downstream consumer of the IP/input BAM pairs
m6a-differential - Downstream differential analysis on peak count matrices; design matrix relies on IP/input pairing recorded here
m6anet-analysis - ONT DRS alternative; uses TRANSCRIPTOME alignment with minimap2, NOT the genome BAMs produced here
modification-visualization - Uses the bigWig output of bamCompare for metagene plots and browser tracks
read-qc/quality-reports - FastQC / MultiQC upstream of trimming
read-alignment/star-alignment - General STAR splice-aware alignment patterns
read-alignment/hisat2-alignment - HISAT2 graph-based alternative; general usage
alignment-files/sam-bam-basics - General BAM mechanics, samtools fundamentals
