Tutorial

Install MethFlowVM

How to install VirtualBox and VirtualBox Extension Pack

Before you can use the virtual machine of MethFlow (MethFlowVM) you need to download and install a hypervisor or a virtual machine monitor, i.e. a computer software that creates and runs virtual machines. There are many hypervisors that are compatible with the Open Virtualization Format (OVF or OVA), like VirtualBox, VMware and XenServer, among others. OVA is the format that we use to distribute MethFlowVM. We recommend using VirtualBox (version 5 or higher) because it is free, open source and user friendly.

To install VirtualBox visit this website.

Download the version for your operating system (Windows, OS X, Linux and Solaris) and the VirtualBox Extension Pack (is the same for all supported platforms). Install VirtualBox and then double-click on the VirtualBox Extension Pack file to install it. If the double click does not work, open VirtualBox and go to

File ➜ Preferences ➜ Extensions

to add the VirtualBox Extension Pack file. Once installed, restart VirtualBox.

↑

Download Virtual Machine

Download the MethFlowVM OVA file here.

Then open VirtualBox and go to

File ➜ Import Appliance

to import the MethFlowVM. Select the OVA file and leave the default options (later we will change these options). Import the OVA file.

↑

Configure Virtual Machine

Once imported, you will see MethFlowVM in the VirtualBox interface (left menu). Now click on MethFlowVM to select it for configuration.

To change the memory allocation, go to

Settings (orange gear) ➜ System ➜ Motherboard ➜ Base Memory

We recommend using 1.5 GB of memory per Gb of the reference assembly. For example, the human genome has roughly 3Gb and therefore the recommended memory would be 4.5.

To change the number of CPUs assigned, go to

Settings (orange gear) ➜ System ➜ Processor ➜ Processor(s)

You must use at least 4 CPUs for proper functioning.

↑

Use a shared folder

It is highly recommended to use a shared folder to store the datasets, assemblies and results. If you do not use a shared folder, the performance of MethFlow will be considerably lower and the virtual disk will grow undesirably. Please use a shared folder.

To select the shared folder in the host machine:

Settings (orange gear) ➜ Shared Folders

In order to ‘activate’ a shared folder, two pieces of information are needed:

The name and path of the folder in the host machine.
The name of the shared folder for the virtual machine.

Note that it is important to activate the auto-mount feature.

Within the virtual machine, the shared folder will be on /media and the name of the folder will be the name you entered, preceded by sf_. For example, /media/sf_mydata. Please note that the path of the shared folder on the host and the virtual machine are different.

↑

Launch the virtual machine

After closing the configuration window, in the VirtualBox interface, select MethFlowVM and click on Start button (green arrow). MethFlowVM will start in a new window. To log into the VM, enter your credentials:

User: methflow
Password: methflow

Optionally, you can change the password, by typing in a terminal:

passwd

Note that the virtual machine can access internet if the host machine is connected to the web. However, the VM is not accessible from outside and therefor a lax password should be not security issue.

↑

Set your working folder

Just after logging in, you will be asked which working folder you want to use. If you have not added any shared folder you will use /home/methflow to store data (strongly not recommended!), else you must choose a shared folder.

At later stages, you can change the working folder to any shared or unshared folder. Open terminal and type set_working_folder .

↑

Change keyboard layout and time zone

To change the keyboard layout, type the following command in a terminal of the MethFlowVM:

sudo dpkg-reconfigure keyboard-configuration

You will be asked the password for your user in the MethFlowVM (default: methflow). Select your keyboard layout and ready.

To change your time zone, type the following command in a terminal of the MethFlowVM:

sudo tzselect

You will be asked the password for your user in the MethFlowVM (default: methflow). Select your time zone and ready.

↑

Quick Start

Get test datasets

The best way to start using MethFlow is with a toy example. To get the test datasets, type the following command in a terminal of MethFlowVM:

get_test_datasets

The data is then downloaded to the shared folder path set in .methflowrc. The home folder is used if you did not set a shared folder (strongly not recommended!). These test datasets contain data from nine samples (three individuals, three tissues), a small assembly (chromosomes 12 and 19 of hg38) and the collection of adapters of Trimmomatic.

↑

Set your adapter database

It is advisable to have a single multi-FASTA file with the sequences of adapters and possible contaminants that may contain the sequencing reads. We refer to this file as adapter database. If specified, MethFlow never ask you where is the adapter databases. Otherwise, you will be asked for it every time you use MethFlow. To specify, type:

set_adapter_database vm_path

where vm_path is the path in the virtual machine for the multi-FASTA file you want to use as adapter database.

Do not worry if you want to use a different adapter database for a particular input. You can always change it using the same command or specify another file for this particular input (see Change parameters and launch options).

Finally, you can always overwrite the adapter file path by means of the command line when launching.

↑

Set results base folder

You can also specify a base output folder and MethFlow will not ask for it. To specify, type:

set_output_folder vm_path

where vm_path is the path in the virtual machine for the base output folder you want to use. For each input directory an output folder will be created within the base output folder, with the original name of the input directory and the timestamp (name_timestamp) as name. In this folder all partial and final results will be written.

↑

Launch MethFlow (default arguments)

If you followed the previous sections, you have almost everything you need to launch MethFlow. Otherwise, the program will ask for the missing mandatory parameters. To use MethFlow, type:

MethFlow run

This command will first ask the location of the input data directory. Note that you are asked for the path of a directory, not of a file. This directory should contain all input files of a sample. Next, the assembly directory must be given. Again, note that this is a directory, not a file. It contains the chromosome sequences in fasta format and the corresponding Bowtie indexes as needed by Bismark. Note that the indexes will be generated only once, it the first run.

If you want to use a second assembly, type:

MethFlow run –use_assembly2

This command allows to use the unassigned or ambiguous (or both) reads from the first round to map them against a second assembly. All uniquely mapped reads will be finally used for methylation profiling by MethylExtract.

↑

Analyze the result files

The output folder of every analyzed input sample directory contains a number of folders:

Methylation_Maps folder. With three folders inside:
- MethylExtract folder. It contains between one to three methylation map files, one for each analyzed methylation context (see Change parameters and launch options): CG.output, CHG.output and CHH.output. These files contain the methylation profiling results at a single cytosine resolution: the methylation context, the position on the genome, the number of reads where this cytosine is methylated, the coverage and the sequencing quality. For a full description of this format visit the manual of MethylExtract.
- methylKit folder. The methylation profiling results in methylKit input format. This format can also be used by MethylSig. For a full description of this format visit the manual of methylKit.
- methylKit_plots folder. The methylation ratio and coverage distributions of the input sample, plotted by methylKit. Files are in PDF format. To get these plots in other formats, see Downstream analysis.
SNVs folder. There is only one file here: SNVs.vcf. This file contains the sequence variants detected in the input sample against the reference genome assembly. The VCF format specifications can be seen here.
Logs folder. Contains a folder for each program used during the pipeline. Each of these folders have two logs for every processed file: one log recording the standard output and the other recording the standard error.
CITE.txt file. A text file within the references that you should cite if you use MethFlow, including all references to third-party software used in a particular process.
FASTQ folder (for SRA input files). It stores the FASTQ files converted from the original SRA files. There will be either one or two FASTQ files for each SRA file, depending whether the sequencing reads are single-end or paired-end. Only if input sample is SRA or FASTQ and –adapters_trimmed is not specified.
trimmed_FASTQ folder. It stores the trimmed FASTQ files, i.e. the Trimmomatic output files. Only if input sample format is SRA or FASTQ.
FastQC folder. With two folders inside:
- FASTQ_FASTQC folder. Contains the quality report of the FASTQ files before trimming. There is one report for each FASTQ file. Only if input sample format is SRA or FASTQ.
- trimmed_FASTQ_FastQC folder. Contains the quality report of the FASTQ files after trimming. There is one report for each FASTQ file. Only if input sample is SRA or FASTQ and –adapters_trimmed is not specified.
ambiguous_FASTQ folder. It stores the FASTQ files ambiguously mapped against the first assembly. There is one file for each FASTQ used during alignment against the first assembly. This folder appears if the input format is SRA or FASTQ and a second assembly is used.
unmapped_FASTQ folder. It stores the FASTQ files with unmapped reads against the first assembly. There is one file for each FASTQ used during alignment against the first assembly. This folder appears if the input format is SRA or FASTQ and a second assembly is used.
BAM folder. Contains BAM files coming from alignment against the first assembly and, if applicate, the second assembly. There is only one file for each dataset (paired-end data no longer have two files). BAM files from second assembly alignment are merged, if applicable. Only if input sample format is SRA or FASTQ.
fixed_SAM folder. Contains SAM files after bisulfite bias fixing. There is only one file for each dataset. Only if –bisulfite_bias_fixed is not set on the command line.
BSeQC_plots folder. Contains plots about the bisulfite bias. There is only one folder for each dataset. Only if –bisulfite_bias_fixed is not set on the command line.

↑

The local database

Local settings

The easiest way to use MethFlow is to set the value of certain parameters by means of a setting file. This file can be found within your home: $HOME/.methflowrc (note the dot at the beginning), where $HOME is your home directory (i.e. /home/methflow). It is not listed with ls, except you add the option -a.

This file is a text file that can be edited with any plain text editor such as vim or nano:

nano $HOME/.methflowrc

It should contain four variables:

shared: the path of the shared folder to be used.
assemblies: the path of the assemblies folder (see Set assemblies folder).
adapters: the path of the adapter database.
output: the path of the base output folder.

You can modify the variables. If the specified path does not exist or if the parameter is missing at all, MethFlow will ask again on the command line.

↑

Input data

It is highly recommended to provide the input data from the shared folder. The data from different samples must go into separate folders. The input files located within the same folder are interpreted as different runs from the same sample. Accepted formats are SRA, FASTQ, SAM and BAM.

The directory with the sample dataset to be used can be specified in a configuration file (not to be confused with the settings file, see Local settings) or on the command line when you launch MethFlow (see Change parameters and launch options). Otherwise, you will be asked.

↑

Set assemblies folder

To set the assemblies folder you can edit the .methflowrc file or type:

set_assemblies_folder vm_path

where vm_path is the path in the virtual machine for the assemblies you want to use. It is highly recommended that the assemblies folder is within the shared folder.

If you do not set the assemblies folder a folder named assemblies within the path established for the shared folder will be used.

↑

Prepare the assemblies

Each assembly must go into a separate folder into the assemblies folder. The assembly may consist of a multi-FASTA file or several FASTA files, all contained in the same directory. It may contain or not Bismark Bowtie2 indexes. If not, Bismark Bowtie2 indexes will be calculated during the first usage of the assembly by MethFlow.

The directory with the assembly can be specified in a configuration file (not to be confused with the settings file) or on the command line when you launch MethFlow (see Change parameters and launch options). Otherwise, you will be asked.

You can download some assemblies (including Bismark Bowtie2 indexes) with this command:

get_assemblies

The data is then downloaded to the assemblies folder set in .methflowrc.

↑

Change parameters and launch options

To run the MethFlow pipeline we execute the command MethFlow run [arguments] together with the relevant arguments. If you do not specify any arguments the program will enter in the quick mode, where you will be asked interactively.
There is a second command, MethFlow configure [arguments], which can be used to create a configuration file (not to be confused with the settings file). This command does not launch MethFlow but it generates a configuration file with the parameters specified on the command.
MethFlow can be used in three ways:

Interactive:MethFlow run. The program asks you the mandatory arguments through dialogs.
Configuration file: you indicate arguments to be used in a configuration file created by MethFlow and edited by you. Type MethFlow run –config configuration_file to use this mode, where configuration_file is the configuration file previously generated by MethFlow or edited by you.
Command line:MethFlow run [arguments]. The arguments are given when launching the program. If any mandatory arguments are missing you will be asked interactively. It can be combined with the configuration file mode. In case of conflict, the command line value of the conflictive argument will be used.

↑

Mandatory Arguments

Some parameters must be indicated by the user:

input: the path of the input data folder. It must be indicated in a configuration file, on the command line or through a dialog. During the pipeline various arguments are detected: format of the input files, if they have single-end or paired-end reads, if they use phred33 or phred64 and the maximum and the minimum read length.
adapters: the path of the adapter database. It must be indicated in the settings file, in a configuration file, on the command line or through a dialog.
assembly: the path of the first assembly folder. It must be indicated in a configuration file, on the command line or through a dialog. During the pipeline it is checked for Bismark Bowtie2 indexes. If there are not indexes within the folder, they will be calculated.
output: the path of the base output folder. It must be indicated in the settings file, in a configuration file, on the command line or through a dialog.

If –use_assembly2 is used there will be two extra mandatory arguments:

assembly2: the path of the second assembly folder. It must be indicated in a configuration file, on the command line or through a dialog. During the pipeline it is checked for Bismark Bowtie2 indexes. If there are not indexes within the folder, they will be calculated.
use_assembly2_for: indicates which kind of reads will be used for the mapping against the second assembly (ambiguously mapped reads against first assembly, unmapped reads or both). It must be indicated in a configuration file, on the command line or through a dialog. For example, to use it at the command line:

MethFlow run –assembly2 vm_path –use_assembly2_for [ambiguous, unmapped or both]

where vm_path is the path in the virtual machine for the assemblies you want to use. It is highly recommended that the assemblies folder is within the shared folder.

In addition, in this command, you have to chooses using ambiguous, unmapped or both kinds of reads.

↑

Optional Arguments

Most arguments are optional. When they are not given, MethFlow either calculates or tries to estimate them (like minimum_read_length and threads) or it uses the default values. Note that when the parameters are use on the command line, ‘–‘ must precede the parameter name. For example, parameter name: ‘adapter_trimmed’ ➜ on command line: –adapter_trimmed.

adapter_trimmed: it is a bool argument (default: off). When on, the adapters trimming is skipped.
bisulfite_bias_fixed: it is a bool argument (default: off). When on, the bisulfite bias fixing is skipped.
library: indicates whether the type of sequencing library is directional, non-directional or PBAT (options: directional, non_directional or pbat;default: directional). Unfortunately, this argument cannot be estimated before aligning. If you observe a high number of unmapped reads, try changing this argument.
rrbs: it is a bool argument (default: off). Indicate that the sequencing technique used is Reduced Representation Bisulfite Sequencing (RRBS). It takes into account when bisulfite bias fixing. It is recommended to use combined with the argument not_remove_duplicate.
not_seed_mismatch: it is a bool argument (default: off). When on, you do not use mismatches in seed during aligning. When off, you use one mismatch.
seed_length: indicate the length of the seed used during aligning (minimum: 8; maximum: 32; default: 32).
not_remove_duplicate: it is a bool argument (default: off). When on, duplicate reads are not remove during profiling. Recommended for RRBS data.
minimum_phred_score: indicate the minimum accepted phred score during trimming and profiling (default: 20). To set separately for both steps, use advanced arguments (see Manipulate the configuration file).
minimum_read_length: indicate the minimum accepted read length during trimming (default: calculated as half of the original length of the reads).
minimum_coverage: indicate the minimum accepted coverage during profiling (default: 1).
methylation_context: indicate the methylation context to analysis during profiling (options: CG, CHG, CHH or ALL; default: CG)
threads: indicate the maximum number of threads to be used (default: calculated as the number of CPUs of the virtual machine; minimum: 2).

↑

Use command line arguments

All arguments described above can be used in the command to run or configure MethFlow, adding a double hyphen before the name of the argument. For example:

–input vm_path, –adapter_trimmed, –library non_directional or –threads 16.

If you run MethFlow without specify a configuration file, you will be asked for all mandatory arguments not specified on the command line or on the settings file. Optional arguments not indicated will take their default value.

↑

Manipulate the configuration file

The configuration file is the most convenient way to specify the arguments of MethFlow. To create a new configuration file, type:

MethFlow configure [optional: command line arguments]

You can add –config configuration_file to indicate where to save the configuration file. Otherwise, you will be asked. All the specified command line arguments will be directly written in the configure file. You will be asked for mandatory arguments not specified.

If you indicate a previous configuration file with option –config configuration_file and add other command line arguments when using MethFlow configure, the indicated configuration file will be updated with these command line arguments.

If you do the same with MethFlow run the configuration file will not be updated, but for arguments that are in both the configuration file and the command line the value of the command line is taken. This way you can easily launch the same configuration file for other input, for example:

MethFlow run –config configuration_file –input input_folder

The configuration file is in YAML format. Each argument is specified as follows:

argument_name: argument_value

(True or False for bool arguments; None if not used)

This configuration file contains a special argument: steps. This argument is a list of pseudo-commands. Each pseudo-command represents a command to be executed during the pipeline, listed in the order they will be used.

Use pseudo-commands to specified advanced arguments. Within each pseudo-command there are arguments that begin with @, # or $. Arguments that begin with @ or # cannot be changed. If modified, the pipeline will not work. The arguments starting with $ may be changed, but the $ character should be removed.

You can also add new arguments as would be used in the external program indicated in the pseudo-commando (e.g. Trimmomatic). If you add new arguments, please read the manuals of the external programs and make sure that they are supported or recommended for the analysis carried out by the pipeline. For example, never use the LEADING argument of Trimmomatic if you want to fix the bisulfite bias later because this option removes bases from 5’ end of the reads and BSeQC will not be able to check properly the methylation per read position distribution.

You can download a configuration file commented in detail here.

↑

External Program Manuals and Documentation

fastq-dump (SRA Toolkit): we use this program to convert files in SRA format to FASTQ format. Documentation.
Trimmomatic: we use this program to remove the adapter and low quality bases at the 3’ end. Manual.
Bismark: we use this program to align reads against three-letter reference assemblies. Manual.
Bowtie2: it is the aligner that we use in Bismark. Manual.
BSeQC: we use this program to fix the bisulfite bias due to technical factors. Documentation.
MethylExtract: the core of MethFlow. We use this program to profile methylations levels and single nucleotide variants. Manual.
FastQC: the program that we use to check the quality of FASTQ and trimmed FASTQ files. Documentation.
methylKit: the main program that we use in downstream analysis. Manual.

↑

Downstream analysis

methylKit downstream analysis

Using methylKit you can do a lot of downstream analysis, as compare methylation maps of different samples by means of a Pearson correlation matrix, get differential methylation data at level of single cytosine (DMCs) and analysis the composition in DMCs of a wide variety of annotations.

MethFlow converts automatically the methylation maps from MethylExtract output format to methylKit input format during the pipeline (see Analyze the result files). Anyway, you can convert MethylExtract output files into methylKit input files anytime by typing this command:

me2mk -i MethylExtract_Output_File -o methylKit_Input_File -c Methylation_Context (CG, CHG or CHH) [–destrand]

-i, -o and -c are mandatory arguments. Optionally, you can use –destrand to merge the data from both Watson and Crick strands (default: off).

To do a quick descriptive analysis using methylKit you can use the following commands:

Methylation ratio distribution:

mk_methRatio_distribution -i methylKit_Input_File -o Image_Ouput_File [-f format (pdf, ps, svg, png, jpeg, bmp or tiff)] [-a assembly] [-m methylation_context]

-i and -o are mandatory arguments. -f takes pdf as default value.

Coverage distribution:

mk_coverage_distribution -i methylKit_Input_File -o Image_Ouput_File [-f format (pdf, ps, svg, png, jpeg, bmp or tiff)] [-a assembly] [-m methylation_context]

-i and -o are mandatory arguments. -f takes pdf as default value.

Pearson Correlation Matrix:

mk_pearson_correlation -i methylKit_Input_Files -o Image_Ouput_File [-f format (pdf, ps, svg, png, jpeg, bmp or tiff)] [-a assembly] [-m methylation_context]

-i and -o are mandatory arguments. -f takes pdf as default value. You should indicate more than one input file, separated by spaces.

Clustering Tree:

mk_clustering -i methylKit_Input_Files -o Image_Ouput_File [-f format (pdf, ps, svg, png, jpeg, bmp or tiff)] [-a assembly] [-m methylation_context]

-i and -o are mandatory arguments. -f takes pdf as default value. You should indicate more than one input file, separated by spaces.

Principal Component Analysis:

mk_pca -i methylKit_Input_Files -o Image_Ouput_File [-x a_PC_for_x-axis] [-y another_PC_for_y-axis] [–screenplot] [-f format (pdf, ps, svg, png, jpeg, bmp or tiff)] [-a assembly] [-m methylation_context]

-i and -o are mandatory arguments. -f takes pdf as default value. -x and -y takes 1 and 2 as default values, respectively. You should indicate more than one input file, separated by spaces. You can add the optional argument –screenplot to get the screenplot. Otherwise you get the PC indicated in x versus PC and indicated.

Finally, you can also calculate DMCs using this command:

mk_diffmeth -i methylKit_Input_Files -o Image_Ouput_File [-p number_of_threads] [-a assembly] [-m methylation_context]

-i and -o are mandatory arguments. -p takes 4 as default value. You should indicate more than one input file, separated by spaces.

For further details on the output, see the manual of methylKit.

↑

Convert to BED and other formats

In addition to methylKit intput format, you can convert MethylExtract output files to other formats such as BedGraph, BED6 or bigWig:

BedGraph, BED6 and BED6+6:

me2bed -i MethylExtract_Output_File -o BED_Output_File -f Output_Format (bedgraph, bed6, bed6+6 or ucsc) -c Methylation_Context (CG, CHG or CHH)

-i, -o, -f and -c are mandatory arguments.

You can find the specifications of BED and BedGraph formats here.

The score column of the BED file and the dataValue column of the BedGraph file contains a numerical value from 0 to 1000. This value is the methylation level, being 0 completely unmethylated and 1000 completely methylated. The six additional columns of the BED6+6 are:

Watson METH: number of reads methylated for this cytosine (referred to the Watson strand).
Watson COVERAGE: reads covering the cytosine in this sequence context (referred to the Watson strand).
Watson QUAL: PHRED score average for the reads covering the cytosine (referred to the Watson strand).
Crick METH: number of reads methylated for this cytosine (referred to the Crick strand).
Crick COVERAGE: reads covering the cytosine in this context (referred to the Crick strand).
Crick QUAL: PHRED score average for the reads covering the cytosine (referred to the Crick strand).

For more details of these values visit the manual of MethylExtract.

bigBed:

First of all, convert your file to BED6 format. Then get the chromosome sizes file from the assembly multi-FASTA file:

faidx multi-FASTA_input_file -i chromsizes > chrom.sizes

Finally, convert the BED6 input format to bigBed:

bedToBigBed -type=bed6 bed6_input_file chrom.sizes bigBed_output_file

You can find the specifications of bigBed format here.

bigWig:

First of all, convert your file to BedGraph format. Then get the chromosome sizes file from the assembly multi-FASTA file:

faidx multi-FASTA_input_file -i chromsizes > chrom.sizes

Finally, convert the BED6 input format to bigWig:

bedGraphToBigWig bedGraph_input_file chrom.sizes bigWig_output_file

You can find the specifications of bigWig format here.

↑

Send to a Galaxy Server

Galaxy gives us the opportunity to do a myriad of analysis. An easy way to send your data to Galaxy from inside virtual machine is using the helper tool upload2galaxy.

To use this tool, first of all you need to get your API key from the Galaxy Server you want to use. To get your API key, open a web browser, go to the URL of the Galaxy Server you want to use and login. In the top menu, go to

User ➜ Preferences ➜ Manage your API keys

Now click on the button Generate a new key now and copy the Current API key. You will use this key to send files to your Galaxy account.

To send a file to your Galaxy account, type:

upload2galaxy [-u URL_of_a_Galaxy_Server] -k API_Key_of_your_Galaxy_Account -i Path_of_the_File_to_Upload [-n Name_of_the_Galaxy_History_to_be_created]

-i and -k are mandatory arguments. By default -u is the URL of the Galaxy Main Server and -n is MethFlow.

To check your uploaded file, in the Galaxy website go to

User ➜ Saved Histories

And click on MethFlow or the name indicated with -n.

↑

Upload to UCSC Genome Browser

One of the best ways to visualize your data is by UCSC Genome Browser. There you can view your data along chromosomes and compare with a myriad of other genomic annotations.

By following these instructions, you can upload your files to UCSC Genome Browser:

Convert your data to UCSC BED6 format by typing:

me2bed -i MethylExtract_Output_File -o BED_Output_File -f ucsc -c Methylation_Context (CG, CHG or CHH)

Note that -f ucsc is required for you to visualize your data correctly.

Open a browser and go to the UCSC Genome Browser website
Go to My Data ➜ My Sessions in the top menu and login or create an account (your data will continue online after logout).
Once logged, go to My Data ➜ Custom Tracks in the top menu.
Select a file from your local disk and submit it.
Once uploaded, select view in Genome Browser and click on go.

Now you can browse your data. If you want to upload more files:

Go to My Data ➜ Custom Tracks in the top menu.
Click on add custom tracks.
Select a file from your local disk and submit it.
Once uploaded, select view in Genome Browser and click on go.

↑