diff --git a/CHANGELOG.md b/CHANGELOG.md index d7fb120..d881ecc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,11 +4,19 @@ 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). -## [updated] -### Changed +## [unreleased] +### Changed +- Improve differential expression outputs. +- Include transcript and gene count tables in DE_final folder. +- If differential expression subworkflow is used a non redundant transcriptome will be output which includes novel transcripts. +- Added wording to the report about how to identify novel transcripts in the DE tables. - Nextflow minimum required version to 23.04.2 - `--minimap_index_opts` parameter has been changed to `minimap2_index_opts` for consistency. +### Added +- An additional gene name column to the differential gene expression results. This is especially handy for transcriptomes where the gene ID is not the same as gene name (e.g. Ensembl). +- Wording to the report about how to identify novel transcripts in the DE tables. + ## [v0.2.1] ### Changed - Any sample aliases that contain spaces will be replaced with underscores. diff --git a/README.md b/README.md index 3b372d5..c0fb19a 100644 --- a/README.md +++ b/README.md @@ -276,6 +276,11 @@ in `${out_dir}/jaffal_output_${sample_id}` you will find: * `de_analysis/results_dtu_gene.tsv`, `de_analysis/results_dtu_transcript.tsv` and `de_analysis/results_dtu.pdf` - results of differential transcript usage by `DEXSeq`. * `de_analysis/results_dtu_stageR.tsv` - results of the `stageR` analysis of the `DEXSeq` output. * `de_analysis/dtu_plots.pdf` - DTU results plot based on the `stageR` results and filtered counts. +* `de_analysis/all_gene_counts.tsv` - Gene counts generated by the `salmon` tool before filtering. +* `de_analysis/de_transcript_counts.tsv` - Transcript counts generated by the `salmon` tool before filtering. +* `de_analysis/de_tpm_transcript_counts.tsv` - To facilitate comparisons across samples this file shows transcript per million (TPM) of the raw counts. +* `de_analysis/all_counts_filtered.tsv` - Transcript counts filtered with input criteria. Used for 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. ### References diff --git a/bin/workflow_glue/de_plots.py b/bin/workflow_glue/de_plots.py index b728686..2cabbf4 100755 --- a/bin/workflow_glue/de_plots.py +++ b/bin/workflow_glue/de_plots.py @@ -237,6 +237,13 @@ def dtu_section(dtu_file, section, gt_dic, ge_dic): section.table(dtu_results.loc[dtu_pvals.index]) +def dge_names(dge_file, geid_gname): + """Add gene name column to DGE tsv.""" + dge_results = pd.read_csv(dge_file, sep='\t') + dge_results["gene_name"] = dge_results.index.map(lambda x: geid_gname.get(x)) + dge_results.to_csv('results_dge.tsv', index=True, index_label="gene_id") + + def dge_section(dge_file, section, ids_dic): """Create DGE table and plot.""" section.markdown('### Differential gene expression') @@ -312,6 +319,7 @@ def get_translations(gtf): fn = open(gtf).readlines() gene_txid = {} gene_geid = {} + geid_gname = {} def get_feature(row, feature): return row.split(feature)[1].split( @@ -347,7 +355,8 @@ def get_translations(gtf): gene_id = gene_name gene_txid[transcript_id] = gene_name gene_geid[gene_id] = gene_reference - return gene_txid, gene_geid + geid_gname[gene_reference] = gene_name + return gene_txid, gene_geid, geid_gname def de_section( @@ -364,6 +373,10 @@ the GTF-format annotation. These counts were used to perform a statistical analysis to identify the genes and isoforms that show differences in abundance between the experimental conditions. +Any novel genes or transcripts that do not have relevant gene or transcript IDs +are prefixed with MSTRG for use in differential expression analysis. +Find the full sequences of any transcripts in the +`final_non_redundant_transcriptome.fasta` file. """) section.markdown("### Alignment summary stats") alignment_stats = pool_csvs("seqkit") @@ -371,8 +384,9 @@ the experimental conditions. alignment_summary_df = alignment_summary_df.fillna(0).applymap(np.int64) section.table(alignment_summary_df, key='alignment-stats', index=True) salmon_table(tpm, section) - gene_txid, gene_name = get_translations(stringtie) + gene_txid, gene_name, geid_gname = get_translations(stringtie) dge_section(dge, section, gene_name) + dge_names(dge, geid_gname) dexseq_section(dexseq, section, gene_name) dtu_section(dtu, section, gene_txid, gene_name) # missing dtu plots at the moment as too many diff --git a/bin/workflow_glue/report.py b/bin/workflow_glue/report.py index 88c75ed..c5a45c3 100755 --- a/bin/workflow_glue/report.py +++ b/bin/workflow_glue/report.py @@ -888,7 +888,8 @@ def de_section(report): dge = os.path.join("de_report", "results_dge.tsv") dtu = os.path.join("de_report", "results_dtu_stageR.tsv") stringtie = os.path.join("de_report", "stringtie_merged.gtf") - tpm = os.path.join("de_report", "tpm_counts.tsv") + tpm = os.path.join("de_report", "de_tpm_transcript_counts.tsv") + # This will also add a gene name column to the "results_dge.tsv" de_plots.de_section( stringtie=stringtie, dexseq=dexseq, diff --git a/docs/quickstart.md b/docs/quickstart.md index 06225b0..484e47d 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -186,6 +186,11 @@ in `${out_dir}/jaffal_output_${sample_id}` you will find: * `de_analysis/results_dtu_gene.tsv`, `de_analysis/results_dtu_transcript.tsv` and `de_analysis/results_dtu.pdf` - results of differential transcript usage by `DEXSeq`. * `de_analysis/results_dtu_stageR.tsv` - results of the `stageR` analysis of the `DEXSeq` output. * `de_analysis/dtu_plots.pdf` - DTU results plot based on the `stageR` results and filtered counts. +* `de_analysis/all_gene_counts.tsv` - Gene counts generated by the `salmon` tool before filtering. +* `de_analysis/de_transcript_counts.tsv` - Transcript counts generated by the `salmon` tool before filtering. +* `de_analysis/de_tpm_transcript_counts.tsv` - To facilitate comparisons across samples this file shows transcript per million (TPM) of the raw counts. +* `de_analysis/all_counts_filtered.tsv` - Transcript counts filtered with input criteria. Used for 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. ### References diff --git a/main.nf b/main.nf index 556e960..71b8e99 100644 --- a/main.nf +++ b/main.nf @@ -309,14 +309,14 @@ process merge_transcriptomes { path ref_annotation path ref_genome output: - path "non_redundant.fasta", emit: fasta + path "final_non_redundant_transcriptome.fasta", emit: fasta path "stringtie.gtf", emit: gtf """ stringtie --merge -G $ref_annotation -p ${task.cpus} -o stringtie.gtf query_annotations/* seqkit subseq --feature "transcript" --gtf-tag "transcript_id" --gtf stringtie.gtf $ref_genome > temp_transcriptome.fasta seqkit rmdup -s < temp_transcriptome.fasta > temp_del_repeats.fasta cat temp_del_repeats.fasta | sed 's/>.* />/' | sed -e 's/_[0-9]* \\[/ \\[/' > temp_rm_empty_seq.fasta - awk 'BEGIN {RS = ">" ; FS = "\\n" ; ORS = ""} \$2 {print ">"\$0}' temp_rm_empty_seq.fasta > non_redundant.fasta + awk 'BEGIN {RS = ">" ; FS = "\\n" ; ORS = ""} \$2 {print ">"\$0}' temp_rm_empty_seq.fasta > "final_non_redundant_transcriptome.fasta" rm temp_transcriptome.fasta rm temp_del_repeats.fasta rm temp_rm_empty_seq.fasta @@ -342,6 +342,9 @@ process makeReport { path "seqkit/*" output: path("wf-transcriptomes-*.html"), emit: report + // for DE analysis, a `gene_name` column will be added to + // `de_report/results_dge.tsv` + path "results_dge.tsv", emit: de_analysis, optional: true script: // Convert the sample_id arrayList. sids = new BlankSeparatedList(sample_ids) @@ -387,7 +390,6 @@ process makeReport { \$OPT_JAFFAL_CSV \ $OPT_DENOVO \ \$dereport - """ } @@ -495,7 +497,7 @@ workflow pipeline { return l } - + software_versions = getVersions() workflow_params = getParams() @@ -578,7 +580,7 @@ workflow pipeline { gtf = merge_transcriptomes.out.gtf } else { - transcriptome = ref_transcriptome + transcriptome = Channel.fromPath(ref_transcriptome) gtf = ref_annotation } de = differential_expression(transcriptome, input_reads, sample_sheet, gtf) @@ -586,6 +588,7 @@ workflow pipeline { count_transcripts_file = de.count_transcripts dtu_plots = de.dtu_plots de_outputs = de.de_outputs + counts = de.counts } else{ de_report = file("$projectDir/data/OPTIONAL_FILE") count_transcripts_file = file("$projectDir/data/OPTIONAL_FILE") @@ -605,6 +608,8 @@ workflow pipeline { count_transcripts_file) report = makeReport.out.report + + results = results.concat(makeReport.out.report) @@ -645,15 +650,19 @@ workflow pipeline { .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){ - results = results.concat(de.dtu_plots, de_outputs) + de_update = makeReport.out.de_analysis + de_results = report.concat(transcriptome, de_outputs.flatten(), counts.flatten(), de_update) + results = results.concat(de_results.map{ [it, "de_analysis"] }) } - results = fastq_ingress_results.map { [it, "fastq_ingress_results"] }.concat(results.map{ [it, null]}) + results.concat(workflow_params.map{ [it, null]}) + emit: results - telemetry = workflow_params } // entrypoint workflow diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf index c96d784..c6d0664 100644 --- a/subworkflows/differential_expression.nf +++ b/subworkflows/differential_expression.nf @@ -32,9 +32,9 @@ process mergeCounts { input: path counts output: - path "all_counts.tsv" + path "de_transcript_counts.tsv" """ - workflow-glue merge_count_tsvs -z -o all_counts.tsv -tsvs ${counts} + workflow-glue merge_count_tsvs -z -o de_transcript_counts.tsv -tsvs ${counts} """ } @@ -43,9 +43,9 @@ process mergeTPM { input: path counts output: - path "tpm_counts.tsv" + path "de_tpm_transcript_counts.tsv" """ - workflow-glue merge_count_tsvs -o tpm_counts.tsv -z -tpm True -tsvs $counts + workflow-glue merge_count_tsvs -o de_tpm_transcript_counts.tsv -z -tpm True -tsvs $counts """ } @@ -104,7 +104,7 @@ process plotResults { output: path "de_analysis/dtu_plots.pdf", emit: dtu_plots path "sample_sheet.tsv", emit: sample_sheet_csv - path "de_analysis", emit: stageR + path "de_analysis/*", emit: stageR """ mkdir merged mv $sample_sheet de_analysis/coldata.tsv @@ -173,9 +173,11 @@ workflow differential_expression { analysis.stageR).combine(plotResults.out.sample_sheet_csv).combine(merged).combine( ref_annotation).combine(merged_TPM) count_transcripts_file = count_transcripts.out.seqkit_stats.collect() + all_counts = merged_TPM.concat(merged, analysis.flt_counts, analysis.gene_counts) emit: all_de = de_report count_transcripts = count_transcripts_file dtu_plots = plotResults.out.dtu_plots de_outputs = plotResults.out.stageR + counts = all_counts } diff --git a/test_data/sample_sheet.csv b/test_data/sample_sheet.csv index 00adea1..61f376c 100644 --- a/test_data/sample_sheet.csv +++ b/test_data/sample_sheet.csv @@ -4,4 +4,4 @@ barcode02,sample02,sample02,untreated barcode03,sample03,sample03,untreated barcode04,sample04,sample04,treated barcode05,sample05,sample05,treated -barcode06,sample06,sample06,treated \ No newline at end of file +barcode06,sample06,sample06,treated