Peak-calling Algorithms (WonderPeaks and PeakStream) as Tools for Improved ChIP-seq and Transcriptomic Analysis in Fungal Pathogens

Summary

This report introduces WonderPeaks, a novel computational tool for analyzing RNA-seq and ChIP-seq data. This tool successfully identifies peaks (read pileups) in sequencing data, enabling the characterization of untranslated region boundaries in RNA-seq and detecting chromatin enrichment in ChIP-seq, providing valuable insights for fungal pathogen research.

Abstract

Characterizing changes in gene expression through transcriptomics and transcription regulator activity has become a fundamental approach for understanding the diverse responses involved in fungal pathogenesis. This paper introduces two computational tools designed to address key challenges in the study of transcriptional regulation in fungal pathogens, particularly non-model with limited genomic annotation. First, we present WonderPeaks, a novel peak-calling algorithm that leverages the first derivative of mapped genomic data from next-generation sequencing (NGS) experiments to identify enriched peaks in Chromatin ImmunoPrecipitation followed by sequencing (ChIP-seq). Second, we introduce PeakStream, an extension of WonderPeaks for annotating 3' untranslated regions (UTRs) in transcriptomic data generated using poly(A)-primed library preparation. Together, these tools provide an end-to-end data analysis pipeline, offering a user-friendly solution for researchers studying transcription regulation in fungi. We demonstrate their effectiveness with data from the fungal pathogen Candida albicans, successfully identifying verified peaks in ChIP-seq data and annotating validated UTRs through comparison with total RNA sequencing data under the same conditions. We also discuss the limitations of WonderPeaks for ChIP-seq data compared to current state-of-the-art methods and propose directions for future improvements. Ultimately, this work provides practical guidance and powerful resources for studying transcriptional regulation, with immediate relevance to pathogenic fungi and potential applications in broader genomic studies.

Introduction

Fungal pathogens are an emerging global health concern, with infections rising in recent years¹. Many of these pathogens exhibit high antifungal resistance and are associated with significant mortality rates². However, compared to model fungal organisms, many pathogenic fungi remain poorly characterized, highlighting the need for further research into their mechanisms of pathogenicity. Next-generation sequencing (NGS) techniques like Chromatin Immunoprecipitation Sequencing (ChIP-Seq) and RNA-Sequencing (RNA-Seq) play a critical role in uncovering the molecular mechanisms of gene expression underlying fungal pathogenicity.

The quality of insights derived from NGS data is highly dependent on the accuracy of the software used for raw data analysis. Among the key challenges to NGS data analysis is peak calling -- accurately identifying regions of enriched NGS reads -- which is particularly complex due to the wide variety of library preparation and sequencing techniques, making a universal solution impractical. The MACS algorithm³, including its more recent version, MACS3, is widely regarded as the gold standard for analyzing ChIP-seq datasets. However, MACS relies on user-defined parameters -- such as minimum peak length and maximum gap -- which may not be universally applicable and are often difficult to determine prior to analysis. Notably, the latest MACS3 release includes a cut-off analysis feature that allows users to estimate parameters before peak calling. For enhanced performance, users can also provide a list of "blacklisted" genomic regions known to introduce bias due to chromatin structure or copy number variation. While MACS remains the most commonly used and trusted peak-calling tool for ChIP-seq data, few alternative algorithms are available, particularly for cases requiring highly tailored parameter settings.

RNA-Seq is an invaluable technique for studying the gene expression responses of pathogenic fungi during in vivo growth, such as in tissue culture or mouse infection models^4,5,6,7. A high sequencing depth is required for accurate differential expression analysis in these conditions, which can be cost- and resource-prohibitive^8,9. Library preparation methods like poly-adenylation (poly(A))-priming Sequencing (3'RNA-Seq), which uses primers designed to anneal to the poly(A)-tails of mRNA for cDNA generation, can help reduce the sequencing depth needed for gene expression analysis¹⁰. However, this approach relies on high-quality genome annotations, particularly of 3' Untranslated Regions (UTRs), where peaks from poly(A)-priming events are typically located¹¹. The genome annotations of many understudied fungal pathogens lack any UTR annotations, making the use of 3'RNA-Seq in these organisms difficult. In addition, the UTR length for a single gene can be dynamic across different growth conditions and cell types^12,13. While a number of new analysis tools have been developed to identify and annotate UTRs, many of them are designed for mammalian datasets, whose gene organization differs greatly from that of fungi, or require data from independent sequencing experiments, such as single-cell or reverse mRNA sequencing, that can increase time and costs for a researcher looking to do transcriptome analysis^12,14,15.

In this paper, we present WonderPeaks, a novel peak-calling software, designed on the principles of the first derivative, that can be used to dynamically call peaks in NGS datasets (Figure 1). WonderPeaks identifies peaks by calculating the first derivative of the coverage signal and using this value -- the slope of the peak -- to define potential peaks. The algorithm searches for instances where the first derivative exhibits a local maximum above a user-supplied or data-inferred slope threshold (indicating an increasing signal), followed by a local minimum above the same threshold (indicating a decreasing signal), thereby detecting all candidate peaks in the dataset. For ChIP-seq applications, WonderPeaks compares all candidate peaks between test and control samples to identify uniquely enriched peaks. By applying WonderPeaks to a previously published ChIP-seq dataset of a transcription factor in the fungal pathogen Candida albicans¹⁶, we demonstrated its ability to successfully identify peaks upstream of key genes highlighted in the original study, while also discussing the algorithm's current limitations in this application.

We also introduce PeakStream, a software tool that leverages WonderPeaks to identify peaks in 3'RNA-Seq datasets. 3'RNA-Seq libraries depend on accurate 3' UTR annotations, as reads generated through poly(A)-priming often extend beyond the stop codon of the coding sequences (CDS) of genes and thus are uncounted when using standard annotations focused solely on coding regions. The PeakStream analysis pipeline was designed to create new genome annotations using 3' RNA-Seq data, focusing on regions downstream of gene coding sequence (CDS) regions. PeakStream assigns these peaks to genes, generating a new genome annotation for use in downstream read counting programs. We show that the use of PeakStream can accurately identify and assign downstream poly(A)-generated peaks to the proper gene in a 3'RNA-Seq C. albicans dataset. PeakStream also annotates peaks unlikely to be associated with any current gene annotations, facilitating the discovery of possible novel transcripts. Together, PeakStream and WonderPeaks represent a powerful suite of user-friendly tools for peak detection in next-generation sequencing (NGS) datasets.

Figure 1: Overview figure of peak-calling by WonderPeaks and PeakStream. Left: Peaking calling using the first derivative. Top Right: Peak-calling on ChIP-Seq datasets using WonderPeaks. Bottom Right: Peak-calling on RNA-Seq datasets using PeakStream. Please click here to view a larger version of this figure.

Protocol

1. Installation (skip if completed)

1. Prerequisites
   1. Install Anaconda or Minconda to load requirements to run preprocessing (section 5) and WonderPeaks (section 6 or section 7).
      NOTE: For a user guide for Anaconda, see the reference¹⁶.
   2. Install Python: Python within Anaconda or Miniconda.
   3. Install Jupyter Notebooks to execute all functions in this method.
      NOTE: A beginner's user guide for Jupyter Notebooks can be found in the reference¹⁷.
      CAUTION: For fungal genomes (≤100 Mbp), ensure a computing environment with at least 20 cores, 8 GB of RAM, and 30 GB of available disk space.
2. Install preprocessing functions.
   1. In the terminal, execute: conda create -n WP_preprocessing
   2. In the terminal, execute: conda activate WP_preprocessing
   3. In the terminal, execute: conda env update --file environment.yml --name
      NOTE: The environment.yml is a file that contains all the package dependencies and should be downloaded from https://github.com/mgarber21/WonderPeaks_preprocessing.git.
   4. In the terminal, execute: pip install WonderPeaks-preprocessing
   5. In the terminal, execute: conda deactivate WP_preprocessing
3. Install WonderPeaks Functions:
   1. In the terminal, execute: conda create -n WonderPeaks
   2. In the terminal, execute: conda activate WonderPeaks
   3. In the terminal, execute: conda env update --file environment.yml --name
      NOTE: The environment.yml is a file that contains all the package dependencies and should be downloaded from https://github.com/mgarber21/WonderPeaks.git.
   4. In the terminal, execute: pip install WonderPeaks
   5. In the terminal, execute: conda deactivate WonderPeaks
      NOTE: Steps 1.2 and 1.3 accomplish the following: They create a dedicated Conda environment for preprocessing (section 5) and WonderPeaks (sections 6 and 7), isolating dependencies to avoid conflicts with other software. They activate the environment, setting it up for installing and running WP_preprocessing or WonderPeaks-specific functions. They install the software dependencies and tools required for preprocessing the data. They deactivate the environment when not in use to prevent accidental modification and free up system resources.
4. Download Jupyter Notebooks and templates from the WonderPeaks GitHub repository. Upload the WonderPeaks downloads to the operating system containing the raw data (directory will be created in section 2).
   NOTE: The Jupyter Notebooks contain prewritten scripts and templates necessary for running preprocessing, WonderPeaks, and PeakStream workflows. Uploading them to the same system as the raw data ensures the paths and directories align correctly.

2. Create a data directory

NOTE: WonderPeaks and PeakStream workflows require that all data (raw and processed) is stored in the same directory. This step explains how to create this new directory ({data_directory} = /path/to/your/data) and how to move experimental raw data (unprocessed sequencing reads) to a folder within this directory called raw_data.

1. Create a data directory.
   1. In the terminal, execute mkdir {data_directory} (e.g. mkdir /path/to/your/data)
2. Create a raw data subdirectory for unprocessed sequencing reads.
   1. In the terminal, execute mkdir {data_directory}/raw_data (e.g. mkdir /path/to/your/data /raw_data)
3. Move the unprocessed sequencing reads to the raw data directory.
   1. In the terminal, execute mv {current_path_to_raw_data}/*fastq* {data_directory}/raw_data (e.g. mv current/data/path/*fastq* /path/to/your/data /raw_data

3. Create the user inputs file (NGS_user_input.csv)

NOTE: The input file specifies the user-generated configurations for running preprocessing and WonderPeaks.

1. Download template NGS_user_inputs.csv from the WonderPeaks GitHub repository.
2. Update fields in NGS_user_inputs.csv. Update the fields as follows:
   Data directory: /path/to/your/data
   Genome directory: /path/to/your/genome
   Genome fasta: genome.fasta
   Genome annotation: genome_annotation.gtf (provide the file name of the genome annotation in GTF format)
3. Save updated NGS_user_inputs.csv to the data directory created in section 2.
   CAUTION: Do not change the file name; WonderPeaks will only recognize this file if it is named NGS_user_input.csv.

4. Create metadata file (NGS_user_metadata.csv)

NOTE: The metadata file is used to store any information relevant to the experiment. Additional columns can be added to describe the conditions of the experiment, but these will not affect subsequent steps.

1. Download template NGS_user_metadata.csv from the WonderPeaks GitHub repository.
2. Update fields in NGS_user_metadata.csv The fields are as follows:
   1. file: Ensure that the file name has no spaces, includes the file handle (e.g., fastq, fastq.gz), and does not include the absolute path.
   2. bedgraph: Specify whether the file should be included in PeakStream by setting this field to either TRUE or FALSE.
      1. Set the bedgraph field to FALSE when the file can reasonably be excluded from PeakStream analysis. For example, in an RNAseq experiment, set the bedgraph parameter bedgraph=FALSE for mutants or other cases where UTR differences between samples are not expected. However, be sure to set the bedgraph field to TRUE for all control files in an RNAseq experiment and for all files in a ChIPseq experiment.
   3. designfactor
      1. designfactor1: Specify a design factor relevant to the experimental design (e.g., treatment or sample_type).
      2. designfactor2: Specify a second design factor relevant to the experimental design (e.g., strain or epitope). For an RNAseq experiment, include treatment and strain as typical design factors. The treatment column lists the treatments applied (e.g., control, drug1), and the strain column lists strain information (e.g., wildtype, mutant). For a ChIPseq experiment, include sample_type and epitope as typical design factors. The sample_type column lists whether the protein was tagged or untagged, and the epitope column lists the name of the epitope used.
         NOTE: designfactors are attributes specific to the experimental design.
         WonderPeaks is compatible with either untagged controls or input controls as the baseline.
      3. Make sure that design factor columns do not include unique replicate number (e.g., sample_type: [tagged, tagged, untagged_control, untagged_control] not sample_type: [tagged_1, tagged_2, untagged_control _1, untagged_control _2]. Providing unique replicate numbers will result in an error during the run.
      4. Be sure to use underscores (_) instead of spaces in design factor names.
      5. For ChIPseq application, ensure that the sample_type (or custom name) column in the metadata file includes terms that contain words tag and control.
         NOTE: For example, valid entries could include tagged and untagged_control. Providing a sample_type column without these terms will result in an error during the run.
         CAUTION: For ChIP application, the user must specify two design factors.
      6. Add design factors to the relevant row of NGS_user_inputs.csv. Ensure that the column names used for design factors are listed as a semicolon separated string (e.g., treatment;strain or sample_type;epitope).
         CAUTION: The design factors in NGS_user_inputs.csv must match columns of NGS_user_metadata.csv exactly. Any mismatched run will result in an error during the run (Figure 2, Supplemental Table S1, and Supplemental Table S2).

Figure 2: Example NGS_user_input.csv and NGS_user_metadata.csv. Examples of the NGS_user_input.csv (top panel) and NGS_user_metadata.csv (bottom panel), highlighting the match between the designfactor and designfactor columns with pink or blue text and arrows. Please click here to view a larger version of this figure.

5. Preprocessing NGS data

NOTE: Skip to section 5 or section 6 if using custom preprocessing.

1. Open NGS preprocessing Jupyter Notebook (NGS_Preprocessing.ipynb).
2. Activate WP_preprocessing environment (created in step 1.2) in the upper-right corner of the notebook interface.
3. Execute the first cell by holding down the Shift key and then pressing Enter (Shift+Enter).
4. In the second cell of the Jupyter Notebook, update the directory path by setting the Directory = "path/to/your/data", where path/to/your/data/ is the directory created in section 2.
5. Generate alignment files. The preprocessing functions will perform Trimming using FastP¹⁸; Quality Control using FastQC¹⁹ and MultiQC²⁰; Alignment using STAR²¹. The output alignment files are saved to a subdirectory named startout within the data directory (e.g., path/to/your/data/starout); Filtering (optional) using samtools view²², filter the alignment file to only keep reads above a threshold specified in NGS_user_inputs.csv.
   NOTE: These functions will process only single reads (e.g., R1) of the dataset at a time. Users can specify run options for FastP and STAR in the NGS_user_inputs.csv (e.g., FastP: adapter_sequence (optional); STAR: genomeDir, genomeFastaFiles, sjdbGTFfil).
6. Execute preprocessing functions in the second cell with Shift + Enter.
   NOTE: Tasks in the second cell may take several hours to complete. If the run is interrupted, repeat steps 5.3-5.6 to restart the run. Progress from previous steps will not be overwritten, and the process will continue where it left off.
7. Generate trace files of the alignment coverage using BamCoverage²³ (see steps 5.7.1 and 5.7.2).
   NOTE: For ChIPseq, WonderPeaks requires single bedgraph files containing coverage for both forward and reverse reads. For RNAseq with Poly(A) priming, PeakStream requires two bedgraph files, one for forward reads (_fwd.bedgraph) and one for reverse reads (_rev.bedgraph). Forward and reverse reads are generated using the filterRNAstrand parameter within BamCoverage²³.
   1. ChIPseq using the following parameters: outfilfeformat="bedgraph", strand=None, binsize=20, smoothLength=60, minMappingQuality=255, normalizeUsing="CPM".
      NOTE: Output: Produces single bedgraph files containing coverage for both forward and reverse reads. Output is stored in /path/to/your/data/bedgraphout (Figure 3).
      1. Execute the BamCoverage function in the third cell using Shift + Enter.
   2. RNAseq using the following parameters: outfilfeformat="bedgraph", strand="forward" or "reverse", binsize=20, smoothLength=60, minMappingQuality=255, normalizeUsing="CPM".
      CAUTION: Be sure to execute function twice with strand set to either forward or reverse to generate files for reads in both directions.
      NOTE: Output: Produces two bedgraph files: one for forward reads (_fwd.bedgraph) and one for reverse reads (_rev.bedgraph). Output is stored in /path/to/your/data/ bedgraphout (Figure 3).
      1. Execute the BamCoverage function in the third cell using Shift + Enter.

Figure 3: File organization for WonderPeaks. A screenshot of the data folder with the bedgraph files in the directory bedgrapghout/normalizeUsingCPM. Please click here to view a larger version of this figure.

6. WonderPeaks for ChIPseq

1. Precheck
   1. Confirm all bedgraph files with filehandle .bedgraph are in a subdirectory within the data directory called bedgraphout (Figure 3).
   2. Confirm that designfactors in user_inputs file (NGS_user_inputs.csv) (Figure 2) match columns in the metadata file (NGS_user_metadata.csv) and that the rows of the designfactor columns are not unique (see caution in step 4.2.4).
2. Open NGS preprocessing Jupyter Notebook (WP4ChIP.ipynb).
3. Activate WonderPeaks (environment created in step 1.3) in the upper-right corner of the notebook interface.
4. Execute the cells using Shift+Enter until the breakpoint to run peak calling. Once finished, a record of the data processed will be saved and stored in a subdirectory within the data directory called WonderPeaks
   1. Look for WOnder_init.csv: a concatenation of all raw coverage and the results from calculating the first derivative.
   2. Note WOnder_unfiltered_peaks.csv: a concatenation of all unfiltered peaks called based on the first derivative.
   3. Observe bedgraph_summary.csv: a summary of score statistics after grouping each file and chromosome.
5. Define run parameters:
   1. Run the 1^st cell below the markdown breakpoint.
      NOTE: A plot displaying the raw data segregated by the specified designfactors and a table showing the design factors will pop up; use the table and plot to determine values in subsequent steps (Figure 4).
   2. In the next cell, specify the values for score_cut, fold_change, and designfactor (Figure 4).
      1. score_cut is the threshold value used to determine whether a peak is to be considered in the output. To determine the score_cut, observe the plot and choose a value close to the median of the tagged data (see hashed line, Figure 4). Enter this value as follows: score_cut= value.
      2. fold_change is the threshold value of the ratio scores of tagged:untagged used to determine whether a peak is considered real. To determine the fold_change, observe the plot and choose a value above the ratio of the medians of the untagged and tagged data. Enter this value as follows: fold_change= value.
      3. designfactor_value is specified as part of the experimental design. Possible design factors are listed in red in the printed table. To determine the designfactor, choose one of the values listed in red. Enter this value between quotations as follows: designfactor_value ="{ value}".
6. Execute the next cells using Shift + Enter to run peak filtering and mapping. The data and summary plots will be stored in a subdirectory within the data directory called WonderPeaks.
   1. Observe {designfactor_value}_taggedVuntagged.csv: A pivot table of all overlapping peaks with column for each of the tagged and untagged samples.
   2. Note {designfactor_value}_all_tagged_peaks.csv: A summary table of all the real peaks, based on user parameters (step 6.5).
   3. Observe {designfactor_value}_peaks2gtf.csv: A mapping of the real peaks, based on user parameters (step 6.5), to genes in the user specified annotations file.
7. Optional: Toggle the parameters in step 6.5 by re-executing steps 6.5-6.6. If using the same designfactor_value, the files generated, as described in step 6.6, will be overwritten.

Figure 4: Screenshot highlighting the designfactor_value and user-specified thresholds in WonderPeaks for ChIP-seq. Screenshot of WonderPeaks jupyter notebook, highlighting the possible designfactor_value options from the displayed table, and how to implement the designfactor_value in the next cell. The top black arrow points to a table displaying possible designfactor_value entries; the Op value is circled and shown as the selected user input for designfactor_value in the options cell (bottom black arrow). In the graph, solid and dashed lines indicate the approximate median peak scores for tagged and untagged samples, respectively, in the opaque cell experiments. These medians are used to define the score_cut (tagged median) and fold_change (the ratio of tagged to untagged medians) parameters. Please click here to view a larger version of this figure.

7. PeakStream for 3'RNAseq

1. Precheck:
   1. Confirm all bedgraph files are in a subdirectory within the data directory called bedgraphout (Figure 2).
   2. Open NGS preprocessing Jupyter Notebook (PeakStream.ipynb)
   3. Activate WonderPeaks (environment created in step 1.3) in the upper-right corner of the notebook interface.
2. Execute the cells using Shift+Enter until the breakpoint to run peak calling and peak mapping. Once finished, a new annotation file with predicted 3' UTRs and FeatureCounts²⁴ read count files will be saved and stored in a subdirectory within your data directory called PeakStream (Figure 5) .
   ​NOTE: By default, the output file will only include annotations for protein coding biotypes, but this can be toggled using the biotype option.

Figure 5: File organization for PeakStream. A screenshot of the data folder with the bedgraph files in the directory bedgrapghout. Please click here to view a larger version of this figure.

Results

WonderPeaks
After performing a ChIP-seq experiment, researchers commonly use peak callers, such as MACS³, to identify genomic regions enriched by an epitope-tagged DNA-binding protein. We developed WonderPeaks as a user-friendly peak caller designed to identify peaks using the method described above.

WonderPeaks identifies peaks by first calculating the first derivative of the coverage and uses this value (the slope of the peak) to define potential ...

Discussion

Next-generation sequencing (NGS) techniques provide unparalleled insight into gene regulation and expression in fungal pathogens. As such, computational tools must be both comprehensive -- capturing all the data generated in an experiment -- and accessible to general users, particularly bench scientists. In this report, we introduced two tools, WonderPeaks and PeakStream, that address these needs for fungal pathogen researchers in ChIP-seq and RNA-seq workflows.

WonderPeaks is a user-friendly ...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
CORALL Total RNA-seq V1 kit	Lexogen	095	Wet lab material
Filamentous-Fungi riboPOOL	siTOOLs	dp-P096-6	Wet lab material
High Sensitivity RNA ScreenTape	Agilent	5067-5579	Wet lab material
High Sensitivity RNA ScreenTape Ladder	Agilent	5067-5581	Wet lab material
High Sensitivity RNA ScreenTape Sample Buffer	Agilent	5067-5580	Wet lab material
https://github.com/mgarber21/WonderPeaks/blob/main/environment.yml			list of dependecies for WonderPeaks
https://github.com/mgarber21/WonderPeaks_preprocessing/blob/main/environment.yml			list of dependecies for WonderPeaks_preprocessing
Monarch Spin RNA Cleanup Ki	NEB	T2040L	Wet lab material
pygenometracks (3.9)
QuantSeq 3′ mRNA-Seq FWD Library Prep Kit V1	Lexogen	015	Wet lab material
Qubit RNA High Sensitivity (HS) Assay Kit	Invitrogen	Q32852	Wet lab material
RNA Clean & Concentrator-5	Zymo Research	R1016	Wet lab material
TURBO DNA-free Kit	ThermoFisher	AM1907	Wet lab material
WonderPeaks(0.1.14)
WonderPeaks_preprocessing(0.2.3)