MBDToolBox

<h3>MBDToolBox</h3><p>MBDToolBox - Docker container to support the data descriptor paper describing an MBD-seq dataset on neuroblastoma samples.</p><h3>References:</h3><ul><li><p>DNA-methylation profiling of primary neuroblastoma tumors using methyl-CpG-binding domain sequencing (MBD-seq)<br>Anneleen Decock, Maté Ongenaert, Wim Van Criekinge, Frank Speleman, Jo Vandesompele</p></li><li><p>Methyl-CpG-binding domain sequencing reveals a prognostic methylation signature in neuroblastoma</p></li></ul><p>Anneleen Decock, Maté Ongenaert, Robrecht Cannoodt, Kimberly Verniers, Bram De Wilde, Geneviève Laureys, Nadine Van Roy, Ana P. Berbegall, Julie Bienertova-Vasku, Nick Bown, Nathalie Clément, Valérie Combaret, Michelle Haber, Claire Hoyoux, Jayne Murray, Rosa Noguera, Gaelle Pierron, Gudrun Schleiermacher, Johannes H. Schulte, Ray L. Stallings, Deborah A. Tweddle for the Children’s Cancer and Leukaemia Group (CCLG), Katleen De Preter, Frank Speleman, Jo Vandesompele.</p><p>Oncotarget, 2015<br>doi: 10.18632/oncotarget.6477<br><a href="http://www.ncbi.nlm.nih.gov/pubmed/26646589">http://www.ncbi.nlm.nih.gov/pubmed/26646589</a></p><h3>Tools:</h3><ul><li>sratoolkit (v2.5.4-1) - NCBI toolkit to work with files from SRA (e.g. convert SRA to FASTQ)</li><li>fastqc (v0.11.4) - QC on raw sequencing reads (FASTQ files)</li><li>bowtie2 (v2.2.6) - maps reads to the reference genome (hg19 index files included as well)</li><li>samtools (v0.1.19) - tools to manipulate SAM files (used to sort/index BAM files)</li><li>picard (v1.140) - tools to manipulate next-generation sequencing formats (BAM files)(used to mark duplicate reads in the BAM file)</li><li>samstat (v1.5.1) - QC statistics on mapped reads (SAM/BAM files)</li><li>BamUtil (v1.0.13) - statistics and operations, QC on mapped reads (BAM files)</li><li>MACS (v1.4.3) - peak caller for ChIP-seq or other capture methods (such as MBD-seq)</li></ul><h3>Usage:</h3><p>This toolbox works with all kinds of enrichment strategies, the example given here is based on a paired-end MBD-seq experiment on a primary neuroblastoma sample</p><h4>Steps:</h4><ol><li>Get the SRA-file (SRA repository, NCBI) and convert it to FASTQ files (paired-end) (SRAtoolkit)</li><li>Perform QC on the FASTQ raw reads (FastQC)</li><li>Map reads to the human hg19 reference genome taking into account the paired-end nature (bowtie2)</li><li>Sort and index the SAM/BAM files (samtools)</li><li>Mark duplicates (PCR duplicates during library prep) (picard)</li><li>Produce statistics and QC on BAM files (samstat, BamUtils)</li><li>Call peaks to identify enriched regions, covered by MBD (MACS)</li></ol><h4>Complete script:</h4><p>1. Get the SRA-file (SRA repository, NCBI) and convert it to FASTQ files (paired-end) (SRAtoolkit)</p><blockquote><p>set some local options</p></blockquote><p><code>sudo locale-gen en_US.UTF-8</code></p><blockquote><p>get the SRA file from SRA archives through FTP</p></blockquote><p><code>curl ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX%2FSRX103%2FSRX1039196/SRR2040898/SRR2040898.sra -o 1509_MBD.sra</code></p><p><code>curl ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX%2FSRX103%2FSRX1039241/SRR2040943/SRR2040943.sra -o 1509_input.sra</code></p><blockquote><p>output: downloaded SRA file</p><p>SRA to FASTQ reads, paired-end (--split-files used) and GZipped (giving two fastq.gz files: one with _1 and the other with _2, representing the paired reads)</p></blockquote><p><code>fastq-dump --split-files 1509_MBD.sra --gzip</code><br><code>fastq-dump --split-files 1509_input.sra --gzip</code></p><blockquote><p>output: 2 fastq files, ending with _1 and _2, containing the paired reads</p></blockquote><p>2. Perform QC on the FASTQ raw reads (FastQC)</p><blockquote><p>make output directory</p></blockquote><p><code>mkdir OUT_fastqc</code></p><blockquote><p>FastQC analysis on one of the pairs, results in the created output directory</p></blockquote><p><code>fastqc 1509_MBD_1.fastq.gz --outdir OUT_fastqc/</code></p><blockquote><p>output: in <em>OUT_fastqc</em> directory: a HTML file and a ZIP file, containing the QC report and the images</p></blockquote><p>3. Map reads to the human hg19 reference genome taking into account the paired-end nature (bowtie2)</p><blockquote><p>read mapping, qualities from SRA are always phred+33 (--phred33), paired-end with max. 500 bp between pairs (-X 500)<br>8 parallel processes (--threads 8), reference genome hg19 (-x hg19; bowtie2 index files included in container)<br>output passed to samtools to create BAM file instead of SAMfile</p></blockquote><p><code>bowtie2 -q --phred33 --sensitive -X 500 --threads 8 -t -x bowtie2/hg19 -1 1509_MBD_1.fastq.gz -2 1509_MBD_2.fastq.gz | samtools view -bhSo 1509_MBD_unsorted.bam -</code></p><p><code>bowtie2 -q --phred33 --sensitive -X 500 --threads 8 -t -x bowtie2/hg19 -1 1509_input_1.fastq.gz -2 1509_input_2.fastq.gz | samtools view -bhSo 1509_input_unsorted.bam -</code></p><blockquote><p>output: unsorted BAM file (Bowtie2 output SAM format, which is directly, on the fly converted to BAM using<em>samtools view</em> command; don't forget the "-" at the end of the command, it is required for this passing (pipe) to samtools to work</p><p>if you want to, you could remove all FASTQ files and the SRA</p></blockquote><p><code>rm 1509*.sra</code></p><p><code>rm 1509*.fastq.gz</code></p><p>4. Sort and index the SAM/BAM files (samtools)</p><blockquote><p>sort BAM file using samtools</p></blockquote><p><code>samtools sort 1509_MBD_unsorted.bam 1509_MBD.sorted</code><br><code>samtools sort 1509_input_unsorted.bam 1509_input.sorted</code></p><blockquote><p>output: sorted BAM file, ready for processing with Picard tools</p><p>if you want to, you could remove the unsorted BAM file</p></blockquote><p><code>rm 1509_MBD_unsorted.bam</code><br><code>rm 1509_input_unsorted.bam</code></p><p>5. Mark duplicates (PCR duplicates during library prep) (picard)</p><blockquote><p>mark duplicate reads</p></blockquote><p><code>java -Xmx8g -jar picard/picard.jar MarkDuplicates INPUT=1509_MBD.sorted.bam OUTPUT=1509_MBD.sorted_nodups.bam ASSUME_SORTED=true METRICS_FILE=Picard_1509_MBD_duplicates.txt</code></p><p><code>java -Xmx8g -jar picard/picard.jar MarkDuplicates INPUT=1509_input.sorted.bam OUTPUT=1509_input.sorted_nodups.bam ASSUME_SORTED=true METRICS_FILE=Picard_1509_input_duplicates.txt</code></p><blockquote><p>output: sorted BAM file, with duplicate reads marked; run summary by Picard (including the read numbers, paired reads, reads marked duplicate)</p><p>if you want to, you could remove the BAM file without duplicates marked</p></blockquote><p><code>rm 1509_MBD.sorted.bam</code><br><code>rm 1509_input.sorted.bam</code></p><blockquote><p>index sorted BAM file with duplicates flagged</p></blockquote><p><code>samtools index 1509_MBD.sorted_nodups.bam</code><br><code>samtools index 1509_input.sorted_nodups.bam</code></p><blockquote><p>output: BAM indes, ending with .bam.bai</p></blockquote><p>6. Produce statistics and QC on BAM files (samstat, BamUtils)</p><blockquote><p>perform samstat on sorted BAM file</p></blockquote><p><code>samstat 1509_MBD.sorted_nodups.bam</code></p><blockquote><p>output: HTML file with samstat quality statistics</p><p>perform bamUtils on sorted BAM file</p></blockquote><p><code>bam stats --in 1509_MBD.sorted_nodups.bam --basic --bamIndex 1509_MBD.sorted_nodups.bam.bai</code></p><p><code>bam stats --in 1509_MBD.sorted_nodups.bam --phred --bamIndex 1509_MBD.sorted_nodups.bam.bai</code></p><blockquote><p>output: BamUtil statistics (outputted on screen)</p></blockquote><p>7. Call peaks to identify enriched regions, covered by MBD (MACS)</p><blockquote><p>make output directory</p></blockquote><p><code>mkdir OUT_macs</code></p><blockquote><p>call peaks and write a single WIG file for visualization in IGV etc. (--wig --single-profile)</p></blockquote><p><code>macs -t 1509_MBD.sorted_nodups.bam -c 1509_input.sorted_nodups.bam --outdir=OUT_macs/ --name=1509_MBD_macspeaks -f BAM --petdist=200 -g hs --wig --single-profile</code></p><blockquote><p>output: MACS peak files (Excel-file and BED files); WIG files for visualization, all in the <em>OUT_macs</em> directory</p></blockquote><h3>Alternative scripts</h3><p>Other included scripts are examples of further analysis.</p><blockquote><p>Example how to run BALM - also 'peak / enriched regions' based</p></blockquote><p><code>script_BALM.sh</code></p><blockquote><p>R scripts and utilities to create a count matrix around transcription start sites (TSS) or in 5 kb genomic windows</p></blockquote><pre><code>script_count2000.R TSS2000.bed </code></pre><blockquote><p>R script demonstrating MEDIPS, an R/BioConductor package, calculating absolute methylation scores for all CG sites</p></blockquote><p><code>script_MEDIPS.R</code></p><p>Alternative scripts</p><p>1. Get the TDF and BED.gz files from the GEO repositories (Cohort I:GSE69224; Cohort II: GSE69243) </p><p>These files are in the "Supplementary file" section of GEO; you can download all these files at the same time. </p><p>For the supplied IGV session example file, you need the files for cohort II, which can be all downloaded from this location: http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE69243&format=file </p><p>Caution: file size: 32 GB</p><p>2. <strong>Decompress the combined archive if needed (if downloading multiple files) and place them in a directory with the supplied IGV session file (example for cohort II, file:"MBD cohort II.xml") and the sample annotations file ("annotations.txt")</strong> </p><p>It is crucial that the files are in the same directory as the session-file as the paths are defined relative</p><p>3. <strong>Download and open IGV (<a href="https://www.broadinstitute.org/igv/">https://www.broadinstitute.org/igv/</a>) and load the session file.</strong> </p><p>Alternatively, individual tracks (TDF files and BED files) can be loaded into IGV.</p><div><br>