From ec3c346fd8ea83c4b969cea9abf25e03c03c8157 Mon Sep 17 00:00:00 2001 From: Sarah Griffiths Date: Wed, 26 Jul 2023 18:14:00 +0000 Subject: [PATCH] fix and test for files with no gene id in attributes --- .gitlab-ci.yml | 36 +++++++++++++++++++--- CHANGELOG.md | 5 ++- README.md | 2 +- bin/de_analysis.R | 12 +++----- bin/workflow_glue/de_plots.py | 41 +++++++++++++++++-------- docs/intro.md | 2 +- nextflow.config | 2 +- nextflow_schema.json | 2 +- subworkflows/differential_expression.nf | 18 +++++++++-- 9 files changed, 86 insertions(+), 34 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 4915790..dd1b9d1 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -45,7 +45,7 @@ docker-run: - MATRIX_NAME: [ "fusions", "differential_expression", "isoforms", "only_differential_expression", "differential_expression_gff3", - "ncbi_gzip", "denovo" + "ncbi_gzip", "denovo", "ncbi_no_gene_id", "ensembl_with_versions" ] rules: # NOTE As we're overriding the rules block for the included docker-run @@ -109,15 +109,41 @@ docker-run: build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam - if: $MATRIX_NAME == "ncbi_gzip" variables: - NF_BEFORE_SCRIPT: wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf differential_expression.tar.gz + NF_BEFORE_SCRIPT: wget -O differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf differential_expression_ncbi.tar.gz NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB \ ---fastq differential_expression/differential_expression_fastq \ + --fastq differential_expression_ncbi/differential_expression_fastq \ --transcriptome-source precomputed \ --de_analysis \ - --ref_genome differential_expression/GRCh38.p14.NCBI_test.fna.gz \ - --ref_annotation differential_expression/GRCh38.p14_NCBI_test.gtf.gz \ + --ref_genome differential_expression_ncbi/GRCh38.p14.NCBI_test.fna.gz \ + --ref_annotation differential_expression_ncbi/GRCh38.p14_NCBI_test.gtf.gz \ --direct_rna --minimap_index_opts '-w 25' \ --transcriptome_assembly false --sample_sheet test_data/sample_sheet.csv" NF_IGNORE_PROCESSES: > preprocess_reads,merge_transcriptomes,assemble_transcripts, build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam + - if: $MATRIX_NAME == "ncbi_no_gene_id" + variables: + NF_BEFORE_SCRIPT: wget -O differential_expression_ncbi.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression_ncbi.tar.gz && tar -xzvf differential_expression_ncbi.tar.gz + NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB \ + --fastq differential_expression_ncbi/differential_expression_fastq \ + --transcriptome-source precomputed --de_analysis \ + --ref_genome differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna.gz \ + --ref_annotation differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.gff.gz \ + --direct_rna --ref_transcriptome differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_rna.fna.gz \ + --transcriptome_assembly false --sample_sheet test_data/sample_sheet.csv" + NF_IGNORE_PROCESSES: > + preprocess_reads,merge_transcriptomes,assemble_transcripts, + build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam + - if: $MATRIX_NAME == "ensembl_with_versions" + variables: + NF_BEFORE_SCRIPT: wget -O differential_expression.tar.gz https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-isoforms/differential_expression.tar.gz && tar -xzvf differential_expression.tar.gz + NF_WORKFLOW_OPTS: "-executor.\\$$local.memory 16GB \ + --fastq differential_expression/differential_expression_fastq \ + --transcriptome-source precomputed --de_analysis \ + --ref_genome differential_expression/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz \ + --ref_annotation differential_expression/Homo_sapiens.GRCh38.109.gtf.gz \ + --direct_rna --ref_transcriptome differential_expression/Homo_sapiens.GRCh38.cdna.all.fa.gz \ + --transcriptome_assembly false --sample_sheet test_data/sample_sheet.csv" + NF_IGNORE_PROCESSES: > + preprocess_reads,merge_transcriptomes,assemble_transcripts, + build_minimap_index,get_transcriptome,merge_gff_bundles,run_gffcompare,build_minimap_index,split_bam diff --git a/CHANGELOG.md b/CHANGELOG.md index 8225d14..9f140bd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,12 +4,15 @@ 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] +## [v0.2.1] ### Changed - Any sample aliases that contain spaces will be replaced with underscores. +- Updated documentation to explain we only support Ensembl, NCBI and ENCODE annotation file types. ### Fixed - Documentation parameter examples corrected. +- Handling for annotation files that use gene as gene_id attribute. +- Handling for Ensembl annotation files. ## [v0.2.0] ### Changed diff --git a/README.md b/README.md index 6a2ec58..1925c14 100644 --- a/README.md +++ b/README.md @@ -82,7 +82,7 @@ Differential gene expression is sensitive to the input data quantity and quality - Directory containing cDNA/direct RNA reads. Or a directory containing subdirectories each with reads from different samples (in fastq/fastq.gz format) - Reference genome in fasta format (required for reference-based assembly). -- Optional reference annotation in GFF2/3 format (extensions allowed are .gtf(.gz), .gff(.gz), .gff3(.gz)) (required for differential expression analysis `--de_analysis`). +- Optional reference annotation in GFF2/3 format (extensions allowed are .gtf(.gz), .gff(.gz), .gff3(.gz)) (required for differential expression analysis `--de_analysis`). 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. - For fusion detection, JAFFAL reference files (see Quickstart) diff --git a/bin/de_analysis.R b/bin/de_analysis.R index 5dceb26..f60dcc8 100755 --- a/bin/de_analysis.R +++ b/bin/de_analysis.R @@ -9,6 +9,7 @@ min_samps_feature_expr <- args[3] min_gene_expr <- args[4] min_feature_expr <- args[5] annotation_type <- args[6] +strip_version <- args[7] cat("Loading counts, conditions and parameters.\n") cts <- as.matrix(read.csv("merged/all_counts.tsv", sep="\t", row.names="Reference", stringsAsFactors=FALSE)) @@ -27,15 +28,10 @@ txdf <- select(txdb, keys(txdb,"GENEID"), "TXNAME", "GENEID") tab <- table(txdf$GENEID) txdf$ntx<- tab[match(txdf$GENEID, names(tab))] -strip_version<-function(x) { - tmp<-data.frame(strsplit(x,".", fixed=TRUE), stringsAsFactors=FALSE) - tmp<-as.vector(tmp[1,]) - colnames(tmp) <- c() - rownames(tmp) <- c() - return(tmp) -} -#rownames(cts) <- strip_version(rownames(cts)) +if (strip_version == "true"){ + rownames(cts) <- lapply(rownames(cts), sub, pattern = "\\.\\d+$", replacement = "") +} cts <- cts[rownames(cts) %in% txdf$TXNAME, ] # FIXME: filter for transcripts which are in the annotation. Why they are not all there? diff --git a/bin/workflow_glue/de_plots.py b/bin/workflow_glue/de_plots.py index 98174ca..b728686 100755 --- a/bin/workflow_glue/de_plots.py +++ b/bin/workflow_glue/de_plots.py @@ -173,8 +173,9 @@ def dexseq_section(dexseq_file, section, id_dic): section.markdown(dexseq_caption) dexseq_results = pd.read_csv(dexseq_file, sep='\t') dexseq_results.index.name = "gene_id:trancript_id" + # Replace gene id with more useful gene name where possible dexseq_results.index = dexseq_results.index.map( - lambda x: id_dic[x.split(':')[0]] + ':' + str(x.split(':')[1])) + lambda x: str(id_dic.get(x.split(':')[0])) + ':' + str(x.split(':')[1])) dexseq_pvals = dexseq_results.sort_values(by='pvalue', ascending=True) section.table(dexseq_results.loc[dexseq_pvals.index], index=True) section.markdown(""" @@ -220,8 +221,10 @@ thresholds defined are shaded as 'Up-' or 'Down-' regulated.""") def dtu_section(dtu_file, section, gt_dic, ge_dic): """Plot dtu section.""" dtu_results = pd.read_csv(dtu_file, sep='\t') - dtu_results["gene_name"] = dtu_results["txID"].apply(lambda x: gt_dic[x]) - dtu_results["geneID"] = dtu_results["geneID"].apply(lambda x: ge_dic[x]) + dtu_results["gene_name"] = dtu_results["txID"].apply( + lambda x: gt_dic.get(x)) + dtu_results["geneID"] = dtu_results["geneID"].apply( + lambda x: ge_dic.get(x)) dtu_pvals = dtu_results.sort_values(by='gene', ascending=True) dtu_caption = '''Table showing gene and transcript identifiers and their FDR corrected probabilities @@ -248,8 +251,8 @@ experimental conditions, the log-scaled counts per million measure of abundance and the false discovery corrected p-value (FDR). This table has not been filtered for genes that satisfy statistical or magnitudinal thresholds""" section.markdown(dge_caption) - dge_results.index = dge_results.index.map(lambda x: ids_dic[x]) - dge_pvals.index = dge_pvals.index.map(lambda x: ids_dic[x]) + dge_results.index = dge_results.index.map(lambda x: ids_dic.get(x)) + dge_pvals.index = dge_pvals.index.map(lambda x: ids_dic.get(x)) section.table(dge_results.loc[dge_pvals.index], index=True) dge = pd.read_csv(dge_file, sep="\t") section.markdown(""" @@ -317,19 +320,31 @@ def get_translations(gtf): for i in fn: if i.startswith("#"): continue - try: - gene_name = get_feature(i, 'gene_name') - except IndexError: + # Different gtf/gff formats contain different attributes + # and different formating (eg. gene_name="xyz" or gene_name "xyz") + if 'gene_name' in i: + gene_name = get_feature(i, "gene_name") + elif 'gene_id' in i: gene_name = get_feature(i, 'gene_id') - try: + elif 'gene' in i: + gene_name = get_feature(i, "gene") + else: + continue + + if 'ref_gene_id' in i: gene_reference = get_feature(i, 'ref_gene_id') - except IndexError: + elif 'gene_id' in i: gene_reference = get_feature(i, 'gene_id') - try: + else: + gene_reference = gene_name + if 'transcript_id' in i: transcript_id = get_feature(i, 'transcript_id') - except IndexError: + else: transcript_id = "unknown" - gene_id = get_feature(i, 'gene_id') + if 'gene_id' in i: + gene_id = get_feature(i, 'gene_id') + else: + gene_id = gene_name gene_txid[transcript_id] = gene_name gene_geid[gene_id] = gene_reference return gene_txid, gene_geid diff --git a/docs/intro.md b/docs/intro.md index a7c0fa3..605d676 100644 --- a/docs/intro.md +++ b/docs/intro.md @@ -70,5 +70,5 @@ Differential gene expression is sensitive to the input data quantity and quality - Directory containing cDNA/direct RNA reads. Or a directory containing subdirectories each with reads from different samples (in fastq/fastq.gz format) - Reference genome in fasta format (required for reference-based assembly). -- Optional reference annotation in GFF2/3 format (extensions allowed are .gtf(.gz), .gff(.gz), .gff3(.gz)) (required for differential expression analysis `--de_analysis`). +- Optional reference annotation in GFF2/3 format (extensions allowed are .gtf(.gz), .gff(.gz), .gff3(.gz)) (required for differential expression analysis `--de_analysis`). 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. - For fusion detection, JAFFAL reference files (see Quickstart) diff --git a/nextflow.config b/nextflow.config index 0d309d2..677ed0f 100644 --- a/nextflow.config +++ b/nextflow.config @@ -115,7 +115,7 @@ manifest { description = 'Transcriptome analysis including gene fusions, differential expression as well as assembly and annotation of cDNA and direct RNA sequencing data.' mainScript = 'main.nf' nextflowVersion = '>=22.10.8' - version = 'v0.2.0' + version = 'v0.2.1' } executor { diff --git a/nextflow_schema.json b/nextflow_schema.json index 5654cf0..e470723 100644 --- a/nextflow_schema.json +++ b/nextflow_schema.json @@ -345,7 +345,7 @@ } }, "docs": { - "intro": "## Introduction\n\nThis workflow identifies RNA isoforms using either cDNA or direct RNA (dRNA) \nOxford Nanopore reads.\n\n### Preprocesing\ncDNA reads are initially preprocessed by [pychopper](https://github.com/epi2me-labs/pychopper) \nfor the identification of full-length reads, as well as trimming and orientation correction (This step is omitted for \n direct RNA reads).\n\n\n### Transcript assembly\n\n#### Reference-aided transcript assembly approach\n* Full length reads are mapped to a supplied reference genome using [minimap2](https://github.com/lh3/minimap2)\n* Transcripts are assembled by [stringtie](http://ccb.jhu.edu/software/stringtie) \nin long read mode (with or without a guide reference annotation) to generate the GFF annotation.\n* The annotation generated by the pipeline is compared to the reference annotation. \nusing [gffcompare](http://ccb.jhu.edu/software/stringtie/gffcompare.shtml)\n\n#### de novo-based transcript assembly (experimental!)\n* Sequence clusters are generated using [isONclust2](https://github.com/nanoporetech/isONclust2)\n * If a reference genome is supplied, cluster quality metrics are determined by comparing \n with clusters generated from a minimap2 alignment.\n* A consensus sequence for each cluster is generated using [spoa](https://github.com/rvaser/spoa)\n* Three rounds of polishing using racon and minimap2 to give a final polished CDS for each gene.\n* Full-length reads are then mapped to these polished CDS.\n* Transcripts are assembled by stringtie as for the reference-based approach.\n* __Note__: This approach is currently not supported with direct RNA reads.\n\n### Fusion gene detection\nFusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFFA), with the JAFFAL extension for use \nwith ONT long reads. \n\n### Differential expression analysis\n\nDifferential gene expression (DGE) and differential transcript usage (DTU) analyses aim to identify genes and/or transcripts that show statistically altered expression patterns in a studied biological system. The results of the differential analyses are presented in a quantitative format and therefore the degree of change (up or down regulation) between experimental conditions can be calculated for each gene identified.\n\nThese differential analyses work by taking a \u201csnapshot\u201d of mRNA abundance and calculating the relative levels of transcripts and isoforms. In this context, expression corresponds to the number of messenger RNAs (mRNA) measured from each gene isoform within the organism / tissue / culture being investigated. In order to determine expression levels across the whole genome, sequence data specifically targeting the mRNA molecules can be generated.\n\nOxford Nanopore Technologies provides a number of sequencing solutions to allow users to generate the required snapshot of gene expression. This can be achieved by both sequencing the mRNA directly, or via a complementary DNA (cDNA) proxy. In contrast to short read sequencing technologies, entire mRNA transcripts can be captured as single reads. The example data provided with this tutorial is from a study based on the PCR-cDNA kit. This is a robust choice for performing differential transcript usage studies. This kit is suitable for preparation of sequence libraries from low mRNA input quantities. The cDNA population is enriched through PCR with low bias; an important prerequisite for the subsequent statistical analysis.\n\n[Workflow-transcriptomes](https://github.com/epi2me-labs/wf-transcriptomes) includes a subworkflow for DGE and DTU. The first step involves using either a reference alignment or _de novo_ assembly approach to create a set of mRNA sequences per sample. These are merged into a non-redundant transcriptome using [stringtie merge](http://ccb.jhu.edu/software/stringtie). The reads are then aligned to the transcriptome using minimap2 in a splice-aware manner. [Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving per transcript counts and then the following R packages are used for analysis.\n\n### Pre-filtering of quantitative data using DRIMSeq\nDRIMSeq (Nowicka and Robinson (2016)) is used to filter the transcript count data from the salmon 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 defined in the config.yaml file. The default parameters defined for this analysis include\n* min_samps_gene_expr = 3 - a transcript must be mapped to a gene in at least this minimum number of samples for the gene be included in the analysis\n*\tmin_samps_feature_expr = 1 - a transcript must be mapped to an isoform in at least this this minimum number of samples for the gene isoform to be included in the analysis\n*\tmin_gene_expr = 10 - the minimum number of total mapped sequence reads for a gene to be considered expressed\n*\tmin_feature_expr = 3 - the minimum number of total mapped sequence reads for a gene isoform to be considered\n\n### edgeR based differential expression analysis\n+A statistical analysis is first performed using edgeR (Robinson, McCarthy, and Smyth (2010), McCarthy et al. (2012)) to identify the subset of differentially expressed genes. The filtered list of gene counts is used as input. A normalisation factor is calculated for each sequence library (using the default TMM method - please see McCarthy et al. (2012) 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))\n\n### Differential transcript usage using DEXSeq\nDifferential transcript usage analysis is performed using the R DEXSeq package (Reyes et al. (2013)). 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. \n\n### StageR stage-wise analysis of DGE and DTU\nThe final component of this isoform analysis is a stage-wise statistical test using the R software package `stageR` (Van den Berge and Clement (2018)). 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.\n\n## Running the workflow\nFor the differential expression analysis section you should have at least 3 repeats for each sample. \nYour fastq data will need to be organised in to 6 directories that represent 3 repeats for each condition. You may also need to provide a condition sheet. \n\n\n## Analysis \nDifferential gene expression is sensitive to the input data quantity and quality. There should be equivalence between samples in the number of sequence reads, mapped reads and quality scores. The sequence and alignment summary plots in the report can be used to assess these metrics. There is also a table that shows the transcript per million(TPM) calculated from the salmon counts. TPM normalizes the data for gene length and then sequencing depth, and makes it easier to compare across samples compared to counts.\n\n### Workflow inputs\n- Directory containing cDNA/direct RNA reads. Or a directory containing subdirectories each with reads from different samples\n (in fastq/fastq.gz format)\n- Reference genome in fasta format (required for reference-based assembly).\n- Optional reference annotation in GFF2/3 format (extensions allowed are .gtf(.gz), .gff(.gz), .gff3(.gz)) (required for differential expression analysis `--de_analysis`).\n- For fusion detection, JAFFAL reference files (see Quickstart) \n", + "intro": "## Introduction\n\nThis workflow identifies RNA isoforms using either cDNA or direct RNA (dRNA) \nOxford Nanopore reads.\n\n### Preprocesing\ncDNA reads are initially preprocessed by [pychopper](https://github.com/epi2me-labs/pychopper) \nfor the identification of full-length reads, as well as trimming and orientation correction (This step is omitted for \n direct RNA reads).\n\n\n### Transcript assembly\n\n#### Reference-aided transcript assembly approach\n* Full length reads are mapped to a supplied reference genome using [minimap2](https://github.com/lh3/minimap2)\n* Transcripts are assembled by [stringtie](http://ccb.jhu.edu/software/stringtie) \nin long read mode (with or without a guide reference annotation) to generate the GFF annotation.\n* The annotation generated by the pipeline is compared to the reference annotation. \nusing [gffcompare](http://ccb.jhu.edu/software/stringtie/gffcompare.shtml)\n\n#### de novo-based transcript assembly (experimental!)\n* Sequence clusters are generated using [isONclust2](https://github.com/nanoporetech/isONclust2)\n * If a reference genome is supplied, cluster quality metrics are determined by comparing \n with clusters generated from a minimap2 alignment.\n* A consensus sequence for each cluster is generated using [spoa](https://github.com/rvaser/spoa)\n* Three rounds of polishing using racon and minimap2 to give a final polished CDS for each gene.\n* Full-length reads are then mapped to these polished CDS.\n* Transcripts are assembled by stringtie as for the reference-based approach.\n* __Note__: This approach is currently not supported with direct RNA reads.\n\n### Fusion gene detection\nFusion gene detection is performed using [JAFFA](https://github.com/Oshlack/JAFFA), with the JAFFAL extension for use \nwith ONT long reads. \n\n### Differential expression analysis\n\nDifferential gene expression (DGE) and differential transcript usage (DTU) analyses aim to identify genes and/or transcripts that show statistically altered expression patterns in a studied biological system. The results of the differential analyses are presented in a quantitative format and therefore the degree of change (up or down regulation) between experimental conditions can be calculated for each gene identified.\n\nThese differential analyses work by taking a \u201csnapshot\u201d of mRNA abundance and calculating the relative levels of transcripts and isoforms. In this context, expression corresponds to the number of messenger RNAs (mRNA) measured from each gene isoform within the organism / tissue / culture being investigated. In order to determine expression levels across the whole genome, sequence data specifically targeting the mRNA molecules can be generated.\n\nOxford Nanopore Technologies provides a number of sequencing solutions to allow users to generate the required snapshot of gene expression. This can be achieved by both sequencing the mRNA directly, or via a complementary DNA (cDNA) proxy. In contrast to short read sequencing technologies, entire mRNA transcripts can be captured as single reads. The example data provided with this tutorial is from a study based on the PCR-cDNA kit. This is a robust choice for performing differential transcript usage studies. This kit is suitable for preparation of sequence libraries from low mRNA input quantities. The cDNA population is enriched through PCR with low bias; an important prerequisite for the subsequent statistical analysis.\n\n[Workflow-transcriptomes](https://github.com/epi2me-labs/wf-transcriptomes) includes a subworkflow for DGE and DTU. The first step involves using either a reference alignment or _de novo_ assembly approach to create a set of mRNA sequences per sample. These are merged into a non-redundant transcriptome using [stringtie merge](http://ccb.jhu.edu/software/stringtie). The reads are then aligned to the transcriptome using minimap2 in a splice-aware manner. [Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving per transcript counts and then the following R packages are used for analysis.\n\n### Pre-filtering of quantitative data using DRIMSeq\nDRIMSeq (Nowicka and Robinson (2016)) is used to filter the transcript count data from the salmon 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 defined in the config.yaml file. The default parameters defined for this analysis include\n* min_samps_gene_expr = 3 - a transcript must be mapped to a gene in at least this minimum number of samples for the gene be included in the analysis\n*\tmin_samps_feature_expr = 1 - a transcript must be mapped to an isoform in at least this this minimum number of samples for the gene isoform to be included in the analysis\n*\tmin_gene_expr = 10 - the minimum number of total mapped sequence reads for a gene to be considered expressed\n*\tmin_feature_expr = 3 - the minimum number of total mapped sequence reads for a gene isoform to be considered\n\n### edgeR based differential expression analysis\n+A statistical analysis is first performed using edgeR (Robinson, McCarthy, and Smyth (2010), McCarthy et al. (2012)) to identify the subset of differentially expressed genes. The filtered list of gene counts is used as input. A normalisation factor is calculated for each sequence library (using the default TMM method - please see McCarthy et al. (2012) 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))\n\n### Differential transcript usage using DEXSeq\nDifferential transcript usage analysis is performed using the R DEXSeq package (Reyes et al. (2013)). 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. \n\n### StageR stage-wise analysis of DGE and DTU\nThe final component of this isoform analysis is a stage-wise statistical test using the R software package `stageR` (Van den Berge and Clement (2018)). 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.\n\n## Running the workflow\nFor the differential expression analysis section you should have at least 3 repeats for each sample. \nYour fastq data will need to be organised in to 6 directories that represent 3 repeats for each condition. You may also need to provide a condition sheet. \n\n\n## Analysis \nDifferential gene expression is sensitive to the input data quantity and quality. There should be equivalence between samples in the number of sequence reads, mapped reads and quality scores. The sequence and alignment summary plots in the report can be used to assess these metrics. There is also a table that shows the transcript per million(TPM) calculated from the salmon counts. TPM normalizes the data for gene length and then sequencing depth, and makes it easier to compare across samples compared to counts.\n\n### Workflow inputs\n- Directory containing cDNA/direct RNA reads. Or a directory containing subdirectories each with reads from different samples\n (in fastq/fastq.gz format)\n- Reference genome in fasta format (required for reference-based assembly).\n- Optional reference annotation in GFF2/3 format (extensions allowed are .gtf(.gz), .gff(.gz), .gff3(.gz)) (required for differential expression analysis `--de_analysis`). 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.\n- For fusion detection, JAFFAL reference files (see Quickstart) \n", "links": "## Useful links\n\n* [nextflow](https://www.nextflow.io/)\n* [docker](https://www.docker.com/products/docker-desktop)\n* [Singularity](https://sylabs.io/singularity/)\n* [racon](https://github.com/isovic/racon)\n* [spoa](https://github.com/rvaser/spoa)\n* [inONclust](https://github.com/ksahlin/isONclust)\n* [isONclust2](https://github.com/nanoporetech/isONclust2)" } } \ No newline at end of file diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf index cb992c1..39c34fe 100644 --- a/subworkflows/differential_expression.nf +++ b/subworkflows/differential_expression.nf @@ -53,7 +53,7 @@ process mergeTPM { process deAnalysis { label "isoforms" errorStrategy "retry" - maxRetries 1 + maxRetries 3 input: path sample_sheet path merged_tsv @@ -68,16 +68,28 @@ process deAnalysis { script: // Just try both annotation file type because a .gff extension may be gff2(gtf) or gff3 String annotation_type = "gtf" + String strip_version = "false" if (task.attempt == 2){ annotation_type = "gff3" + strip_version = "false" + log.info("Retry deAnalysis with gff format setting.") } + else if (task.attempt == 3){ + annotation_type = "gff3" + strip_version = "true" + log.info("Retry deAnalysis with gff format setting and version removal.") + } + else if (task.attempt == 4){ + strip_version = "true" + log.info("Retry deAnalysis with gtf format setting and version removal.") + } + """ mkdir merged mkdir de_analysis mv $merged_tsv merged/all_counts.tsv mv $sample_sheet de_analysis/coldata.tsv - de_analysis.R annotation.gtf $params.min_samps_gene_expr $params.min_samps_feature_expr $params.min_gene_expr $params.min_feature_expr $annotation_type - + de_analysis.R annotation.gtf $params.min_samps_gene_expr $params.min_samps_feature_expr $params.min_gene_expr $params.min_feature_expr $annotation_type $strip_version """ }