-
Notifications
You must be signed in to change notification settings - Fork 10
Bulk RNA Seq Nanopore tutorial
This tutorial will go through an end to end analysis for RNA-Seq analysis with Nanopore data using Eoulsan. Before you start, you will need an installed copy of Eoulsan (see the installation guide).
This tutorial will use an example mouse RNA-Seq Nanopore dataset generated for the "Nanopore long-read RNAseq reveals widespread transcriptional variation among the surface receptors of individual B cells" article by Byrne et al. This analysis do not produce any differentialy expressed genes as the sample are biological replicates. The purpose of this tutorial is only to show how to handle RNA-Seq Nanopore data with Eoulsan.
The analysis of Nanopore RNA-Seq data is very similar to Illumina RNA-Seq data. There is very few changes between the two workflows: we replace the STAR mapper by Minimap2 and the parameter of the filterread, filtersam and expression steps are just slightly different. Most of the step
Warning: You need about 40 GB of RAM to perform this tutorial.
Outline of the process:
| Step | Module used | Inputs | Output |
|---|---|---|---|
| Create a genome index for Minimap2 | minimap2indexgenerator | FASTA | ZIP |
| Quality control of the raw FASTQ files with FastQC | fastqc | FASTQ | HTML |
| Filter raw FASTQ files | filterreads | FASTQ | FASTQ |
| Map the reads using Minimap2 | mapreads | FASTQ | SAM |
| Remove unmap and multimatches SAM entries | filtersam | SAM | SAM |
| Convert SAM files into sorted and indexed BAM files | sam2bam | SAM | BAM |
| Convert BAM files int BedGraph files | bam2bedgraph | BAM | BedGraph |
| Convert BAM files int BED files | splice2bed | BAM | BED |
| Count alignments with HTSeq-count | expression | SAM | TSV |
| Convert expression files into XLSX files | expressionresultsannotation | TSV | XLSX |
| Quality control of the analysis with MultiQC | multiqc | FASTQ | HTML |
| Normalisation and differential analysis with DESeq 2 | normalization | TSV | TSV + graphics |
| Convert differential analysis output files into XLSX files | expressionresultsannotation | TSV | XLSX |
# Step 1: get data (Raw
$ wget http://outils.genomique.biologie.ens.fr/leburon/downloads/eoulsan-samples/rnaseq-nanopore/SRR5286960.fastq.bz2 && \
wget http://outils.genomique.biologie.ens.fr/leburon/downloads/eoulsan-samples/rnaseq-nanopore/SRR5286961.fastq.bz2 && \
wget http://outils.genomique.biologie.ens.fr/leburon/downloads/eoulsan-samples/rnaseq-nanopore/SRR5286962.fastq.bz2 && \
wget http://outils.genomique.biologie.ens.fr/leburon/downloads/eoulsan-samples/rnaseq-nanopore/SRR5286963.fastq.bz2 && \
wget http://outils.genomique.biologie.ens.fr/leburon/downloads/eoulsan-samples/genomes/mm10ens91.fasta.bz2 && \
wget http://outils.genomique.biologie.ens.fr/leburon/downloads/eoulsan-samples/genomes/mm10ens91.gff.bz2 && \
wget http://outils.genomique.biologie.ens.fr/leburon/downloads/eoulsan-samples/genomes/mm10ens91.tsv.bz2
# Step 2: create the design file and set the biological replicates
$ eoulsan.sh createdesign *.fastq.bz2 mm10ens91.fasta.bz2 mm10ens91.gff.bz2 mm10ens91.tsv.bz2 && vim design.txt
# Step 2 bis: download a ready to use design file
$ wget -O design.txt https://github.com/GenomicParisCentre/eoulsan/wiki/files/design-rnaseq-nanopore.txt
# Step 3: create the workflow file or use an existing one
$ wget https://github.com/GenomicParisCentre/eoulsan/wiki/files/workflow-rnaseq-nanopore.xml
# Step 4: launch the analysis
$ eoulsan.sh exec workflow-rnaseq-nanopore.xml design.txtWe provide in this section some samples files to test Eoulsan. These files have been produced during a mouse RNA-Seq experiment sequenced on an Illumina GAIIx.
-
Read files from two different conditions
-
Mouse genome file
-
Feature annotation file
-
Additional annotation file
-
Worflow file
- workflow-rnaseq-nanopore.xml. This file is an already created workflow file to analyze previous reads files.
In an empty directory, copy the reads, genome and annotation files, then you can create a design file with the next command:
$ eoulsan.sh createdesign *.fastq.bz2 mm10ens91.fasta.bz2 mm10ens91.gff.bz2 mm10ens91.tsv.bz2You can now modify the design file to add additional information. To avoid failure of the deseq2 step, biological replicate must be defined. To do that in the column named "Exp.exp1.Condition" set "condition1" for samples SRR5286960/SRR5286961 and "condition2" for samples SRR5286962/SRR5286963. Note that Eoulsan handle automatically compressed files.
You can also download a ready to use design file with the following command:
$ wget -O design.txt https://github.com/GenomicParisCentre/eoulsan/wiki/files/design-rnaseq-nanopore.txtTo create a workflow file, the best solution is to reuse an existing workflow file and adapt it to your needs. You can download an existing workflow to analyse this data with the following command:
$ wget https://github.com/GenomicParisCentre/eoulsan/wiki/files/workflow-rnaseq-nanopore.xmlThe workflow file contains the list of the steps that will be executed by Eoulsan. Each step have parameters and it is related to a module to execute. Some step modules are only available in local mode (like differential analysis) or in distributed mode, see the modules page for more details. For each step you can change, add or remove parameters. Parameters are specific to a module, consult the documentation of the built-in steps for the list of available parameters of each step.
At least, there is a global section in the workflow file that override the values of Eoulsan configuration file. This section is useful to set for example the path of the temporary directory to use.
Unlike other mapper (e.g. Bowtie/BWA), Minimap2 creation index can be customized using parameters. This step allow to create a custom index for Minimap2. If you don't need a custom Minimap2 index, you can remove this step from your workflow. The following configuration example will enable splicing index mode.
<!-- Create a custom Minimap2 index feature annotation -->
<step skip="false">
<module>minimap2indexgenerator</module>
<parameters>
<!-- Create an index for splicing mode of Minimap2 -->
<parameter>
<name>indexer.arguments</name>
<value>-x splice</value> <!-- cDNA -->
<!--value>-x splice -uf -k14</value--> <!-- Direct RNA -->
</parameter>
</parameters>
</step>This step allow to launch FastQC on the raw data. To enable this step, just add the following lines in your workflow file. More parameters are available in the FastQC module, please see the FastQC module page in the reference manual for more informations.
<!-- FastQC of raw reads -->
<step id="rawfastqc" skip="false">
<module>fastqc</module>
</step>This step allow to remove bad quality reads.
The following configuration example will remove reads that do not pass Illumina filter, trim polyN reads and remove all the reads with a length lower than 40 and with a mean quality lower to 30.
For more information about the filterreads module see the reference documentation.
<!-- Filter reads -->
<step skip="false" discardoutput="asap">
<module>filterreads</module>
<parameters>
<!-- Remove polyN tails of the reads -->
<parameter>
<name>trimpolynend</name>
<value></value>
</parameter>
</parameters>
</step>This step allow to map FASTQ read. The following configuration example will use STAR with default parameters.
STAR is only parameted to keep unmapped read in output SAM files. Eoulsan is bundled with many mappers (Bowtie, Bowtie2, BWA, Minimap, Gmap/Gsnap, STAR...). For more information about the mapreads module see the reference documentation.
Note that with STAR or Minimap2, custom indexes can used.
To create this custom indexes, you need to add a dedicated step (the starindexgenerator module for STAR and the minimap2indexgenerator module for Minimap2.
<!-- Mapping of the reads -->
<step skip="false" discardoutput="asap">
<module>mapreads</module>
<parameters>
<!-- Use use Minimap2 as mapper -->
<parameter>
<name>mapper</name>
<value>minimap2</value>
</parameter>
<!-- The version of Minimap2 to use is 2.12 -->
<parameter>
<name>mapper.version</name>
<value>2.12</value>
</parameter>
<!-- We use Minimap2 in splice mode -->
<parameter>
<name>mapper.arguments</name>
<value>-x splice</value> <!-- cDNA -->
<!--value>-x splice -uf -k14</value--> <!-- Direct RNA -->
</parameter>
</parameters>
</step>In this step, we will remove all SAM entries that are unmapped or multimatches.
Like the filterreads module, this step is highly configurable with additional plugins.
For more information about the filtersam module see the reference documentation.
<!-- Filtering of the SAM alignments -->
<step skip="false" discardoutput="asap">
<module>filtersam</module>
<parameters>
<!-- Remove the unmap entries in the files -->
<parameter>
<name>removeunmapped</name>
<value>true</value>
</parameter>
<!-- Remove alignmments with poor quality -->
<parameter>
<name>quality.threshold</name>
<value>1</value>
</parameter>
<!-- Remove supplementary alignments -->
<parameter>
<name>removesupplementary</name>
<value>true</value>
</parameter>
<!-- Remove the multimaches alignments -->
<parameter>
<name>removemultimatches</name>
<value>true</value>
</parameter>
</parameters>
</step>With this optional step, we will convert SAM files converted in the previous step in BAM files.
The generated BAM files are sorted by coordinates. For each created BAM, an index file (BAI file) is also created.
For more information about the sam2bam module see the reference documentation.
<!-- Convert SAM to BAM file and sort the BAM file by coordinate -->
<step skip="false">
<module>sam2bam</module>
<parameters>
<!-- Compression level of the BAM file -->
<parameter>
<name>compression.level</name>
<value>5</value>
</parameter>
</parameters>
</step>In this optional step, BAM files will be converted in BedGraph files.
Note that this step require the bam-to-bedgraph.xml Galaxy tool wrapper that is available on the eoulsan-tool GitHub repository.
<!-- Create Bedgraph files from the BAM files -->
<step dataproduct="match" skip="false">
<module>bam2bedgraph</module>
<parameters/>
</step>With this optional step, positions of the alignments will be converted in BED 12 files.
For more information about the splice2bed module see the reference documentation.
<!-- Create Bed files from the SAM files -->
<step dataproduct="match" skip="false">
<module>splice2bed</module>
<parameters/>
</step>This step allow to compute the expression of the gene/transcripts.
Currently only the htseq-count counter is implemented in Eoulsan.
Most of the option of this module are the standard options of HTSeq-count.
For more information about this expressionmodule, see the HTSeq-count and module documentation.
<!-- Compute expression -->
<step skip="false">
<module>expression</module>
<parameters>
<!-- We use the fast Java implementation of htseq-count to perform the counting -->
<parameter>
<name>counter</name>
<value>htseq-count</value>
</parameter>
<!-- The output file file will be in TSV format -->
<parameter>
<name>output.file.format</name>
<value>tsv</value>
</parameter>
<!-- The input feature annotation file is in GFF format -->
<parameter>
<name>features.file.format</name>
<value>gff</value>
</parameter>
<!-- The feature type to use is "exon", we count by exon -->
<parameter>
<name>genomic.type</name>
<value>exon</value>
</parameter>
<!-- We aggregate the counts by transcript -->
<parameter>
<name>attribute.id</name>
<value>Parent</value>
</parameter>
<!-- The library was not oriented -->
<parameter>
<name>stranded</name>
<value>no</value>
</parameter>
<!-- The use the union mode as overlap mode -->
<parameter>
<name>overlap.mode</name>
<value>union</value>
</parameter>
<!-- We do not remove ambiguous cases -->
<parameter>
<name>remove.ambiguous.cases</name>
<value>false</value>
</parameter>
</parameters>
</step>The expression files generated in the previous step are simple tabulated text file with only two columns: the identifiers of the feature and the counts.
This files are often opened with spreadsheet software like Microsoft Excel.
Unfortunately this tools can misinterpret the content of the file.
That's why with this step the TSV files are converted in XLSX files (or ODS for LibreOffice).
Moreover, this step add several column with useful information (e.g. gene name, gene description...) that is coming from the additional annotation file.
For more information about this expressionresultsannotation module, see the module documentation.
<!-- Convert expression files to XLSX files -->
<step skip="false">
<module>expressionresultsannotation</module>
<parameters>
<!-- Select the output of the step, ODS and TSV format are also supported -->
<parameter>
<name>output.format</name>
<value>xlsx</value>
</parameter>
</parameters>
</step>MultiQC is tool that has been developed to aggregate all the FastQC reports of an analysis in one report.
Then, it has been enhanced to support many other tools.
With this step, we will generate a MultiQC report with data that come from the FastQC, mapping and expression steps.
For more information about this multiqc module, see the module documentation.
<!-- MultiQC -->
<step skip="false">
<module>multiqc</module>
<parameters>
<!-- MultiQC will contain reports from FastQC, STAR and HTSeq-count -->
<parameter>
<name>reports</name>
<value>fastqc,mapreads,expression</value>
</parameter>
<!-- Docker is required to launch this step -->
<parameter>
<name>use.docker</name>
<value>true</value>
</parameter>
</parameters>
</step>With this step, the expression data will be normalized and then the differential analysis will be performed.
This module allow to use both standard and complex experimental design.
The design file documentation, explain how to define the experimental design of your experiment.
This module can only handle experimental design with biological replicates.
If you don't have any biological replicates, please use instead the DESeq1 module.
The r.execution.mode parameter allow to choose the method to use for launching R (local, RServe or Docker).
We recommend you to use docker.
For more information about this deseq2 module, see the module documentation.
Note: You must define the main.docker.uri parameter in your global section of your workflow file (or in your configuration file) to enable docker.
Usually the value of this parameter id unix:///var/run/docker.sock.
<!-- Normalization and differential analysis with DESeq 2 -->
<step skip="false">
<module>deseq2</module>
<parameters>
<!-- Generate normalization plots -->
<parameter>
<name>norm.fig</name>
<value>true</value>
</parameter>
<!-- Genearte differential analysis plots -->
<parameter>
<name>diffana.fig</name>
<value>true</value>
</parameter>
<!-- Enable Normalization -->
<parameter>
<name>norm.diffana</name>
<value>true</value>
</parameter>
<!-- Enable differential analysis -->
<parameter>
<name>diffana</name>
<value>true</value>
</parameter>
<!-- Set the size factor type to "ratio" -->
<parameter>
<name>size.factors.type</name>
<value>ratio</value>
</parameter>
<!-- Set the fit type to "parametric" -->
<parameter>
<name>fit.type</name>
<value>parametric</value>
</parameter>
<!-- Use Wald as statistical test -->
<parameter>
<name>statistical.test</name>
<value>Wald</value>
</parameter>
<!-- Use Docker to execute this step -->
<parameter>
<name>r.execution.mode</name>
<value>docker</value>
</parameter>
<!-- Define the Docker image to use (optional) -->
<parameter>
<name>docker.image</name>
<value>bioconductor/release_sequencing:3.1</value>
</parameter>
</parameters>
</step>The text files generated by the previous step are simple tabulated text files.
This files are often opened with spreadsheet software like Microsoft Excel.
Unfortunately this tools can misinterpret the content of the text file.
That's why with this step the TSV files are converted in XLSX files (or ODS for LibreOffice).
Moreover, this step add several column with useful information (e.g. gene name, gene description...) that is coming from the additional annotation file.
For more information about this diffanaresultsannotation module, see the module documentation.
<!-- Convert ourput of the differential analysis step to XLSX files -->
<step skip="false">
<module>diffanaresultsannotation</module>
<parameters>
<!-- Select the output of the step, ODS and TSV format are also supported -->
<parameter>
<name>output.format</name>
<value>xlsx</value>
</parameter>
</parameters>
</step>This section of the workflow file define the global configuration of the analysis. For this analysis, path to additional data format, Galaxy tools and Docker connection URI are set.
Note: You can also define this global parameters in a configuration file like the default Eoulsan configuration file ~/.eoulsan, see the configuration file documentation for more information.
<!-- This section contains globals parameters that can also be define in a configuration file. -->
<globals>
<!-- Define the location of the Galaxy tools/modules to use -->
<parameter>
<name>main.galaxy.tool.path</name>
<value>https://raw.githubusercontent.com/GenomicParisCentre/eoulsan-tools/master/galaxytools</value>
</parameter>
<!-- Define the location of additional formats to use -->
<parameter>
<name>main.format.path</name>
<value>https://raw.githubusercontent.com/GenomicParisCentre/eoulsan-tools/master/formats</value>
</parameter>
<!-- Define the location of the Docker connection -->
<parameter>
<name>main.docker.uri</name>
<value>unix:///var/run/docker.sock</value>
</parameter>
</globals>Once your design file and workflow file are ready, you can launch Eoulsan analysis with the following command:
$ eoulsan.sh exec workflow-local.xml design.txtWarning: To perform the normalization and differential analysis steps of this workflow, this demo requires Docker. If you want to run this demo without Docker, you must must install R (or a RServe server) and the related packages (See differential analysis module documentation for more information) and change the execution.mode parameter of the normalization and diffana steps.
Once started and before starting analysis, Eoulsan will check if:
- The design file and workflow file are valid
- All the modules related to the steps exist
- The workflow of all the steps are valid
- The order of the steps is correct (a step can not use data generated after its end)
- The input files are valid (fasta and annotation files)
- A genome mapper index is needed. In this case, a step to generate this index is automatically added
If successful, you will obtain a new directory name like eoulsan-20110310-101235 with log files about the analysis.
Results files are stored in the current directory of the user.