- Updates
- Project Set-Up
- 1.
FASTQ_Me2 - 2.
CpG_Me2 - 3.
DMRichR - 4.
Comethyl - Interpretting Outputs
- Acknowledgements
- Citation
The original CpG_Me, written by Dr. Ben Laufer, was intended to merge WGBS sequencing lanes and generate cytosine reports to be used in further analysis, such as in DMRichR and Comethyl. Upon trying the pipeline for the first time, I ran into some complications specific to my data (which was not downloaded from SLIMS), and I had the idea to make it more modular and a little easier to use for first-time users. Here are some of the updates I made:
FASTQ_Me2has been rewritten as a python script- Locally downloaded sequencing data can now be used in
FASTQ_Me2 - Users will only interact with command line prompts instead of adapting any of the scripts
CpG_Me2can be run locally (i.e. on screen or at your terminal) as a script instead of only through SLURMCpG_Me2has been written as a snakemake file, yielding the following advantages:- Multi-threading of jobs can be handled locally or on SLURM
- Jobs remove corrupted intermediate files if they fail
- Snakemake's "memory" prevents re-running of samples and files that have already been generated
Clone the repository using your project name in the directory you plan to host the project:
git clone https://github.com/vhaghani26/epigenerator {project_name}
Enter the directory
cd {project_name}
Run the following command in your project directory. It will clone the conda environment with all dependencies needed in order to run the workflow outlined here. This creates an environment called epigenerator. If you would like to change the name, feel free to do so where the command says epigenerator. Please note that this may take quite a few minutes to run.
conda env create -f 00_software/environment.yml --name epigenerator
Activate your environment using
conda activate epigenerator
Run everything downstream of this point in this conda environment. Note that you must activate this environment every time you restart your terminal.
For LaSalle Lab
The environment is already installed and ready in a shared space, so all you will need to do is run:
conda activate /quobyte/lasallegrp/programs/.conda/epigenerator
1. 01_genomes/ Setup
Several steps require genomes for alignment. One of the steps, FastQ-Screen, allows you to align your reads to multiple genomes to check for sources of contamination. In your project directory, create a subdirectory called 01_genomes/:
mkdir 01_genomes/
Within 01_genomes/, create subdirectories corresponding to each genome of interest. The directory structure should look something like the following:
project_directory/
01_genomes/
hg38/
Lamba/
mm10/
PhiX/
rheMac10/
rn6/
Each subdirectory containing your genome of interest should contain the appropriate genome files as described by the FastQ-Screen documentation. Please do not worry about downloading or installing FastQ-Screen. The environment you created has all the software you will need. Activate the environment you created in the previous section, download the reference genome, and index the reference genome using Bowtie2. If you have already processed and indexed a genome(s), then you can instead alias those files/directories into 01_genomes as shown in the "For LaSalle Lab" section below.
2. fastq_screen.conf Setup
In the directory 00_software, locate and open the file fastq_screen.conf. Change lines 34+ to reflect the genomes and location of the genomes that you are aligning to. Provided paths can be absolute or relative paths.
For LaSalle Lab
Create 01_genomes/:
mkdir 01_genomes
Change into the directory:
cd 01_genomes
Run the following commands to link our standard genome directories into this directory:
ln -s /quobyte/lasallegrp/genomes/hg19 .
ln -s /quobyte/lasallegrp/genomes/hg38 .
ln -s /quobyte/lasallegrp/genomes/Lambda .
ln -s /quobyte/lasallegrp/genomes/mm10 .
ln -s /quobyte/lasallegrp/genomes/PhiX .
ln -s /quobyte/lasallegrp/genomes/rheMac10 .
ln -s /quobyte/lasallegrp/genomes/rn6 .
You do not need to update fastq_screen.conf.
If you have more than one lane of data, the sequencing lanes will need to be merged. FASTQ_Me2 handles the merging of lanes whether the data is maintained on SLIMS or already downloaded on your system (i.e. local).
SLIMS Data
If your data is on SLIMS, you will be prompted for your SLIMS string and SLIMS directory. This triggers the data download in your project directory. It also computes the MD5 checksum and compares it to the expected checksum hashes to ensure that the files were not corrupted during download.
Local Data
If you've already downloaded your data, you will be prompted to give the absolute path to the raw data.
task_samples.yaml
FASTQ_Me2 no longer merges lanes, as the pipeline has been modified to align, then merge. This increases alignment efficiency and allows for lower resource requests. FASTQ_Me2 will generate a configuration file, task_samples.yaml, that gets read into CpG_Me2. This task_samples.yaml file is required before running CpG_Me2. In order to generate this file, it will ask for your confirmation along the way to make sure that your samples are being handled correctly. If you already have local data instead of data from SLIMS, it is recommended to put all your raw sequence data in the same directory and name it 01_raw_sequences/. Note that you can use symbolic links within this directory to link to the original files as well. If you would like to skip the FASTQ_Me2 step, you can alternatively make your own configuration file. It must be named task_samples.yaml and be formatted as follows:
---
sequence_data: /quobyte/lasallegrp/Viki/project_name/01_raw_sequences
genome: mm10
samples:
- FA114_FKDN220207669-1A_H5GY3DSX3_L3
- FA114_FKDN220207669-1A_H5GY3DSX3_L2
- FA114_FKDN220207669-1A_H5JF2DSX3_L1
The samples list can go on for as many as you need, but they should include the full file name EXCEPT the read orientation and file extension (i.e. _1.fq.gz, _2.fq.gz, _1.fastq.gz, _2.fastq.gz, _R1.fq.gz, _R2.fq.gz, _R1.fastq.gz, _R2.fastq.gz, _R1_001.fq.gz, _R2_001.fq.gz, _R1_001.fastq.gz, _R2_001.fastq.gz).
Running FASTQ_Me2
To run FASTQ_Me2, all you need to do is run the following in your project directory:
python3 01_FASTQ_Me2.py
CpG_Me2 carries out a number of steps to process your sample, which can be seen in the figure below.
Note that 02_CpG_Me2_PE is for use with paired-end data and 02_CpG_Me2_SE is for use with single-end data. Please substitute the script names as needed in the below commands. I have gone with the deafult of paired-end data in this tutorial since paired-end data is more commonly used than single-end data in our lab.
There are two ways that you can run CpG_Me2.
Since all the information setup was done when you run FASTQ_Me2, all you have to do now is run the following to run CpG_Me2 from start to finish:
snakemake -j 1 -p -s 02_CpG_Me2_PE
The -j option for jobs means how many jobs are able to be run in parallel. The most resource-intensive step uses ~50-100 GB of RAM per sample, so be mindful of the resources you have available if you choose to run more than one job at a time. If, for some reason, you are disconnected from your terminal or your job fails, re-run the above command and it will pick up where it left off. It even deletes possibly corrupted files from where it was cut off to ensure ALL outputs are properly generated.
Notice some issues here.
- Because the alignment step uses so many resources, you are limited by the number of jobs you can run at the same time due to required resource allocations
- If you get disconnected from your terminal, your job fails. This means that you need to be logged in and have an active terminal for days at a time. This can be circumvented by running
CpG_Me2inscreen, but you will still need the proper resource allocation for that length of time.
For the above reasons, it is HIGHLY recommended that you instead run it via SLURM. This is similar to what was initially written by Ben, where individual jobs get submitted to SLURM.
In the directory 00_slurm/, there is a file named config.yaml. You will need to modify two things:
- Update your SLURM partition for the two lines (line 18 and line 19) containing
--accountand--partitionby inputting a string. This will look something like--partition=production - Change the
conda_prefix(line 33). It should look something like/software/anaconda3/4.8.3/lssc0-linux/,/home/vhaghani/anaconda3/, or/quobyte/lasallegrp/programs/.conda/
Once you have updated config.yaml, go back to your project directory. Snakemake manages the submission of jobs, so wherever you run it, it will need to stay open. As such, I recommend running it in screen. Activate the conda environment (confirm you are in the environment if you are using screen). When you are ready, run:
snakemake -s 02_CpG_Me2_PE --profile 00_slurm/
For LaSalle Lab
There are some bugs with the way your home directory is mounted on the Genome Center cluster. This will result in Snakemake erroring out when it tries to cache files. To fix it, we will need to make some changes. Briefly, you are creating a cache directory within your project directory. Since this is in a shared lab space, it is properly mounted for use. You are giving full permissions to the cache directory to ensure that Snakemake can store files in the directory. Then, you are setting the Snakemake cache variable to the proper directory so Snakemake knows to use that directory for cache. After these changes are made, you can submit to SLURM. Please run the following in your project directory (ideally in screen with your conda environment activated):
cd 00_slurm/
mkdir .cache
chmod -R 777 .cache
cd ..
Change the paths to whatever path you are using and put the following in your .bashrc or .profile (or whatever you're using):
export SNAKEMAKE_OUTPUT_CACHE=/quobyte/lasallegrp/{your_directory}/{your_project}/00_slurm/.cache
export XDG_CACHE_HOME=/quobyte/lasallegrp/{your_directory}/{your_project}/00_slurm/.cache
Note that these will have to be updated every time you run CpG_Me2, but it is currently the only workaround I'm aware of for dealing with the issues the Genome Center cluster has.
Then, run the following:
snakemake -s 02_CpG_Me2_PE --profile 00_slurm/ --cache 00_slurm/.cache
I am hoping to one day better incorporate DMRichR here, but in the meantime. please refer to Dr. Ben Laufer's DMRichR package for further analysis. This package will aid in the identification of differentially methylated regions within your dataset.
If you are interested in running DMRichR on SLURM, there is a template SLURM script in 00_slurm/ titled DMRichR_SLURM_template.sh. If you are using it, please ensure that you change the appropriate variables, including the #SBATCH parameter section and the path to the DM.R script.
If you are interested in weighted region comethylation network analysis, please refer to Dr. Charles Mordaunt's Comethyl package for further analysis.
The most important outputs you will find are 01_raw_sequences/ and 08_cytosine_reports/. You can feel free to delete the intermediate files, but they are included for your reference in case something goes wrong when you try to run CpG_Me2. The 03_screened and 09_multiqc directories may also be good to keep so you can check on the qualities of your samples at various stages. Below is a description of each of the output directories and what is contained in each. Note that the numbers preceding the directory names are reflective of the order they were generated.
This directory contains the standard error per job per sample. If, for some reason, you encounter an error for a sample, this allows you to view the full error message and debug.
This directory documents the time it took each job to run. This could help inform future decisions about time allocations in SLURM should you need to adjust the time for a submitted job.
This directory contains the primary genome you want to align your data to as well as whatever genomes you want to screen your samples for using FastQ-Screen.
This directory contains the merged sequences for each sample, separated into forward ({sample}_1.fq.gz) and reverse ({sample}_2.fq.gz) reads.
This contains the reads after they have been trimmed, meaning that the adapters and low quality base pairs at sequence ends have been removed.
All samples are screened against some different background genomes to determine sample read origins. This could be helpful if you are concerned about possible contamination of your samples.
This directory contains the aligned sequences, meaning that samples are mapped against your reference genome.
This contains your reads with PCR duplicates removed.
This contains your aligned and deduplicated sequence reads sorted in chromosome coordinate order.
Nucleotide coverage of reads are calculated for input BAM files and output here.
This contains the various outputs of Bismark associated with extracting methylation information from each sample.
According to the documentation, this "provides useful metrics for validating library construction including the insert size distribution and read orientation of paired-end libraries."
This directory contains the cytosine reports that are required for downstream analysis.
This contains quality read outs for the eligible files maintained in the other directories. It can be helpful to view if any of your samples are being problematic.
This directory is generated when you use SLURM to run CpG_Me2. It contains the err and out files typically generated by SLURM per job submission, which corresponds to individual rules run per sample.
This work was largely adapted from Dr. Ben Laufer's original CpG_Me Program. These updates could not have been made without the help of Jules Mouat and Aron Mendiola.
If you are using epigenerator, please cite the following (noting that it may be subject to change upon further publication):
-
Laufer BI*, Neier KE*, Valenzuela AE, Yasui DH, Lein PJ, LaSalle JM. Placenta and Fetal Brain Share a Neurodevelopmental Disorder DNA Methylation Profile in a Mouse Model of Prenatal PCB Exposure. Cell Reports, 2022. doi: 10.1016/j.celrep.2022.110442
-
Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics, 2011. doi: 10.1093/bioinformatics/btr167
-
Matrin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal, 2011. doi: 10.14806/ej.17.1.200
-
Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 2016. doi: 10.1093/bioinformatics/btw354
-
Haghani V, Mouat JS. 2023. Epigenerator. GitHub. https://github.com/vhaghani26/epigenerator.
If you are using Comethyl, also make sure to cite the following:
- Mordaunt CE, Mouat JS, Schmidt RJ, and LaSalle JM. (2022) Comethyl: a network-based methylome approach to investigate the multivariate nature of health and disease. Briefings in Bioinformatics bbab554.