kfukuta/wf-transcriptomes-v202
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Remove fusion detection [CW-4552]

Browse Source
This commit is contained in:
Neil Horner 2024-08-06 20:20:19 +00:00 committed by Sam Nicholls
parent a037538744
commit 0200514536
20 changed files with 40 additions and 450 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
.gitlab-ci.yml
View File

@ -10,7 +10,6 @@ variables:
--transcriptome-source ${CI_PROJECT_NAME}/data/reference-guided \
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/gencode.v22.annotation.chr20.gtf \
--direct_rna --minimap2_index_opts '-k 15' --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \
--jaffal_refBase ${CI_PROJECT_NAME}/data/differential_expression/chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
CI_FLAVOUR: "new"
@ -50,7 +49,7 @@ docker-run:
parallel:
matrix:
- MATRIX_NAME: [
"fusions", "differential_expression", "isoforms", "isoforms_bam",
"differential_expression", "isoforms", "isoforms_bam",
"only_differential_expression", "differential_expression_gff3",
"ncbi_gzip", "ncbi_no_gene_id", "ensembl_with_versions",
"differential_expression_mouse", "no_ref_annotation",
@ -83,14 +82,6 @@ docker-run:
--ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
NF_IGNORE_PROCESSES: run_gffcompare,check_annotation_strand,preprocess_reads,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "fusions"
variables:
NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/wf-isoforms_test_data.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/wf-isoforms_test_data.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config
NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB --fastq ${CI_PROJECT_NAME}/data/ERR6053095_chr20.fastq --transcriptome-source reference-guided \
--ref_genome ${CI_PROJECT_NAME}/data/chr20/hg38_chr20.fa --ref_annotation ${CI_PROJECT_NAME}/data/chr20/gencode.v22.annotation.chr20.gtf \
--jaffal_refBase ${CI_PROJECT_NAME}/data/chr20/ --jaffal_genome hg38_chr20 --jaffal_annotation genCode22 \
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
NF_IGNORE_PROCESSES: preprocess_reads,check_annotation_strand,merge_transcriptomes,decompress_annotation,decompress_ref,decompress_transcriptome,preprocess_ref_transcriptome
- if: $MATRIX_NAME == "differential_expression"
variables:
NF_BEFORE_SCRIPT: mkdir -p ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf ${CI_PROJECT_NAME}/data/differential_expression.tar.gz -C ${CI_PROJECT_NAME}/data/ && wget -O ${CI_PROJECT_NAME}/data/demo.nextflow.config https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/demo.nextflow.config

6
CHANGELOG.md
View File

@ -4,6 +4,12 @@ All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
## [Unreleased]
### Removed
- Fusion detection subworkflow, as the functionality is not robust enough for general use at this time.
### Changed
- Updated pychopper to 2.7.10
## [v1.2.1]
### Changed
- Increase some memory and CPU allocations.

41
README.md
View File

@ -1,6 +1,6 @@
# Workflow Transcriptomes
Transcriptome analysis including assembly and annotation of cDNA and direct RNA sequencing data, gene fusions and differential expression.
Transcriptome analysis including assembly and annotation of cDNA and direct RNA sequencing data and differential expression analysis.
@ -11,7 +11,6 @@ This workflow can be used for the following:
+ Identify RNA transcripts using either cDNA or direct RNA reads.
+ Reference aided transcriptome assembly.
+ Annotation of assembled transcripts.
+ Gene fusions detection.
+ Differential gene expression analysis using a pre-computed or assembled reference transcriptome.
+ Differential transcript usage analysis using a precomputed or assembled reference transcriptome.
@ -89,9 +88,6 @@ nextflow run epi2me-labs/wf-transcriptomes \
--de_analysis \
--direct_rna \
--fastq 'wf-transcriptomes-demo/differential_expression_fastq' \
--jaffal_annotation 'genCode22' \
--jaffal_genome 'hg38_chr20' \
--jaffal_refBase 'wf-transcriptomes-demo/chr20' \
--minimap2_index_opts '-k15' \
--ref_annotation 'wf-transcriptomes-demo/gencode.v22.annotation.chr20.gtf' \
--ref_genome 'wf-transcriptomes-demo/hg38_chr20.fa' \
@ -144,7 +140,7 @@ input_reads.fastq ─── input_directory ─── input_directory
|--------------------------|------|-------------|------|---------|
| fastq | string | FASTQ files to use in the analysis. | This accepts one of three cases: (i) the path to a single FASTQ file; (ii) the path to a top-level directory containing FASTQ files; (iii) the path to a directory containing one level of sub-directories which in turn contain FASTQ files. In the first and second case, a sample name can be supplied with `--sample`. In the last case, the data is assumed to be multiplexed with the names of the sub-directories as barcodes. In this case, a sample sheet can be provided with `--sample_sheet`. | |
| bam | string | BAM or unaligned BAM (uBAM) files to use in the analysis. | This accepts one of three cases: (i) the path to a single BAM file; (ii) the path to a top-level directory containing BAM files; (iii) the path to a directory containing one level of sub-directories which in turn contain BAM files. In the first and second case, a sample name can be supplied with `--sample`. In the last case, the data is assumed to be multiplexed with the names of the sub-directories as barcodes. In this case, a sample sheet can be provided with `--sample_sheet`. | |
| transcriptome_source | string | Select how the transcriptome used for analysis should be prepared. | To analyse only gene fusions and differential expression use of an existing transcriptome may be preferred and so 'precomputed' should be selected. In this case the 'ref_transcriptome' parameter should be specified. To create a reference transcriptome using an existing reference genome, select 'reference guided' and specify the 'ref_genome' parameter. | reference-guided |
| transcriptome_source | string | Select how the transcriptome used for analysis should be prepared. | For differential expression analysis, use of an existing transcriptome may be preferred and so 'precomputed' should be selected. In this case the 'ref_transcriptome' parameter should be specified. To create a reference transcriptome using an existing reference genome, select 'reference guided' and specify the 'ref_genome' parameter. | reference-guided |
| ref_genome | string | Path to reference genome sequence [.fa/.fq/.fa.gz/fq.gz]. Required for reference-based workflow. | A reference genome is required for reference-based assembly of a transcriptome. | |
| ref_transcriptome | string | Transcriptome reference file. Required for precomputed transcriptome calculation and for differential expression analysis. | A reference transcriptome related to the sample under study. Must be supplied when the 'Transcriptome source' parameter has been set to 'precomputed' or to perform differential expression. | |
| ref_annotation | string | A reference annotation in GFF2 or GFF3 format (extensions .gtf(.gz), .gff(.gz), .gff3(.gz)). Only annotation files from [Encode](https://www.encodeproject.org), [Ensembl](https://www.ensembl.org/index.html) and [NCBI](https://www.ncbi.nlm.nih.gov/) are supported. | This will be used for guiding the transcriptome assembly and to label transcripts with their corresponding gene identifiers. Note: If in de_analysis mode transcript strands must be only + or -. | |
@ -179,15 +175,6 @@ input_reads.fastq ─── input_directory ─── input_directory
| stringtie_opts | string | Extra command-line options for stringtie transcript assembly. | For additional String tie options see [here](https://github.com/gpertea/stringtie#stringtie-options). | --conservative |
### Gene Fusion Detection Options
| Nextflow parameter name | Type | Description | Help | Default |
|--------------------------|------|-------------|------|---------|
| jaffal_refBase | string | JAFFAl reference genome directory. | JAFFAL human hg38 reference data directory can be downloaded from here: https://figshare.com/ndownloader/files/25410494 or see the README for alternative instructions. If custom gemome files are required, see the instructions here: https://github.com/Oshlack/JAFFA/wiki/FAQandTroubleshooting#how-can-i-generate-the-reference-files-for-a-non-supported-genome. | |
| jaffal_genome | string | Genome reference prefix. e.g. hg38. | JAFFAL reference files are prefixed with the genome reference file name and need to be supplied . If using the human reference data provided by JAFFAL, this can be left at `hg38`. | hg38 |
| jaffal_annotation | string | Annotation suffix. | JAFFAL reference files are suffixed with the annotation filename and this needs to be supplied. For the human hg38 reference data supplied by JAFFAL, this is `genCode22`. | genCode22 |
### Differential Expression Options
| Nextflow parameter name | Type | Description | Help | Default |
@ -245,8 +232,6 @@ Output files may be aggregated including information for all samples or provided
| Transcript counts filtered | de_analysis/filtered_transcript_counts_with_genes.tsv | Filtered transcript counts, used for differential transcript usage analysis. Includes a reference to the associated gene ID. | aggregated |
| Transcript info table | {{ alias }}_transcripts_table.tsv | This file details each isoform that was reconstructed from the input reads. It contains a subset of columns from the .tmap output from [gffcompare](https://ccb.jhu.edu/software/stringtie/gffcompare.shtml) | per-sample |
| Final non redundant transcriptome | de_analysis/final_non_redundant_transcriptome.fasta | Transcripts that were used for differential expression analysis including novel transcripts with the identifiers used for DE analysis. | aggregated |
| Fusion transcript sequences | jaffal_output_{{ alias }}/jaffa_results.fasta | Fusion transcript sequences output by Jaffa. | per-sample |
| Fusion transcript sequence summary file | jaffal_output_{{ alias }}/jaffa_results.csv | Fusion transcript sequences summary file output by Jaffa. | per-sample |
@ -280,10 +265,7 @@ Transcript GFF files from the chunks with the same sample aliases will then be m
#### 3.6 Create transcriptomes
[Gffread](https://github.com/gpertea/gffread) is used to create a transcriptome FASTA file from the final GFF as well as a merged transcriptome that includes annotations in the FASTA headers where available.
### 4. Find gene fusions
If gene fusion options are provided, fusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFFA), with the JAFFAL extension. To enable this provide the gene fusion detection options: `jaffal_refBase`, `jaffal_genome` and `jaffal_annotation`.
### 5. Differential expression analysis
### 4. Differential expression analysis
Differential gene expression (DGE) and differential transcript usage (DTU) analyses aim to identify genes and transcripts that show statistically altered expression patterns.
@ -306,33 +288,30 @@ barcode05,sample05,treated
barcode06,sample06,treated
```
#### 5.1 Merge cross sample transcriptomes
#### 4.1 Merge cross sample transcriptomes
If a `ref_transcriptome` is not provided, the transcriptomes created by the workflow will be used for DE analysis. To do this, the GFF outputs of GffCompare are merged using StringTie. A final non redundant FASTA file of the transcripts is created using the merged GFF file and the reference genome using seqkit.
#### 5.2 Create a final non redundant transcriptome
#### 4.2 Create a final non redundant transcriptome
The reads from all the samples will be aligned with the final non redundant transcriptome using Minimap2 in a splice aware manner.
#### 5.3 Count genes and transcripts
#### 4.3 Count genes and transcripts
[Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving gene and transcript counts.
#### 5.4 edgeR based differential expression analysis
#### 4.4 edgeR based differential expression analysis
A statistical analysis is first performed using [edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html) to identify the subset of differentially expressed genes using the gene counts as input. A normalisation factor is calculated for each sequence library using the default TMM method (see [McCarthy et al. (2012)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3378882/) for further details). The defined experimental design is used to calculate estimates of dispersion for each of the gene features. Statistical tests are calculated using the contrasts defined in the experimental design. The differentially expressed genes are corrected for false discovery (FDR) using the method of Benjamini & Hochberg ([Benjamini and Hochberg (1995)](https://www.jstor.org/stable/2346101))
#### 5.5 Pre-filtering of quantitative data using DRIMSeq
#### 4.5 Pre-filtering of quantitative data using DRIMSeq
[DRIMSeq](https://bioconductor.org/packages/release/bioc/html/DRIMSeq.html) is used to filter the transcript count data from the Salmon analysis for differential transcript usage (DTU) analysis. The filter step will be used to select for genes and transcripts that satisfy rules for the number of samples in which a gene or transcript must be observed, and minimum threshold levels for the number of observed reads. The parameters used for filtering are `min_samps_gene_expr`, `min_samps_feature_expr`, `min_gene_expr`, and `min_feature_expr`. By default, any transcripts with zero expression or one transcript in all samples are filtered out at this stage.
#### 5.6 Differential transcript usage using DEXSeq
#### 4.6 Differential transcript usage using DEXSeq
Differential transcript usage analysis is performed using the R [DEXSeq](https://bioconductor.org/packages/release/bioc/html/DEXSeq.html) package ([Anders et al. (2012)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3460195/)). Similar to the edgeR package, DEXSeq estimates the variance between the biological replicates and applies generalised linear models for the statistical testing. The key difference is that the DEXSeq method looks for differences at the exon count level. DEXSeq uses the filtered transcript count data prepared earlier in this analysis.
#### 5.7 StageR stage-wise analysis of DGE and DTU
#### 4.7 StageR stage-wise analysis of DGE and DTU
The final component of this isoform analysis is a stage-wise statistical test using the R software package [stageR](https://bioconductor.org/packages/release/bioc/html/stageR.html)([Van den Berge and Clement (2018)](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1277-0)). stageR uses (1) the raw p-values for DTU from the DEXSeq analysis in the previous section and (2) a false-discovery corrected set of p-values from testing whether individual genes contain at least one exon showing DTU. A hierarchical two-stage statistical testing evaluates the set of genes for DTU.
## Troubleshooting
+ If the workflow fails please run it with the demo data set to ensure the workflow itself is working. This will help us determine if the issue is related to the environment, input parameters or a bug.

51
bin/workflow_glue/report.py
View File

@ -62,9 +62,6 @@ def argparser():
parser.add_argument(
"--isoform_table_nrows", required=False, type=int, default=5000,
help="Maximum rows to display in isoforms table")
parser.add_argument(
"--jaffal_csv", required=False, type=str, default=None,
help="Path to JAFFAL results csv")
parser.add_argument(
"--de_report", required=False, type=str, default=None,
help="Differential expression report optional")
@ -762,51 +759,6 @@ def seq_stats_tabs(report, stats):
section.plot(Tabs(tabs=[tabs.get(x) for x in sorted(tabs)]))
def jaffal_table(report, result_csv):
"""Make a table of fusion transcripts identified by JAFFAL."""
cols = [
'sample_id', 'fusion genes', 'chrom1', 'chrom2', 'spanning reads',
'classification', 'known']
section = report.add_section()
try:
df = pd.read_csv(result_csv)
except pd.errors.EmptyDataError:
section.markdown("""
### JAFFAL fusion transcript summary
This section summarizes putative fusion transcripts identified
by [JAFFAL](https://github.com/Oshlack/JAFFA/).
No fusion transcripts were detected for any of the samples.
""")
else:
sid_col = df.pop('sample_id')
df.insert(0, 'sample_id', sid_col)
df = df[cols]
df['chroms'] = df.chrom1.astype(str) + ':' + df.chrom2.astype(str)
df.rename(columns={
'spanning reads': 'nreads',
'fusion genes': 'genes'}, inplace=True)
df.drop(columns=['chrom1', 'chrom2'], inplace=True)
section.markdown("""
### JAFFAL fusion transcript summary
This table summarizes putative fusion transcripts identified
by [JAFFAL](https://github.com/Oshlack/JAFFA/).
* genes: the gene symbols of the fusion partners
* nreads: The number of reads supporting the fusion
* classification: JAFFAL's classification
* known: whether this fusion is in the given set of known gene fusions
* chroms: the respective, original chromosome location of the
two partner genes
""")
section.table(df)
def de_section(report):
"""Make differential transcript expression section."""
dexseq = os.path.join("de_report", "results_dexseq.tsv")
@ -874,9 +826,6 @@ def main(args):
if args.de_report:
de_section(report)
if args.jaffal_csv:
jaffal_table(report, args.jaffal_csv)
# Arguments and software versions
report.add_section(
section=scomponents.version_table(args.versions))

2
docs/01_brief_description.md
View File

@ -1 +1 @@
Transcriptome analysis including assembly and annotation of cDNA and direct RNA sequencing data, gene fusions and differential expression.
Transcriptome analysis including assembly and annotation of cDNA and direct RNA sequencing data and differential expression analysis.

1
docs/02_introduction.md
View File

@ -3,6 +3,5 @@ This workflow can be used for the following:
+ Identify RNA transcripts using either cDNA or direct RNA reads.
+ Reference aided transcriptome assembly.
+ Annotation of assembled transcripts.
+ Gene fusions detection.
+ Differential gene expression analysis using a pre-computed or assembled reference transcriptome.
+ Differential transcript usage analysis using a precomputed or assembled reference transcriptome.

3
docs/04_install_and_run.md
View File

@ -49,9 +49,6 @@ nextflow run epi2me-labs/wf-transcriptomes \
--de_analysis \
--direct_rna \
--fastq 'wf-transcriptomes-demo/differential_expression_fastq' \
--jaffal_annotation 'genCode22' \
--jaffal_genome 'hg38_chr20' \
--jaffal_refBase 'wf-transcriptomes-demo/chr20' \
--minimap2_index_opts '-k15' \
--ref_annotation 'wf-transcriptomes-demo/gencode.v22.annotation.chr20.gtf' \
--ref_genome 'wf-transcriptomes-demo/hg38_chr20.fa' \

11
docs/06_input_parameters.md
View File

@ -4,7 +4,7 @@
|--------------------------|------|-------------|------|---------|
| fastq | string | FASTQ files to use in the analysis. | This accepts one of three cases: (i) the path to a single FASTQ file; (ii) the path to a top-level directory containing FASTQ files; (iii) the path to a directory containing one level of sub-directories which in turn contain FASTQ files. In the first and second case, a sample name can be supplied with `--sample`. In the last case, the data is assumed to be multiplexed with the names of the sub-directories as barcodes. In this case, a sample sheet can be provided with `--sample_sheet`. | |
| bam | string | BAM or unaligned BAM (uBAM) files to use in the analysis. | This accepts one of three cases: (i) the path to a single BAM file; (ii) the path to a top-level directory containing BAM files; (iii) the path to a directory containing one level of sub-directories which in turn contain BAM files. In the first and second case, a sample name can be supplied with `--sample`. In the last case, the data is assumed to be multiplexed with the names of the sub-directories as barcodes. In this case, a sample sheet can be provided with `--sample_sheet`. | |
| transcriptome_source | string | Select how the transcriptome used for analysis should be prepared. | To analyse only gene fusions and differential expression use of an existing transcriptome may be preferred and so 'precomputed' should be selected. In this case the 'ref_transcriptome' parameter should be specified. To create a reference transcriptome using an existing reference genome, select 'reference guided' and specify the 'ref_genome' parameter. | reference-guided |
| transcriptome_source | string | Select how the transcriptome used for analysis should be prepared. | For differential expression analysis, use of an existing transcriptome may be preferred and so 'precomputed' should be selected. In this case the 'ref_transcriptome' parameter should be specified. To create a reference transcriptome using an existing reference genome, select 'reference guided' and specify the 'ref_genome' parameter. | reference-guided |
| ref_genome | string | Path to reference genome sequence [.fa/.fq/.fa.gz/fq.gz]. Required for reference-based workflow. | A reference genome is required for reference-based assembly of a transcriptome. | |
| ref_transcriptome | string | Transcriptome reference file. Required for precomputed transcriptome calculation and for differential expression analysis. | A reference transcriptome related to the sample under study. Must be supplied when the 'Transcriptome source' parameter has been set to 'precomputed' or to perform differential expression. | |
| ref_annotation | string | A reference annotation in GFF2 or GFF3 format (extensions .gtf(.gz), .gff(.gz), .gff3(.gz)). Only annotation files from [Encode](https://www.encodeproject.org), [Ensembl](https://www.ensembl.org/index.html) and [NCBI](https://www.ncbi.nlm.nih.gov/) are supported. | This will be used for guiding the transcriptome assembly and to label transcripts with their corresponding gene identifiers. Note: If in de_analysis mode transcript strands must be only + or -. | |
@ -39,15 +39,6 @@
| stringtie_opts | string | Extra command-line options for stringtie transcript assembly. | For additional String tie options see [here](https://github.com/gpertea/stringtie#stringtie-options). | --conservative |
### Gene Fusion Detection Options
| Nextflow parameter name | Type | Description | Help | Default |
|--------------------------|------|-------------|------|---------|
| jaffal_refBase | string | JAFFAl reference genome directory. | JAFFAL human hg38 reference data directory can be downloaded from here: https://figshare.com/ndownloader/files/25410494 or see the README for alternative instructions. If custom gemome files are required, see the instructions here: https://github.com/Oshlack/JAFFA/wiki/FAQandTroubleshooting#how-can-i-generate-the-reference-files-for-a-non-supported-genome. | |
| jaffal_genome | string | Genome reference prefix. e.g. hg38. | JAFFAL reference files are prefixed with the genome reference file name and need to be supplied . If using the human reference data provided by JAFFAL, this can be left at `hg38`. | hg38 |
| jaffal_annotation | string | Annotation suffix. | JAFFAL reference files are suffixed with the annotation filename and this needs to be supplied. For the human hg38 reference data supplied by JAFFAL, this is `genCode22`. | genCode22 |
### Differential Expression Options
| Nextflow parameter name | Type | Description | Help | Default |

2
docs/07_outputs.md
View File

@ -26,5 +26,3 @@ Output files may be aggregated including information for all samples or provided
| Transcript counts filtered | de_analysis/filtered_transcript_counts_with_genes.tsv | Filtered transcript counts, used for differential transcript usage analysis. Includes a reference to the associated gene ID. | aggregated |
| Transcript info table | {{ alias }}_transcripts_table.tsv | This file details each isoform that was reconstructed from the input reads. It contains a subset of columns from the .tmap output from [gffcompare](https://ccb.jhu.edu/software/stringtie/gffcompare.shtml) | per-sample |
| Final non redundant transcriptome | de_analysis/final_non_redundant_transcriptome.fasta | Transcripts that were used for differential expression analysis including novel transcripts with the identifiers used for DE analysis. | aggregated |
| Fusion transcript sequences | jaffal_output_{{ alias }}/jaffa_results.fasta | Fusion transcript sequences output by Jaffa. | per-sample |
| Fusion transcript sequence summary file | jaffal_output_{{ alias }}/jaffa_results.csv | Fusion transcript sequences summary file output by Jaffa. | per-sample |

22
docs/08_pipeline_overview.md
View File

@ -25,10 +25,7 @@ Transcript GFF files from the chunks with the same sample aliases will then be m
#### 3.6 Create transcriptomes
[Gffread](https://github.com/gpertea/gffread) is used to create a transcriptome FASTA file from the final GFF as well as a merged transcriptome that includes annotations in the FASTA headers where available.
### 4. Find gene fusions
If gene fusion options are provided, fusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFFA), with the JAFFAL extension. To enable this provide the gene fusion detection options: `jaffal_refBase`, `jaffal_genome` and `jaffal_annotation`.
### 5. Differential expression analysis
### 4. Differential expression analysis
Differential gene expression (DGE) and differential transcript usage (DTU) analyses aim to identify genes and transcripts that show statistically altered expression patterns.
@ -51,26 +48,23 @@ barcode05,sample05,treated
barcode06,sample06,treated
```
#### 5.1 Merge cross sample transcriptomes
#### 4.1 Merge cross sample transcriptomes
If a `ref_transcriptome` is not provided, the transcriptomes created by the workflow will be used for DE analysis. To do this, the GFF outputs of GffCompare are merged using StringTie. A final non redundant FASTA file of the transcripts is created using the merged GFF file and the reference genome using seqkit.
#### 5.2 Create a final non redundant transcriptome
#### 4.2 Create a final non redundant transcriptome
The reads from all the samples will be aligned with the final non redundant transcriptome using Minimap2 in a splice aware manner.
#### 5.3 Count genes and transcripts
#### 4.3 Count genes and transcripts
[Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving gene and transcript counts.
#### 5.4 edgeR based differential expression analysis
#### 4.4 edgeR based differential expression analysis
A statistical analysis is first performed using [edgeR](https://bioconductor.org/packages/release/bioc/html/edgeR.html) to identify the subset of differentially expressed genes using the gene counts as input. A normalisation factor is calculated for each sequence library using the default TMM method (see [McCarthy et al. (2012)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3378882/) for further details). The defined experimental design is used to calculate estimates of dispersion for each of the gene features. Statistical tests are calculated using the contrasts defined in the experimental design. The differentially expressed genes are corrected for false discovery (FDR) using the method of Benjamini & Hochberg ([Benjamini and Hochberg (1995)](https://www.jstor.org/stable/2346101))
#### 5.5 Pre-filtering of quantitative data using DRIMSeq
#### 4.5 Pre-filtering of quantitative data using DRIMSeq
[DRIMSeq](https://bioconductor.org/packages/release/bioc/html/DRIMSeq.html) is used to filter the transcript count data from the Salmon analysis for differential transcript usage (DTU) analysis. The filter step will be used to select for genes and transcripts that satisfy rules for the number of samples in which a gene or transcript must be observed, and minimum threshold levels for the number of observed reads. The parameters used for filtering are `min_samps_gene_expr`, `min_samps_feature_expr`, `min_gene_expr`, and `min_feature_expr`. By default, any transcripts with zero expression or one transcript in all samples are filtered out at this stage.
#### 5.6 Differential transcript usage using DEXSeq
#### 4.6 Differential transcript usage using DEXSeq
Differential transcript usage analysis is performed using the R [DEXSeq](https://bioconductor.org/packages/release/bioc/html/DEXSeq.html) package ([Anders et al. (2012)](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3460195/)). Similar to the edgeR package, DEXSeq estimates the variance between the biological replicates and applies generalised linear models for the statistical testing. The key difference is that the DEXSeq method looks for differences at the exon count level. DEXSeq uses the filtered transcript count data prepared earlier in this analysis.
#### 5.7 StageR stage-wise analysis of DGE and DTU
#### 4.7 StageR stage-wise analysis of DGE and DTU
The final component of this isoform analysis is a stage-wise statistical test using the R software package [stageR](https://bioconductor.org/packages/release/bioc/html/stageR.html)([Van den Berge and Clement (2018)](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1277-0)). stageR uses (1) the raw p-values for DTU from the DEXSeq analysis in the previous section and (2) a false-discovery corrected set of p-values from testing whether individual genes contain at least one exon showing DTU. A hierarchical two-stage statistical testing evaluates the set of genes for DTU.

49
main.nf
View File

@ -11,7 +11,6 @@ nextflow.enable.dsl = 2
include { fastq_ingress; xam_ingress } from './lib/ingress'
include { reference_assembly } from './subworkflows/reference_assembly'
include { gene_fusions } from './subworkflows/JAFFAL/gene_fusions'
include { differential_expression } from './subworkflows/differential_expression'
OPTIONAL_FILE = file("$projectDir/data/OPTIONAL_FILE")
@ -395,7 +394,6 @@ process makeReport {
path versions
path "params.json"
path "pychopper_report/*"
path"jaffal_csv/*"
path "per_read_stats/?.gz"
path "aln_stats/*"
path "gffcmp_dir/*"
@ -430,11 +428,6 @@ process makeReport {
else
OPT_GFFCMP_DIR="--gffcompare_dir gffcmp_dir/"
fi
if [ -f "jaffal_csv/OPTIONAL_FILE" ]; then
OPT_JAFFAL_CSV=""
else
OPT_JAFFAL_CSV="--jaffal_csv jaffal_csv/*"
fi
if [ -f "aln_stats/OPTIONAL_FILE" ]; then
OPT_ALN=""
else
@ -460,7 +453,6 @@ process makeReport {
${OPT_ISO_TABLE} \
${OPT_GFFCMP_DIR} \
--isoform_table_nrows !{params.isoform_table_nrows} \
${OPT_JAFFAL_CSV} \
${dereport}
'''
}
@ -500,7 +492,7 @@ process collectFastqIngressResultsInDir {
// publish files from a workflow whilst decoupling the publish from the process steps.
// The process takes a tuple containing the filename and the name of a sub-directory to
// put the file into. If the latter is `null`, puts it into the top-level directory.
process output {
process publish_results {
// publish inputs to output directory
label "isoforms"
cpus 1
@ -541,9 +533,6 @@ workflow pipeline {
reads
ref_genome
ref_annotation
jaffal_refBase
jaffal_genome
jaffal_annotation
ref_transcriptome
use_ref_ann
main:
@ -655,13 +644,6 @@ workflow pipeline {
use_ref_ann = false
}
if (jaffal_refBase){
gene_fusions(full_len_reads, jaffal_refBase, jaffal_genome, jaffal_annotation)
jaffal_out = gene_fusions.out.results_csv.map{ alias, csv -> csv}.collectFile(keepHeader: true, name: 'jaffal.csv')
} else{
jaffal_out = OPTIONAL_FILE
}
if (params.de_analysis){
sample_sheet = file(params.sample_sheet, type:"file")
// check ref annotation contains only + or - strand as DE analysis will error on .
@ -698,7 +680,6 @@ workflow pipeline {
software_versions,
workflow_params,
pychopper_report,
jaffal_out,
per_read_stats,
assembly_stats,
gff_compare,
@ -728,14 +709,9 @@ workflow pipeline {
.concat(results)
}
if (params.jaffal_refBase){
results = results
.concat(gene_fusions.out.results
.map {it -> it[1]})
}
results = results.map{ [it, null] }.concat(fastq_ingress_results.map { [it, "fastq_ingress_results"] })
if (params.de_analysis){
de_results = report.concat(
transcriptome, de_outputs.flatten(), counts.flatten(),
@ -761,6 +737,9 @@ workflow {
error = null
if (params.containsValue("jaffal_refBase")) {
error = "JAFFAL fusion detection has been removed from this workflow."
}
if (params.containsKey("minimap_index_opts")) {
error = "`--minimap_index_opts` parameter is deprecated. Use parameter `--minimap2_index_opts` instead."
}
@ -794,14 +773,6 @@ workflow {
ref_annotation= OPTIONAL_FILE
use_ref_ann = false
}
if (params.jaffal_refBase){
jaffal_refBase = file(params.jaffal_refBase, type: "dir")
if (!jaffal_refBase.exists()) {
error = "--jaffa_refBase: Directory doesn't exist, check path."
}
}else{
jaffal_refBase = null
}
ref_transcriptome = OPTIONAL_FILE
if (params.ref_transcriptome){
log.info("Reference Transcriptome provided will be used for differential expression.")
@ -847,12 +818,8 @@ workflow {
"per_read_stats": true])
}
pipeline(samples, ref_genome, ref_annotation,
jaffal_refBase, params.jaffal_genome, params.jaffal_annotation,
ref_transcriptome, use_ref_ann)
output(pipeline.out.results)
pipeline(samples, ref_genome, ref_annotation, ref_transcriptome, use_ref_ann)
publish_results(pipeline.out.results)
}
workflow.onComplete {

17
nextflow.config
View File

@ -74,14 +74,6 @@ params {
disable_ping = false
////// Fusion detection parameters
jaffal_refBase = null
jaffal_genome = "hg38"
jaffal_annotation = "genCode22"
// The default location of the JAFFA src directory when running in EPI2ME-Labs environment
// This needs overriding if running elsewhere
jaffal_dir = "/home/epi2melabs/JAFFA"
// de options
de_analysis = false
ref_transcriptome = null
@ -96,16 +88,13 @@ params {
"--de_analysis",
"--direct_rna",
"--fastq 'wf-transcriptomes-demo/differential_expression_fastq'",
"--jaffal_annotation 'genCode22'",
"--jaffal_genome 'hg38_chr20'",
"--jaffal_refBase 'wf-transcriptomes-demo/chr20'",
"--minimap2_index_opts '-k15'",
"--ref_annotation 'wf-transcriptomes-demo/gencode.v22.annotation.chr20.gtf'",
"--ref_genome 'wf-transcriptomes-demo/hg38_chr20.fa'",
"--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv'",
]
agent = null
container_sha = "shae7c9f184996a384e99be68e790f0612f0c732867"
container_sha = "shafb1e2372e1535f0b42891ed2c68ffdac2ca1d658"
common_sha = "sha8b5843d549bb210558cbb676fe537a153ce771d6"
}
}
@ -114,7 +103,7 @@ manifest {
name = 'epi2me-labs/wf-transcriptomes'
author = 'Oxford Nanopore Technologies'
homePage = 'https://github.com/epi2me-labs/wf-transcriptomes'
description = 'Transcriptome analysis including gene fusions, differential expression as well as assembly and annotation of cDNA and direct RNA sequencing data.'
description = 'Transcriptome analysis including differential expression as well as assembly and annotation of cDNA and direct RNA sequencing data.'
mainScript = 'main.nf'
nextflowVersion = '>=23.04.2'
version = 'v1.2.1'
@ -171,7 +160,7 @@ profiles {
queue = "${params.aws_queue}"
memory = '8G'
withLabel:isoforms {
container = "${params.aws_image_prefix}-wf-transcriptomes:${params.wf.container_sha}-root"
container = "${params.aws_image_prefix}-wf-transcriptomes:${params.wf.container_sha}"
}
withLabel:wf_common {
container = "${params.aws_image_prefix}-wf-common:${params.wf.common_sha}"

42
nextflow_schema.json
View File

@ -3,7 +3,7 @@
"$id": "https://raw.githubusercontent.com/./master/nextflow_schema.json",
"title": "epi2me-labs/wf-transcriptomes",
"workflow_title": "Workflow Transcriptomes",
"description": "Transcriptome analysis including assembly and annotation of cDNA and direct RNA sequencing data, gene fusions and differential expression.",
"description": "Transcriptome analysis including assembly and annotation of cDNA and direct RNA sequencing data and differential expression analysis.",
"demo_url": "https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/wf-transcriptomes-demo.tar.gz",
"aws_demo_url": "https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-transcriptomes/wf-transcriptomes-demo/aws.nextflow.config",
"url": "https://github.com/epi2me-labs/wf-transcriptomes",
@ -37,7 +37,7 @@
],
"default": "reference-guided",
"description": "Select how the transcriptome used for analysis should be prepared.",
"help_text": "To analyse only gene fusions and differential expression use of an existing transcriptome may be preferred and so 'precomputed' should be selected. In this case the 'ref_transcriptome' parameter should be specified. To create a reference transcriptome using an existing reference genome, select 'reference guided' and specify the 'ref_genome' parameter."
"help_text": "For differential expression analysis, use of an existing transcriptome may be preferred and so 'precomputed' should be selected. In this case the 'ref_transcriptome' parameter should be specified. To create a reference transcriptome using an existing reference genome, select 'reference guided' and specify the 'ref_genome' parameter."
},
"ref_genome": {
"type": "string",
@ -195,41 +195,6 @@
}
}
},
"gene_fusion_detection_options": {
"title": "Gene Fusion Detection Options",
"type": "object",
"description": "Parameters for gene fusion detection",
"properties": {
"jaffal_refBase": {
"type": "string",
"title": "JAFFAL reference genome directory",
"format": "directory-path",
"description": "JAFFAl reference genome directory.",
"help_text": "JAFFAL human hg38 reference data directory can be downloaded from here: https://figshare.com/ndownloader/files/25410494 or see the README for alternative instructions. If custom gemome files are required, see the instructions here: https://github.com/Oshlack/JAFFA/wiki/FAQandTroubleshooting#how-can-i-generate-the-reference-files-for-a-non-supported-genome."
},
"jaffal_genome": {
"type": "string",
"title": "JAFFAL genome reference prefix",
"description": "Genome reference prefix. e.g. hg38.",
"help_text": "JAFFAL reference files are prefixed with the genome reference file name and need to be supplied . If using the human reference data provided by JAFFAL, this can be left at `hg38`.",
"default": "hg38"
},
"jaffal_annotation": {
"type": "string",
"title": "JAFFAL annotation suffix",
"description": "Annotation suffix.",
"help_text": "JAFFAL reference files are suffixed with the annotation filename and this needs to be supplied. For the human hg38 reference data supplied by JAFFAL, this is `genCode22`.",
"default": "genCode22"
},
"jaffal_dir": {
"type": "string",
"format": "directory-path",
"description": "Path to the JAFFAL code directory. If running within EPI2ME-Labs, the default path of /home/epi2melabs/JAFFA within the application container will be used. If using outside of EPI2ME-Labs, the path to the code directory downloaded from github should be supplied.",
"default": "/home/epi2melabs/JAFFA",
"hidden": true
}
}
},
"differential_expression_options": {
"title": "Differential Expression Options",
"type": "object",
@ -368,9 +333,6 @@
{
"$ref": "#/definitions/options_for_reference_based_workflow"
},
{
"$ref": "#/definitions/gene_fusion_detection_options"
},
{
"$ref": "#/definitions/differential_expression_options"
},

17
output_definition.json
View File

@ -191,23 +191,6 @@
"mime-type": "text/x-fasta",
"optional": true,
"type": "aggregated"
},
"jaffa_fasta": {
"filepath": "jaffal_output_{{ alias }}/jaffa_results.fasta",
"title": "Fusion transcript sequences",
"description": "Fusion transcript sequences output by Jaffa.",
"mime-type": "text/x-fasta",
"optional": true,
"type": "per-sample"
},
"jaffa_results": {
"filepath": "jaffal_output_{{ alias }}/jaffa_results.csv",
"title": "Fusion transcript sequence summary file",
"description": "Fusion transcript sequences summary file output by Jaffa.",
"mime-type": "text/csv",
"optional": true,
"type": "per-sample"
}
}
}

2
subworkflows/JAFFAL/build_jaffal_ref.sh
View File

@ -1,2 +0,0 @@
#!/bin/sh

5
subworkflows/JAFFAL/download_jaffal_references.sh
View File

@ -1,5 +0,0 @@
#!/bin/sh
#Download the data. We should we move the data out of Figshare?
wget -O JAFFA_REFERENCE_FILES_HG38_GENCODE22.V2.tar.gz https://figshare.com/ndownloader/files/25410494
tar -zxvf JAFFA_REFERENCE_FILES_HG38_GENCODE22.V2.tar.gz

66
subworkflows/JAFFAL/gene_fusions.nf
View File

@ -1,66 +0,0 @@
process jaffal{
label "isoforms"
cpus params.threads
memory "31 GB"
input:
tuple val(sample_id), path(fastq)
path refBase
val genome
val annotation
output:
tuple val(sample_id), path("jaffal_output_$sample_id"), emit: results
tuple val(sample_id), path("jaffal_output_$sample_id/*jaffa_results.csv"), emit: results_csv
script:
"""
# Jaffa requires a virtual environment
# so override JAVA_TOOL_OPTIONS env variable.
JAVA_TOOL_OPTIONS=""
JAFFAOUT=jaffal_output_$sample_id
# JAFFAL exists with status code 1 when there's 0 fusion hits. Prevent this with '||:'
$params.jaffal_dir/tools/bin/bpipe run \
-n "${task.cpus}" \
-p jaffa_output="\$JAFFAOUT/" \
-p refBase=$refBase \
-p genome=$genome \
-p annotation=$annotation \
-p fastqInputFormat="*.fastq" \
$params.jaffal_dir/JAFFAL.groovy \
$fastq || :
summary="\$JAFFAOUT/all/all.summary"
if [ -f \$summary ]; then
# The summary is writtten so assume JAFFAL completed.
if [ ! -s \$summary ]; then
echo "JAFFAL failed to find any fusion transcripts for ${sample_id}"
touch "\$JAFFAOUT/${sample_id}_jaffa_results.csv"
else
echo JAFFAL found fusion transcripts for ${sample_id}
mv "\$JAFFAOUT/jaffa_results.csv" "\$JAFFAOUT/${sample_id}_jaffa_results.csv"
# Add sample id column and header
sed "s/\$/,${sample_id}/" \$JAFFAOUT/${sample_id}_jaffa_results.csv \
| sed "1 s/${sample_id}/sample_id/" > tmp
mv tmp \$JAFFAOUT/${sample_id}_jaffa_results.csv
fi
else
echo JAFFAL encountered an error while prosessing ${sample_id}
fi
"""
}
// workflow module
workflow gene_fusions {
take:
fastq
refBase
genome
annotation
main:
jaffal(fastq, refBase, genome, annotation)
emit:
results_csv = jaffal.out.results_csv
results = jaffal.out.results
}

21
subworkflows/JAFFAL/install_jaffa.sh
View File

@ -1,21 +0,0 @@
#!/bin/bash
set -e
git clone https://github.com/Oshlack/JAFFA.git &&
cd JAFFA
git checkout 24b1c3b
cp ../subworkflows/JAFFAL/install_linux64.sh .
./install_linux64.sh
# JAFFA uses tools.groovy to locate binaries. We modify it as we only need a small subset or they are already
# included in our env
# Tools to be compiled from src/
bin=$(realpath tools/bin)
echo "//Tools built locally" >> tools.groovy
declare -a tools=("reformat" "extract_seq_from_fasta" "make_simple_read_table" "process_transcriptome_align_table" "make_3_gene_fusion_table")
for b in "${tools[@]}"; do
echo "$b=\"$bin/$b\"" >> tools.groovy
done
echo "minimap2=\"minimap2\"" >> tools.groovy

111
subworkflows/JAFFAL/install_linux64.sh
View File

@ -1,111 +0,0 @@
#!/bin/bash
# 21/06/22: This script has been modified to install only those applications needed for epi2melabs/wf-transcriptomes
## This script will install the tools required for the JAFFA pipeline.
## It will fetched each tool from the web and placed into the tools/ subdirectory.
## Paths to all installed tools can be found in the file tools.groovy at the
## end of execution of this script. These paths can be changed if a different
## version of software is required. Note that R must be installed manually
##
## Last Modified: Sep. 2021 by Nadia Davidson
mkdir -p tools/bin
cd tools
#a list of which programs need to be installed
commands="bpipe reformat extract_seq_from_fasta make_simple_read_table process_transcriptome_align_table make_3_gene_fusion_table dedupe"
#installation methods
function bpipe_install {
wget -O bpipe-0.9.9.2.tar.gz https://github.com/ssadedin/bpipe/releases/download/0.9.9.2/bpipe-0.9.9.2.tar.gz
tar -zxvf bpipe-0.9.9.2.tar.gz ; rm bpipe-0.9.9.2.tar.gz
ln -s $PWD/bpipe-0.9.9.2/bin/* $PWD/bin/
}
function make_3_gene_fusion_table_install {
g++ -std=c++11 -O3 -o bin/make_3_gene_fusion_table ../src/make_3_gene_fusion_table.c++
}
function extract_seq_from_fasta_install {
g++ -std=c++11 -O3 -o bin/extract_seq_from_fasta ../src/extract_seq_from_fasta.c++
}
function make_simple_read_table_install {
g++ -std=c++11 -O3 -o bin/make_simple_read_table ../src/make_simple_read_table.c++
}
function process_transcriptome_align_table_install {
g++ -std=c++11 -O3 -o bin/process_transcriptome_align_table ../src/process_transcriptome_align_table.c++
}
function make_count_table_install {
g++ -O3 -o bin/make_count_table ../src/make_count_table.c++
}
function dedupe_install {
wget --no-check-certificate https://sourceforge.net/projects/bbmap/files/BBMap_36.59.tar.gz
tar -zxvf BBMap_36.59.tar.gz
rm BBMap_36.59.tar.gz
for script in `ls $PWD/bbmap/*.sh` ; do
s=`basename $script`
s_pre=`echo $s | sed 's/.sh//g'`
echo "$PWD/bbmap/$s \$@" > $PWD/bin/$s_pre
chmod +x $PWD/bin/$s_pre
done
}
#function bypass_genomic_alignment_install {
# g++ -std=c++11 -O3 -o bin/bypass_genomic_alignment ../src/bypass_genomic_alignment.c++
#}
#Check if the version of gcc is >= 4.9
gcc_version=`gcc -dumpversion`
gcc_check=`echo -e "$gcc_version\n4.9" | sort -n | tail -n1`
if [[ $gcc_chek = "4.9" ]]
then
echo "Your version of gcc is $gcc_version."
echo "gcc must be >= 4.9 to install JAFFA. Exiting..."
exit 1
fi
echo "gcc check passed"
echo "// Path to tools used by the JAFFA pipeline" > ../tools.groovy
for c in $commands ; do
c_path=`which $PWD/bin/$c 2>/dev/null`
if [ -z $c_path ] ; then
echo "$c not found, fetching it"
${c}_install
c_path=`which $PWD/bin/$c 2>/dev/null`
fi
echo "$c=\"$c_path\"" >> ../tools.groovy
done
#finally check that R is install
R_path=`which R 2>/dev/null`
if [ -z $R_path ] ; then
echo "R not found!"
echo "Please go to http://www.r-project.org/ and follow the installation instructions."
echo "Note that the IRanges R package must be installed."
fi
echo "R=\"$R_path\"" >> ../tools.groovy
#loop through commands to check they are all installed
echo "Checking that all required tools were installed:"
Final_message="All commands installed successfully!"
for c in $commands ; do
c_path=`which $PWD/bin/$c 2>/dev/null`
if [ -z $c_path ] ; then
echo -n "WARNING: $c could not be found!!!! "
echo "You will need to download and install $c manually, then add its path to tools.groovy"
Final_message="WARNING: One or more command did not install successfully. See warning messages above. \
You will need to correct this before running JAFFA."
else
echo "$c looks like it has been installed"
fi
done
echo "**********************************************************"
echo $Final_message

10
subworkflows/JAFFAL/tools.groovy
View File

@ -1,10 +0,0 @@
// Path to tools used by the JAFFA pipeline
// Conda-installable tools
bpipe="bpipe"
//trimmomatic="trimmomatic"
R="/usr/bin/R"
minimap2="minimap2"
dedupe="dedupe"