Merge branch 'CW-3075' into 'dev'
CW-3075 improve counts output files Closes CW-3075 See merge request epi2melabs/workflows/wf-transcriptomes!149
This commit is contained in:
Sarah Griffiths
commit
077ad86a92
4
.github/ISSUE_TEMPLATE/bug_report.yml
vendored
4
.github/ISSUE_TEMPLATE/bug_report.yml
vendored
@ -128,8 +128,8 @@ body:
|
||||
label: Were you able to successfully run the latest version of the workflow with the demo data?
|
||||
description: For CLI execution, were you able to successfully run the workflow using the demo data available in the [Install and run](./README.md#install-and-run) section of the `README.md`? For execution in the EPI2ME application, were you able to successfully run the workflow via the "Use demo data" button?
|
||||
options:
|
||||
- yes
|
||||
- no
|
||||
- 'yes'
|
||||
- 'no'
|
||||
- other (please describe below)
|
||||
validations:
|
||||
required: true
|
||||
|
||||
11
CHANGELOG.md
11
CHANGELOG.md
@ -4,7 +4,16 @@ 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]
|
||||
## [unreleased]
|
||||
### Changed
|
||||
- Improve documentation around filtering of transcripts done before DTU analysis.
|
||||
- Renamed files:
|
||||
- `de_analysis/all_counts_filtered.tsv` to `de_analysis/filtered_transcript_counts_with_genes.tsv`
|
||||
- `de_analysis/de_tpm_transcript_counts.tsv` to `de_analysis/unfiltered_tpm_transcript_counts.tsv`
|
||||
### Added
|
||||
- Output additional `de_analysis/cpm_gene_counts.tsv` with counts per million gene counts.
|
||||
- Output additional `de_analysis/unfiltered_transcript_counts_with_genes.tsv` with unfiltered transcript counts with associated gene IDs.
|
||||
- Add gene name column to the de_analysis counts TSV files.
|
||||
### Fixed
|
||||
- Mapping stage using a single thread only.
|
||||
|
||||
|
||||
25
README.md
25
README.md
@ -156,10 +156,10 @@ input_reads.fastq ─── input_directory ─── input_directory
|
||||
| Nextflow parameter name | Type | Description | Help | Default |
|
||||
|--------------------------|------|-------------|------|---------|
|
||||
| de_analysis | boolean | Run DE anaylsis | Running this requires you to provide at least two replicates for a control and treated sample as well as a sample sheet param. | False |
|
||||
| min_gene_expr | integer | Minimum gene counts | The minimum number of total mapped sequence reads for a gene to be considered expressed. | 10 |
|
||||
| min_feature_expr | integer | Minimum transcript counts | The minimum number of total mapped sequence reads for a transcript to be considered. | 3 |
|
||||
| min_samps_gene_expr | integer | Genes expressed in a minimum of this many samples will be included in the differential expression analysis. | A gene must be mapped to at least this minimum number of samples for the gene be included in the analysis. | 3 |
|
||||
| min_samps_feature_expr | integer | Transcripts expressed in minimum this many samples | A transcript must be mapped in at least this this minimum number of samples to be included in the analysis. | 1 |
|
||||
| min_gene_expr | integer | The minimum number of total mapped sequence reads required for a gene to be considered in differential transcript usage analysis. | Filtering at the gene level ensures that the observed transcript ratios are calculated with a minimum number of counts per gene. | 10 |
|
||||
| min_feature_expr | integer | The minimum number of reads assigned to a transcript for it to be considered in differential transcript usage analysis. | Filter out transcripts that do not have this minimum number of transcript expression, reducing noise. | 3 |
|
||||
| min_samps_gene_expr | integer | Set the minimum number of samples in which a gene is expressed to be included in the differential transcript usage analysis. | A gene must be expressed in at least this number of samples for the gene be included in the differential transcript usage analysis. Filtering at the gene level improves the reliability of the observed transcript ratios. | 3 |
|
||||
| min_samps_feature_expr | integer | Set the minimum number of samples in which a transcript is expressed to be included in the differential transcript usage analysis. | A transcript must expressed in at least this minimum number of samples to be included in the analysis. Should be equal to the number of replicates per sample you have. | 1 |
|
||||
|
||||
|
||||
### Advanced Options
|
||||
@ -201,7 +201,7 @@ Output files may be aggregated including information for all samples or provided
|
||||
| Annotated assembled transcriptome | {{ alias }}_merged_transcriptome.fas | Per sample annotated assembled transcriptome. | per-sample |
|
||||
| Alignment summary statistics | {{ alias }}_read_aln_stats.tsv | Per sample alignment summary statistics. | per-sample |
|
||||
| GFF compare results. | {{ alias }}_gffcompare | All GFF compare output files. | per-sample |
|
||||
| Differential gene expression results | /de_analysis/results_dge.tsv | This is a gene-level result file that describes genes and the probability that they show differential expression between experimental conditions . | aggregated |
|
||||
| Differential gene expression results | /de_analysis/results_dge.tsv | This is a gene-level result file that describes genes and their probability of showing differential expression between experimental conditions. | aggregated |
|
||||
| Differential gene expression report | /de_analysis/results_dge.pdf | Summary report of differential gene expression analysis as a PDF. | aggregated |
|
||||
| Differential transcript usage gene TSV | /de_analysis/results_dtu_gene.tsv | This is a gene-level result file from DEXSeq that lists annotated genes and their probabilities of differential expression. | aggregated |
|
||||
| Differential transcript usage report | /de_analysis/results_dtu.pdf | Summary report of differential transcript usage results as a PDF. | aggregated |
|
||||
@ -209,9 +209,10 @@ Output files may be aggregated including information for all samples or provided
|
||||
| Differential transcript usage stageR TSV | /de_analysis/results_dtu_stageR.tsv | This is the output from StageR and it shows both gene and transcript probabilities of differential expression | aggregated |
|
||||
| Differential transcript usage DEXSeq TSV | /de_analysis/results_dexseq.tsv | The complete output from the DEXSeq-analysis, shows both gene and transcript probabilities of differential expression. | aggregated |
|
||||
| Gene counts | /de_analysis/all_gene_counts.tsv | Raw gene counts created by the Salmon tool, before filtering. | aggregated |
|
||||
| Transcript counts | /de_analysis/all_transcript_counts.tsv | Raw transcript counts created by the Salmon tool, before filtering. | aggregated |
|
||||
| Transcript counts filtered | /de_analysis/all_counts_filtered.tsv | Filtered transcript counts, used for DE_analysis. | aggregated |
|
||||
| Transcript per million counts | /de_analysis/de_tpm_transcript_counts.tsv | This file shows transcript per million (TPM) of the raw counts to facilitate comparisons across sample. | aggregated |
|
||||
| Gene counts per million | /de_analysis/cpm_gene_counts.tsv | This file shows counts per million (CPM) of the raw gene counts to facilitate comparisons across samples. | aggregated |
|
||||
| Transcript counts | /de_analysis/unfiltered_transcript_counts_with_genes.tsv | Raw transcript counts created by the Salmon tool, before filtering. Includes reference to the associated gene ID. | aggregated |
|
||||
| Transcript per million counts | /de_analysis/unfiltered_tpm_transcript_counts.tsv | This file shows transcripts per million (TPM) of the raw counts to facilitate comparisons across samples. | aggregated |
|
||||
| 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 |
|
||||
| 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 |
|
||||
@ -283,11 +284,11 @@ The reads from all the samples will be aligned with the final non redundant tran
|
||||
#### 5.3 Count genes and transcripts
|
||||
[Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving gene and transcript counts.
|
||||
|
||||
#### 5.4 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. 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`.
|
||||
#### 5.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 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. The filtered list of gene counts is used as input. A normalisation factor is calculated for each sequence library (using the default TMM method described by [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
|
||||
[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
|
||||
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.
|
||||
|
||||
@ -2,6 +2,7 @@
|
||||
|
||||
suppressMessages(library("DRIMSeq"))
|
||||
suppressMessages(library("GenomicFeatures"))
|
||||
suppressMessages(library("edgeR"))
|
||||
args <- commandArgs(trailingOnly=TRUE)
|
||||
ref_annotation <- args[1]
|
||||
min_samps_gene_expr <- args[2]
|
||||
@ -12,11 +13,11 @@ 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))
|
||||
cts <- as.matrix(read.csv("all_counts.tsv", sep="\t", row.names="Reference", stringsAsFactors=FALSE))
|
||||
|
||||
# Set up sample data frame:
|
||||
#changed this to sample_id
|
||||
coldata <- read.csv("de_analysis/coldata.tsv", row.names="alias", sep=",", stringsAsFactors=TRUE)
|
||||
coldata <- read.csv("sample_sheet.csv", row.names="alias", sep=",", stringsAsFactors=TRUE)
|
||||
|
||||
coldata$sample_id <- rownames(coldata)
|
||||
# check if control condition exists, sets as reference
|
||||
@ -46,6 +47,9 @@ rownames(txdf) <- NULL
|
||||
# Create counts data frame:
|
||||
counts<-data.frame(gene_id=txdf$GENEID, feature_id=txdf$TXNAME, cts)
|
||||
|
||||
# output unfiltered version of the counts table now we have paired transcripts with gene ids
|
||||
write.table(counts, file="de_analysis/unfiltered_transcript_counts_with_genes.tsv", sep="\t", row.names = FALSE, quote=FALSE)
|
||||
|
||||
cat("Filtering counts using DRIMSeq.\n")
|
||||
|
||||
d <- dmDSdata(counts=counts, samples=coldata)
|
||||
@ -64,15 +68,22 @@ suppressMessages(library("dplyr"))
|
||||
# Sum transcript counts into gene counts:
|
||||
cat("Sum transcript counts into gene counts.\n")
|
||||
trs_cts <- counts(d)
|
||||
write.table(trs_cts, file="merged/all_counts_filtered.tsv",sep="\t")
|
||||
write.table(trs_cts, file="merged/filtered_transcript_counts_with_genes.tsv", sep="\t", row.names = FALSE, quote=FALSE)
|
||||
|
||||
gene_cts <- trs_cts_unfiltered %>% dplyr::select(c(1, 3:ncol(trs_cts))) %>% group_by(gene_id) %>% summarise_all(tibble::lst(sum)) %>% data.frame()
|
||||
rownames(gene_cts) <- gene_cts$gene_id
|
||||
gene_cts$gene_id <- NULL
|
||||
write.table(gene_cts, file="merged/all_gene_counts.tsv",sep="\t")
|
||||
write.table(gene_cts, file="merged/all_gene_counts.tsv", sep="\t", quote=FALSE)
|
||||
|
||||
# Output count per million of the gene counts using edgeR CPM
|
||||
cpm_gene_counts <- cpm(gene_cts)
|
||||
# Add gene_id as index column header
|
||||
cpm_gene_counts <- cbind(var_name = rownames(cpm_gene_counts), cpm_gene_counts)
|
||||
rownames(cpm_gene_counts) <- NULL
|
||||
colnames(cpm_gene_counts)[1] <- "gene_id"
|
||||
write.table(cpm_gene_counts, file="de_analysis/cpm_gene_counts.tsv", sep="\t", quote=FALSE, row.names = FALSE)
|
||||
|
||||
# Differential gene expression using edgeR:
|
||||
suppressMessages(library("edgeR"))
|
||||
cat("Running differential gene expression analysis using edgeR.\n")
|
||||
|
||||
y <- DGEList(gene_cts)
|
||||
|
||||
@ -5,18 +5,17 @@ suppressMessages(library(ggplot2))
|
||||
suppressMessages(library(tidyr))
|
||||
|
||||
# Set up sample data frame:
|
||||
coldata <- read.csv("de_analysis/coldata.tsv", row.names="alias", sep=",")
|
||||
coldata$sample_id <- rownames(coldata)
|
||||
coldata <- read.csv("sample_sheet.tsv", row.names="alias", sep=",")
|
||||
coldata$condition <- factor(coldata$condition, levels=rev(levels(coldata$condition)))
|
||||
coldata$type <-NULL
|
||||
coldata$patient <-NULL
|
||||
|
||||
# Read stageR results:
|
||||
stageR <- read.csv("de_analysis/results_dtu_stageR.tsv", sep="\t")
|
||||
stageR <- read.csv("results_dtu_stageR.tsv", sep="\t")
|
||||
names(stageR) <- c("gene_id", "transcript_id", "p_gene", "p_transcript");
|
||||
|
||||
# Read filtered counts:
|
||||
counts <- read.csv("merged/all_counts_filtered.tsv", sep="\t");
|
||||
counts <- read.csv("filtered_transcript_counts_with_genes.tsv", sep="\t");
|
||||
names(counts)[2]<-"transcript_id"
|
||||
|
||||
# Join counts and stageR results:
|
||||
|
||||
@ -237,11 +237,29 @@ 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."""
|
||||
def add_gene_names_dge(dge_file, geid_gname):
|
||||
"""Add gene name column to DGE TSV with gene ID to gene Name dict."""
|
||||
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")
|
||||
column_to_move = dge_results.index.map(
|
||||
lambda x: geid_gname.get(x))
|
||||
dge_results.insert(0, "gene_name", column_to_move)
|
||||
dge_results.to_csv(dge_file, index=True, index_label="gene_id", sep="\t")
|
||||
|
||||
|
||||
def add_gene_names_filtered(unfiltered_file, geid_gname):
|
||||
"""Use gene id to gene name dict to add name column to counts TSV."""
|
||||
unfiltered = pd.read_csv(unfiltered_file, sep='\t')
|
||||
unfiltered.insert(1, "gene_name", unfiltered.gene_id.map(
|
||||
lambda x: geid_gname.get(x)))
|
||||
unfiltered.to_csv(unfiltered_file, index=False, sep='\t')
|
||||
|
||||
|
||||
def add_tpm_names(tpm_file, txid_gname):
|
||||
"""Use transcript id to gene name dict to add gname column to TPM TSV."""
|
||||
tpm = pd.read_csv(tpm_file, sep='\t')
|
||||
tpm.insert(1, "gene_name", tpm.Reference.map(
|
||||
lambda x: txid_gname.get(x)))
|
||||
tpm.to_csv(tpm_file, index=False, sep='\t')
|
||||
|
||||
|
||||
def dge_section(dge_file, section, ids_dic):
|
||||
@ -361,7 +379,8 @@ def get_translations(gtf):
|
||||
|
||||
def de_section(
|
||||
stringtie, dge, dexseq, dtu,
|
||||
tpm, report):
|
||||
tpm, report, filtered, unfiltered,
|
||||
gene_counts):
|
||||
"""Differential expression sections."""
|
||||
section = report.add_section()
|
||||
section.markdown("# Differential expression.")
|
||||
@ -385,8 +404,16 @@ Find the full sequences of any transcripts in the
|
||||
section.table(alignment_summary_df, key='alignment-stats', index=True)
|
||||
salmon_table(tpm, section)
|
||||
gene_txid, gene_name, geid_gname = get_translations(stringtie)
|
||||
|
||||
# Use dictionaries to add gene names to the counts tsv files to help users
|
||||
add_gene_names_dge(dge, geid_gname)
|
||||
add_gene_names_dge(gene_counts, geid_gname)
|
||||
add_gene_names_filtered(filtered, geid_gname)
|
||||
add_gene_names_filtered(unfiltered, geid_gname)
|
||||
add_tpm_names(tpm, gene_txid)
|
||||
|
||||
# Add tables to report
|
||||
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
|
||||
|
||||
@ -882,15 +882,23 @@ 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", "de_tpm_transcript_counts.tsv")
|
||||
# This will also add a gene name column to the "results_dge.tsv"
|
||||
tpm = os.path.join("de_report", "unfiltered_tpm_transcript_counts.tsv")
|
||||
filtered = os.path.join(
|
||||
"de_report", "filtered_transcript_counts_with_genes.tsv")
|
||||
unfiltered = os.path.join(
|
||||
"de_report", "unfiltered_transcript_counts_with_genes.tsv")
|
||||
gene_counts = os.path.join("de_report", "all_gene_counts.tsv")
|
||||
# This will also add a gene name column to the above counts tsv files
|
||||
de_plots.de_section(
|
||||
stringtie=stringtie,
|
||||
dexseq=dexseq,
|
||||
dge=dge,
|
||||
dtu=dtu,
|
||||
tpm=tpm,
|
||||
report=report)
|
||||
report=report,
|
||||
filtered=filtered,
|
||||
unfiltered=unfiltered,
|
||||
gene_counts=gene_counts)
|
||||
|
||||
|
||||
def main(args):
|
||||
|
||||
@ -45,10 +45,10 @@
|
||||
| Nextflow parameter name | Type | Description | Help | Default |
|
||||
|--------------------------|------|-------------|------|---------|
|
||||
| de_analysis | boolean | Run DE anaylsis | Running this requires you to provide at least two replicates for a control and treated sample as well as a sample sheet param. | False |
|
||||
| min_gene_expr | integer | Minimum gene counts | The minimum number of total mapped sequence reads for a gene to be considered expressed. | 10 |
|
||||
| min_feature_expr | integer | Minimum transcript counts | The minimum number of total mapped sequence reads for a transcript to be considered. | 3 |
|
||||
| min_samps_gene_expr | integer | Genes expressed in a minimum of this many samples will be included in the differential expression analysis. | A gene must be mapped to at least this minimum number of samples for the gene be included in the analysis. | 3 |
|
||||
| min_samps_feature_expr | integer | Transcripts expressed in minimum this many samples | A transcript must be mapped in at least this this minimum number of samples to be included in the analysis. | 1 |
|
||||
| min_gene_expr | integer | The minimum number of total mapped sequence reads required for a gene to be considered in differential transcript usage analysis. | Filtering at the gene level ensures that the observed transcript ratios are calculated with a minimum number of counts per gene. | 10 |
|
||||
| min_feature_expr | integer | The minimum number of reads assigned to a transcript for it to be considered in differential transcript usage analysis. | Filter out transcripts that do not have this minimum number of transcript expression, reducing noise. | 3 |
|
||||
| min_samps_gene_expr | integer | Set the minimum number of samples in which a gene is expressed to be included in the differential transcript usage analysis. | A gene must be expressed in at least this number of samples for the gene be included in the differential transcript usage analysis. Filtering at the gene level improves the reliability of the observed transcript ratios. | 3 |
|
||||
| min_samps_feature_expr | integer | Set the minimum number of samples in which a transcript is expressed to be included in the differential transcript usage analysis. | A transcript must expressed in at least this minimum number of samples to be included in the analysis. Should be equal to the number of replicates per sample you have. | 1 |
|
||||
|
||||
|
||||
### Advanced Options
|
||||
|
||||
@ -12,7 +12,7 @@ Output files may be aggregated including information for all samples or provided
|
||||
| Annotated assembled transcriptome | {{ alias }}_merged_transcriptome.fas | Per sample annotated assembled transcriptome. | per-sample |
|
||||
| Alignment summary statistics | {{ alias }}_read_aln_stats.tsv | Per sample alignment summary statistics. | per-sample |
|
||||
| GFF compare results. | {{ alias }}_gffcompare | All GFF compare output files. | per-sample |
|
||||
| Differential gene expression results | /de_analysis/results_dge.tsv | This is a gene-level result file that describes genes and the probability that they show differential expression between experimental conditions . | aggregated |
|
||||
| Differential gene expression results | /de_analysis/results_dge.tsv | This is a gene-level result file that describes genes and their probability of showing differential expression between experimental conditions. | aggregated |
|
||||
| Differential gene expression report | /de_analysis/results_dge.pdf | Summary report of differential gene expression analysis as a PDF. | aggregated |
|
||||
| Differential transcript usage gene TSV | /de_analysis/results_dtu_gene.tsv | This is a gene-level result file from DEXSeq that lists annotated genes and their probabilities of differential expression. | aggregated |
|
||||
| Differential transcript usage report | /de_analysis/results_dtu.pdf | Summary report of differential transcript usage results as a PDF. | aggregated |
|
||||
@ -20,9 +20,10 @@ Output files may be aggregated including information for all samples or provided
|
||||
| Differential transcript usage stageR TSV | /de_analysis/results_dtu_stageR.tsv | This is the output from StageR and it shows both gene and transcript probabilities of differential expression | aggregated |
|
||||
| Differential transcript usage DEXSeq TSV | /de_analysis/results_dexseq.tsv | The complete output from the DEXSeq-analysis, shows both gene and transcript probabilities of differential expression. | aggregated |
|
||||
| Gene counts | /de_analysis/all_gene_counts.tsv | Raw gene counts created by the Salmon tool, before filtering. | aggregated |
|
||||
| Transcript counts | /de_analysis/all_transcript_counts.tsv | Raw transcript counts created by the Salmon tool, before filtering. | aggregated |
|
||||
| Transcript counts filtered | /de_analysis/all_counts_filtered.tsv | Filtered transcript counts, used for DE_analysis. | aggregated |
|
||||
| Transcript per million counts | /de_analysis/de_tpm_transcript_counts.tsv | This file shows transcript per million (TPM) of the raw counts to facilitate comparisons across sample. | aggregated |
|
||||
| Gene counts per million | /de_analysis/cpm_gene_counts.tsv | This file shows counts per million (CPM) of the raw gene counts to facilitate comparisons across samples. | aggregated |
|
||||
| Transcript counts | /de_analysis/unfiltered_transcript_counts_with_genes.tsv | Raw transcript counts created by the Salmon tool, before filtering. Includes reference to the associated gene ID. | aggregated |
|
||||
| Transcript per million counts | /de_analysis/unfiltered_tpm_transcript_counts.tsv | This file shows transcripts per million (TPM) of the raw counts to facilitate comparisons across samples. | aggregated |
|
||||
| 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 |
|
||||
| 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 |
|
||||
|
||||
@ -60,11 +60,11 @@ The reads from all the samples will be aligned with the final non redundant tran
|
||||
#### 5.3 Count genes and transcripts
|
||||
[Salmon](https://github.com/COMBINE-lab/salmon) is used for transcript quantification, giving gene and transcript counts.
|
||||
|
||||
#### 5.4 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. 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`.
|
||||
#### 5.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 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. The filtered list of gene counts is used as input. A normalisation factor is calculated for each sequence library (using the default TMM method described by [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
|
||||
[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
|
||||
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.
|
||||
|
||||
2
main.nf
2
main.nf
@ -378,7 +378,7 @@ process makeReport {
|
||||
|
||||
label "isoforms"
|
||||
cpus 2
|
||||
memory "2 GB"
|
||||
memory "4 GB"
|
||||
|
||||
input:
|
||||
path versions
|
||||
|
||||
@ -227,29 +227,29 @@
|
||||
"type": "integer",
|
||||
"title": "Minimum gene expression",
|
||||
"default": 10,
|
||||
"description": "Minimum gene counts",
|
||||
"help_text": "The minimum number of total mapped sequence reads for a gene to be considered expressed."
|
||||
"description": "The minimum number of total mapped sequence reads required for a gene to be considered in differential transcript usage analysis.",
|
||||
"help_text": "Filtering at the gene level ensures that the observed transcript ratios are calculated with a minimum number of counts per gene."
|
||||
},
|
||||
"min_feature_expr": {
|
||||
"type": "integer",
|
||||
"title": "Minimum feature expression",
|
||||
"default": 3,
|
||||
"description": "Minimum transcript counts",
|
||||
"help_text": "The minimum number of total mapped sequence reads for a transcript to be considered."
|
||||
"description": "The minimum number of reads assigned to a transcript for it to be considered in differential transcript usage analysis.",
|
||||
"help_text": "Filter out transcripts that do not have this minimum number of transcript expression, reducing noise."
|
||||
},
|
||||
"min_samps_gene_expr": {
|
||||
"type": "integer",
|
||||
"title": "Minimum samples with gene expression",
|
||||
"description": "Genes expressed in a minimum of this many samples will be included in the differential expression analysis.",
|
||||
"description": "Set the minimum number of samples in which a gene is expressed to be included in the differential transcript usage analysis.",
|
||||
"default": 3,
|
||||
"help_text": "A gene must be mapped to at least this minimum number of samples for the gene be included in the analysis."
|
||||
"help_text": "A gene must be expressed in at least this number of samples for the gene be included in the differential transcript usage analysis. Filtering at the gene level improves the reliability of the observed transcript ratios."
|
||||
},
|
||||
"min_samps_feature_expr": {
|
||||
"type": "integer",
|
||||
"title": "Minimum samples with feature expression",
|
||||
"default": 1,
|
||||
"description": "Transcripts expressed in minimum this many samples",
|
||||
"help_text": "A transcript must be mapped in at least this this minimum number of samples to be included in the analysis."
|
||||
"description": "Set the minimum number of samples in which a transcript is expressed to be included in the differential transcript usage analysis.",
|
||||
"help_text": "A transcript must expressed in at least this minimum number of samples to be included in the analysis. Should be equal to the number of replicates per sample you have."
|
||||
}
|
||||
}
|
||||
},
|
||||
|
||||
@ -83,7 +83,7 @@
|
||||
"dge-results-tsv": {
|
||||
"filepath": "/de_analysis/results_dge.tsv",
|
||||
"title": "Differential gene expression results",
|
||||
"description": "This is a gene-level result file that describes genes and the probability that they show differential expression between experimental conditions .",
|
||||
"description": "This is a gene-level result file that describes genes and their probability of showing differential expression between experimental conditions.",
|
||||
"mime-type": "text/tab-separated-values",
|
||||
"optional": true,
|
||||
"type": "aggregated"
|
||||
@ -144,26 +144,34 @@
|
||||
"optional": true,
|
||||
"type": "aggregated"
|
||||
},
|
||||
"transcipt_counts": {
|
||||
"filepath": "/de_analysis/all_transcript_counts.tsv",
|
||||
"title": "Transcript counts",
|
||||
"description": "Raw transcript counts created by the Salmon tool, before filtering.",
|
||||
"gene_counts_per_million": {
|
||||
"filepath": "/de_analysis/cpm_gene_counts.tsv",
|
||||
"title": "Gene counts per million",
|
||||
"description": "This file shows counts per million (CPM) of the raw gene counts to facilitate comparisons across samples.",
|
||||
"mime-type": "text/tab-separated-values",
|
||||
"optional": true,
|
||||
"type": "aggregated"
|
||||
},
|
||||
"transcipt_counts_filtered": {
|
||||
"filepath": "/de_analysis/all_counts_filtered.tsv",
|
||||
"title": "Transcript counts filtered",
|
||||
"description": "Filtered transcript counts, used for DE_analysis.",
|
||||
"transcript_counts": {
|
||||
"filepath": "/de_analysis/unfiltered_transcript_counts_with_genes.tsv",
|
||||
"title": "Transcript counts",
|
||||
"description": "Raw transcript counts created by the Salmon tool, before filtering. Includes reference to the associated gene ID.",
|
||||
"mime-type": "text/tab-separated-values",
|
||||
"optional": true,
|
||||
"type": "aggregated"
|
||||
},
|
||||
"tpm_transcript_counts": {
|
||||
"filepath": "/de_analysis/de_tpm_transcript_counts.tsv",
|
||||
"filepath": "/de_analysis/unfiltered_tpm_transcript_counts.tsv",
|
||||
"title": "Transcript per million counts",
|
||||
"description": "This file shows transcript per million (TPM) of the raw counts to facilitate comparisons across sample.",
|
||||
"description": "This file shows transcripts per million (TPM) of the raw counts to facilitate comparisons across samples.",
|
||||
"mime-type": "text/tab-separated-values",
|
||||
"optional": true,
|
||||
"type": "aggregated"
|
||||
},
|
||||
"transcipt_counts_filtered": {
|
||||
"filepath": "/de_analysis/filtered_transcript_counts_with_genes.tsv",
|
||||
"title": "Transcript counts filtered",
|
||||
"description": "Filtered transcript counts, used for differential transcript usage analysis. Includes a reference to the associated gene ID.",
|
||||
"mime-type": "text/tab-separated-values",
|
||||
"optional": true,
|
||||
"type": "aggregated"
|
||||
|
||||
@ -37,9 +37,9 @@ process mergeCounts {
|
||||
input:
|
||||
path counts
|
||||
output:
|
||||
path "de_transcript_counts.tsv"
|
||||
path "unfiltered_transcript_counts.tsv"
|
||||
"""
|
||||
workflow-glue merge_count_tsvs -z -o de_transcript_counts.tsv -tsvs ${counts}
|
||||
workflow-glue merge_count_tsvs -z -o unfiltered_transcript_counts.tsv -tsvs ${counts}
|
||||
"""
|
||||
}
|
||||
|
||||
@ -50,9 +50,10 @@ process mergeTPM {
|
||||
input:
|
||||
path counts
|
||||
output:
|
||||
path "de_tpm_transcript_counts.tsv"
|
||||
path "unfiltered_tpm_transcript_counts.tsv"
|
||||
// Use tpm parameter with merge_counts_tsvs.py to out transcript per million file
|
||||
"""
|
||||
workflow-glue merge_count_tsvs -o de_tpm_transcript_counts.tsv -z -tpm True -tsvs $counts
|
||||
workflow-glue merge_count_tsvs -o unfiltered_tpm_transcript_counts.tsv -z -tpm True -tsvs $counts
|
||||
"""
|
||||
}
|
||||
|
||||
@ -64,12 +65,13 @@ process deAnalysis {
|
||||
cpus 4
|
||||
memory "16 GB"
|
||||
input:
|
||||
path sample_sheet
|
||||
path merged_tsv
|
||||
path "sample_sheet.csv"
|
||||
path "all_counts.tsv"
|
||||
path "annotation.gtf"
|
||||
output:
|
||||
path "de_analysis/results_dtu_stageR.tsv", emit: stageR
|
||||
path "merged/all_counts_filtered.tsv", emit: flt_counts
|
||||
path "merged/filtered_transcript_counts_with_genes.tsv", emit: flt_counts
|
||||
path "de_analysis/unfiltered_transcript_counts_with_genes.tsv", emit: unflt_counts
|
||||
path "merged/all_gene_counts.tsv", emit: gene_counts
|
||||
path "de_analysis/results_dge.tsv", emit: dge
|
||||
path "de_analysis/results_dexseq.tsv", emit: dexseq
|
||||
@ -96,8 +98,6 @@ process deAnalysis {
|
||||
"""
|
||||
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 $strip_version
|
||||
"""
|
||||
}
|
||||
@ -108,20 +108,16 @@ process plotResults {
|
||||
cpus 2
|
||||
memory "2 GB"
|
||||
input:
|
||||
path flt_count
|
||||
path res_dtu
|
||||
path sample_sheet
|
||||
path "filtered_transcript_counts_with_genes.tsv"
|
||||
path "results_dtu_stageR.tsv"
|
||||
path "sample_sheet.tsv"
|
||||
path de_analysis
|
||||
output:
|
||||
path "de_analysis/dtu_plots.pdf", emit: dtu_plots
|
||||
path "sample_sheet.tsv", emit: sample_sheet_csv
|
||||
path "de_analysis/*", emit: stageR
|
||||
"""
|
||||
mkdir merged
|
||||
mv $sample_sheet de_analysis/coldata.tsv
|
||||
mv $flt_count merged/all_counts_filtered.tsv
|
||||
plot_dtu_results.R
|
||||
mv de_analysis/coldata.tsv sample_sheet.tsv
|
||||
"""
|
||||
}
|
||||
|
||||
@ -182,11 +178,10 @@ workflow differential_expression {
|
||||
merged_TPM = mergeTPM(count_transcripts.out.counts.collect())
|
||||
analysis = deAnalysis(sample_sheet, merged, ref_annotation)
|
||||
plotResults(analysis.flt_counts, analysis.stageR, sample_sheet, analysis.de_analysis)
|
||||
de_report = analysis.flt_counts.combine(analysis.gene_counts).combine(analysis.dge).combine(analysis.dexseq).combine(
|
||||
analysis.stageR).combine(plotResults.out.sample_sheet_csv).combine(merged).combine(
|
||||
ref_annotation).combine(merged_TPM)
|
||||
de_report = analysis.flt_counts.concat(analysis.gene_counts, analysis.dge, analysis.dexseq,
|
||||
analysis.stageR, plotResults.out.sample_sheet_csv, merged, ref_annotation, merged_TPM, analysis.unflt_counts).collect()
|
||||
count_transcripts_file = count_transcripts.out.seqkit_stats.collect()
|
||||
all_counts = merged_TPM.concat(merged, analysis.flt_counts, analysis.gene_counts)
|
||||
all_counts = merged_TPM.concat(analysis.flt_counts, analysis.gene_counts)
|
||||
emit:
|
||||
all_de = de_report
|
||||
count_transcripts = count_transcripts_file
|
||||
|
||||
Loading…
Reference in New Issue
Block a user