nf-core/methylong

Edit

Introduction

nf-core/methylong is a bioinformatics pipeline that is tailored for long-read methylation calling. This pipeline requires a genome reference as input, and can take either modification-basecalled ONT reads, PacBio HiFi reads (modBam), raw sequencing Pod5 reads or raw Bam reads. The ONT workflow includes modcalling (optional), preprocessing (trim and repair) of reads, genome alignment and methylation calling. The PacBio HiFi workflow includes modcalling (optional), genome alignment and methylation calling. Methylation calls are extracted into BED/BEDGRAPH format, readily for direct downstream analysis. The downstream workflow includes SNV calling, phasing and DMR analysis.

ONT workflow:

1. modcalling (optional)

   • basecall pod5 reads to modBam - dorado basecaller
   • optional: m6A call - fibertools add-nucleosomes

2. trim and repair tags of input modBam

   • trim and repair workflow:
     1. sort modBam - samtools sort
     2. convert modBam to fastq - samtools fastq
     3. trim barcode and adapters - porechop
     4. convert trimmed modfastq to modBam - samtools import
     5. repair MM/ML tags of trimmed modBam - modkit repair

3. align to reference (plus sorting and indexing) - dorado aligner( default) / minimap2

   • optional: remove previous alignment information before running dorado aligner using samtools reset
   • include alignment summary - samtools flagstat

4. create bedMethyl - modkit pileup, 5x base coverage minimum.

   • optional: extract m6A information into bedMethyl - fibertools extract

5. create bedgraphs (optional)

PacBio workflow:

1. modcalling (optional)

   • modcall bam reads to modBam - jasmine (default) or ccsmeth
   • optional: m6A call - fibertools predict-m6a

2. align to reference - pbmm2 (default) or minimap2

   • minimap workflow:

     1. convert modBam to fastq - samtools convert
     2. alignment - minimap2
     3. sort and index - samtools sort
     4. alignment summary - samtools flagstat

   • pbmm2 workflow:

     1. alignment and sorting - pbmm2
     2. index - samtools index
     3. alignment summary - samtools flagstat

3. create bedMethyl - pb-CpG-tools (default) or modkit pileup

   • notes about using pb-CpG-tools pileup:

     • 5x base coverage minimum.
     • 2 pile up methods available from pb-CpG-tools:
       1. default using model
       2. or count (differences described here: https://github.com/PacificBiosciences/pb-CpG-tools)
     • pb-CpG-tools by default merge mC signals on CpG into forward strand. To ‘force’ strand specific signal output, I followed the suggestion mentioned in this issue (PacificBiosciences/pb-CpG-tools#37) which uses HP tags to tag forward and reverse reads, so they were output separately.

   • optional: extract m6A information into bedMethyl - fibertools extract

4. create bedgraph (optional)

Downstream workflow:

1. SNV calling - clair3

2. phasing - whatshap phase

3. DMR analysis

   • includes DMR haplotype level and population scale:

     1. tag reads by haplotype - whatshap haplotype
     2. create bedMethyl - modkit pileup
     3. DMR - DSS (default) or modkit dmr

Usage

Required input:

• ONT or PacBio HiFi reads
  • unaligned modification basecalled bam (modBam)
  • if input modBam was aligned, remove previous alignment information using --reset
  • raw ONT pod5
  • raw bam
• reference genome

First, prepare a samplesheet with your input data that looks as follows:

group,sample,path,ref,method
test,Col_0,ont_modbam.bam,Col_0.fasta,ont
 

Now, you can run the pipeline using:

nextflow run nf-core/methylong \
   -profile <docker/singularity/.../institute> \
   --input samplesheet.csv \
   --outdir <OUTDIR>

For more details and further functionality, please refer to the usage documentation and the parameter documentation.

Pipeline output

To see the results of an example test run with a full size dataset refer to the results tab on the nf-core website pipeline page. For more details about the output files and reports, please refer to the output documentation.

Folder stuctures of the outputs:

 
├── ont/sampleName
│   │
│   ├── fastqc
│   │
│   ├── basecall
│   │   └── calls.bam
│   │
│   ├── trim
│   │   ├── trimmed.fastq.gz
│   │   └── trimmed.log
│   │
│   ├── repair
│   │   ├── repaired.bam
│   │   └── repaired.log
│   │
│   ├── alignment
│   │   ├── aligned.bam
│   │   ├── aligned.bai
│   │   └── aligned.flagstat
│   │
│   ├── snvcall
│   │   ├── merge_output.vcf.gz
│   │   ├── merge_output.vcf.gz.tbi
│   │   └── SNV_PASS.vcf
│   │
│   ├── phase
│   │   ├── phased.vcf
│   │   ├── haplotagged.bam
│   │   └── haplotagged.readlist
│   │
│   ├── pileup/modkit
│   │   ├── pileup.bed.gz
│   │   └── pileup.log
│   │
│   ├── bedgraph
│   │   └── bedgraphs
│   │
│   ├── dmr_haplotype_level/dss
│   │   ├── preprocessed_<1|2|etc>.bed
│   │   ├── DSS_DMLtest.txt
│   │   ├── DSS_callDML.txt
│   │   ├── DSS_callDMR.txt
│   │   └── DSS.log
│   │
│   └── dmr_population_scale
│       ├── population_scale_DMLtest.txt
│       ├── population_scale_callDML.txt
│       ├── population_scale_callDMR.txt
│       └── population_scale.log
│
│
├── pacbio/sampleName
│   │
│   ├── fastqc
│   │
│   ├── modcall
│   │   ├── modbam.bam
│   │   ├── m6a_predicted.bam
│   │   └── m6a.bed
│   │
│   ├── aligned_minimap2/ aligned_pbmm2
│   │   ├── aligned.bam
│   │   ├── aligned.bai/csi
│   │   └── aligned.flagstat
│   │
│   ├── pileup: modkit/pb_cpg_tools
│   │   ├── pileup.bed.gz
│   │   ├── pileup.log
│   │   └── pileup.bw (only pb_cpg_tools)
│   │
│   ├── snvcall
│   │   ├── merge_output.vcf.gz
│   │   └── SNV_PASS.vcf
│   │
│   ├── phase
│   │   ├── phased.vcf
│   │   ├── haplotagged.bam
│   │   └── haplotagged.readlist
│   │
│   ├── bedgraph
│   │   └── bedgraphs
│   │
│   ├── dmr_haplotype_level/dss
│   │   ├── preprocessed_1.bed
│   │   ├── preprocessed_2.bed
│   │   ├── DSS_DMLtest.txt
│   │   ├── DSS_callDML.txt
│   │   ├── DSS_callDMR.txt
│   │   └── DSS.log
│   │
│   └── dmr_population_scale
│       ├── population_scale_DMLtest.txt
│       ├── population_scale_callDML.txt
│       ├── population_scale_callDMR.txt
│       └── population_scale.log
│
│
└── multiqc
    │
    ├── fastqc
    └── flagstat
 

bedgraph outputs all have min. 5x base coverage.

Credits

nf-core/methylong was originally written by Jin Yan Khoo, from the Faculty of Biology of the Ludwig-Maximilians University (LMU) in Munich, Germany, further contributions were made by Yi Jin Xiong, from Central South University (CSU) in Changsha, China.

I thank the following people for their extensive assistance in the development of this pipeline:

• Felix Lenner
• Júlia Mir Pedrol
• Matthias Hörtenhuber
• Sateesh Peri
• Niklas Schandry

Contributions and Support

If you would like to contribute to this pipeline, please see the contributing guidelines.

For further information or help, don’t hesitate to get in touch on the Slack #methylong channel (you can join with this invite).

Citations

If you use nf-core/methylong for your analysis, please cite it using the following doi: 10.5281/zenodo.15366448

An extensive list of references for the tools used by the pipeline can be found in the CITATIONS.md file.

You can cite the nf-core publication as follows:

The nf-core framework for community-curated bioinformatics pipelines.

Philip Ewels, Alexander Peltzer, Sven Fillinger, Harshil Patel, Johannes Alneberg, Andreas Wilm, Maxime Ulysse Garcia, Paolo Di Tommaso & Sven Nahnsen.

Nat Biotechnol. 2020 Feb 13. doi: 10.1038/s41587-020-0439-x.