Account for stringtie multigene transcript artefacts
This commit is contained in:
parent
112e87ad94
commit
0a5c39cdf4
@ -12,6 +12,8 @@ variables:
|
|||||||
--direct_rna --minimap2_index_opts '-k 15' --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \
|
--direct_rna --minimap2_index_opts '-k 15' --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression/sample_sheet.csv \
|
||||||
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
||||||
CI_FLAVOUR: "new"
|
CI_FLAVOUR: "new"
|
||||||
|
PYTEST_CONTAINER_NAME: "wf-common"
|
||||||
|
PYTEST_CONTAINER_CONFIG_KEY: "common_sha"
|
||||||
|
|
||||||
macos-run:
|
macos-run:
|
||||||
# Let's avoid those ARM64 runners for now
|
# Let's avoid those ARM64 runners for now
|
||||||
@ -146,7 +148,7 @@ docker-run:
|
|||||||
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna.gz \
|
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna.gz \
|
||||||
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.gff.gz \
|
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.gff.gz \
|
||||||
--direct_rna --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_rna.fna.gz \
|
--direct_rna --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_rna.fna.gz \
|
||||||
--transcriptome_assembly false --sample_sheet test_data/sample_sheet.csv \
|
--sample_sheet test_data/sample_sheet.csv \
|
||||||
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
||||||
NF_IGNORE_PROCESSES: >
|
NF_IGNORE_PROCESSES: >
|
||||||
preprocess_reads,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,
|
preprocess_reads,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,
|
||||||
@ -160,7 +162,7 @@ docker-run:
|
|||||||
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz \
|
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz \
|
||||||
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.109.gtf.gz \
|
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.109.gtf.gz \
|
||||||
--direct_rna --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.cdna.all.fa.gz \
|
--direct_rna --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression/Homo_sapiens.GRCh38.cdna.all.fa.gz \
|
||||||
--transcriptome_assembly false --sample_sheet test_data/sample_sheet.csv \
|
--sample_sheet test_data/sample_sheet.csv \
|
||||||
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
||||||
NF_IGNORE_PROCESSES: >
|
NF_IGNORE_PROCESSES: >
|
||||||
preprocess_reads,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,
|
preprocess_reads,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,
|
||||||
@ -174,7 +176,7 @@ docker-run:
|
|||||||
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_mouse/GRCm39.genome.fa.gz \
|
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_mouse/GRCm39.genome.fa.gz \
|
||||||
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_mouse/gencode.vM33.annotation.gtf \
|
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_mouse/gencode.vM33.annotation.gtf \
|
||||||
--direct_rna --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_mouse/gencode.vM33.transcripts.fa.gz \
|
--direct_rna --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_mouse/gencode.vM33.transcripts.fa.gz \
|
||||||
--transcriptome_assembly false --sample_sheet ${CI_PROJECT_NAME}/data/differential_expression_mouse/sample_sheet.csv \
|
--sample_sheet ${CI_PROJECT_NAME}/data/differential_expression_mouse/sample_sheet.csv \
|
||||||
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
||||||
NF_IGNORE_PROCESSES: >
|
NF_IGNORE_PROCESSES: >
|
||||||
preprocess_reads,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,
|
preprocess_reads,faidx,gz_faidx,merge_transcriptomes,assemble_transcripts,decompress_annotation,
|
||||||
@ -215,7 +217,7 @@ docker-run:
|
|||||||
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna.gz \
|
--ref_genome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.fna.gz \
|
||||||
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.gff.gz \
|
--ref_annotation ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_genomic.gff.gz \
|
||||||
--direct_rna --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_rna.fna.gz \
|
--direct_rna --ref_transcriptome ${CI_PROJECT_NAME}/data/differential_expression_ncbi/GCF_000001405.40_GRCh38.p14_rna.fna.gz \
|
||||||
--transcriptome_assembly false --sample_sheet test_data/sample_sheet.csv \
|
--sample_sheet test_data/sample_sheet.csv \
|
||||||
--igv \
|
--igv \
|
||||||
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
-c ${CI_PROJECT_NAME}/data/demo.nextflow.config"
|
||||||
NF_IGNORE_PROCESSES: >
|
NF_IGNORE_PROCESSES: >
|
||||||
|
|||||||
@ -4,6 +4,10 @@ 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/),
|
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).
|
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
|
||||||
|
|
||||||
|
## [Unreleased]
|
||||||
|
### Fixed
|
||||||
|
- Bug that led to incorrect gene_id being assigned in the DE plots.
|
||||||
|
|
||||||
## [v1.5.0]
|
## [v1.5.0]
|
||||||
### Updated
|
### Updated
|
||||||
- Workflow report updated to use `ezcharts`.
|
- Workflow report updated to use `ezcharts`.
|
||||||
|
|||||||
@ -1,6 +1,5 @@
|
|||||||
#!/usr/bin/env python
|
#!/usr/bin/env python
|
||||||
"""Create de report section."""
|
"""Create de report section."""
|
||||||
|
|
||||||
import os
|
import os
|
||||||
|
|
||||||
from dominate.tags import h5, p
|
from dominate.tags import h5, p
|
||||||
@ -52,7 +51,7 @@ def create_summary_table(df):
|
|||||||
avg_acc, avg_mapq])
|
avg_acc, avg_mapq])
|
||||||
|
|
||||||
|
|
||||||
def dexseq_section(dexseq_file, id_dic, pval_thresh):
|
def dexseq_section(dexseq_file, tr_id_to_gene_name, tr_id_to_gene_id, pval_thresh):
|
||||||
"""Add gene isoforms table and plot."""
|
"""Add gene isoforms table and plot."""
|
||||||
h5("Differential Isoform usage")
|
h5("Differential Isoform usage")
|
||||||
p("""Table showing gene isoforms, ranked by adjusted
|
p("""Table showing gene isoforms, ranked by adjusted
|
||||||
@ -64,9 +63,17 @@ def dexseq_section(dexseq_file, id_dic, pval_thresh):
|
|||||||
|
|
||||||
dexseq_results = pd.read_csv(dexseq_file, sep='\t')
|
dexseq_results = pd.read_csv(dexseq_file, sep='\t')
|
||||||
dexseq_results.index.name = "gene_id:transcript_id"
|
dexseq_results.index.name = "gene_id:transcript_id"
|
||||||
# Replace gene id with more useful gene name where possible
|
|
||||||
|
# Replace any occurrences of stringtie-generated MSTRG gene ids with
|
||||||
|
# reference gene_ids.
|
||||||
dexseq_results.index = dexseq_results.index.map(
|
dexseq_results.index = dexseq_results.index.map(
|
||||||
lambda x: str(id_dic.get(x.split(':')[0])) + ':' + str(x.split(':')[1]))
|
lambda ge_tr: str( # lookup gene_id from transcript_id [1]
|
||||||
|
f"{tr_id_to_gene_id.get(ge_tr.split(':')[1])}: {str(ge_tr.split(':')[1])}")
|
||||||
|
)
|
||||||
|
|
||||||
|
# Add gene name column.
|
||||||
|
dexseq_results.insert(0, "gene_name", dexseq_results.index.map(
|
||||||
|
lambda x: tr_id_to_gene_name.get(x.split(':')[1])))
|
||||||
|
|
||||||
DataTable.from_pandas(
|
DataTable.from_pandas(
|
||||||
dexseq_results.sort_values(by='pvalue', ascending=True), use_index=True)
|
dexseq_results.sort_values(by='pvalue', ascending=True), use_index=True)
|
||||||
@ -100,13 +107,12 @@ def dexseq_section(dexseq_file, id_dic, pval_thresh):
|
|||||||
EZChart(plot)
|
EZChart(plot)
|
||||||
|
|
||||||
|
|
||||||
def dtu_section(dtu_file, gt_dic, ge_dic):
|
def dtu_section(dtu_file, txid_to_gene_name):
|
||||||
"""Plot dtu section."""
|
"""Plot dtu section."""
|
||||||
dtu_results = pd.read_csv(dtu_file, sep='\t')
|
dtu_results = pd.read_csv(dtu_file, sep='\t')
|
||||||
dtu_results["gene_name"] = dtu_results["txID"].apply(
|
dtu_results["gene_name"] = dtu_results["txID"].apply(
|
||||||
lambda x: gt_dic.get(x))
|
lambda x: txid_to_gene_name.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_pvals = dtu_results.sort_values(by='gene', ascending=True)
|
||||||
raw("""Table showing gene and transcript identifiers
|
raw("""Table showing gene and transcript identifiers
|
||||||
and their FDR-corrected (False discovery rate - Benjamini-Hochberg) probabilities
|
and their FDR-corrected (False discovery rate - Benjamini-Hochberg) probabilities
|
||||||
@ -120,11 +126,10 @@ def dtu_section(dtu_file, gt_dic, ge_dic):
|
|||||||
raw("""View dtu_plots.pdf file to see plots of differential isoform usage""")
|
raw("""View dtu_plots.pdf file to see plots of differential isoform usage""")
|
||||||
|
|
||||||
|
|
||||||
def dge_section(dge_file, ids_dic, pval_thresh):
|
def dge_section(df, pval_thresh):
|
||||||
"""Create DGE table and MA plot."""
|
"""Create DGE table and MA plot."""
|
||||||
h5("Differential gene expression")
|
h5("Differential gene expression")
|
||||||
dge_results = pd.read_csv(dge_file, sep='\t')
|
df[['logFC', 'logCPM', 'F']] = df[
|
||||||
dge_results[['logFC', 'logCPM', 'F']] = dge_results[
|
|
||||||
['logFC', 'logCPM', 'F']].round(2)
|
['logFC', 'logCPM', 'F']].round(2)
|
||||||
|
|
||||||
p("""Table showing the genes from the edgeR analysis.
|
p("""Table showing the genes from the edgeR analysis.
|
||||||
@ -134,10 +139,9 @@ def dge_section(dge_file, ids_dic, pval_thresh):
|
|||||||
This table has not been
|
This table has not been
|
||||||
filtered for genes that satisfy statistical or magnitudinal thresholds""")
|
filtered for genes that satisfy statistical or magnitudinal thresholds""")
|
||||||
|
|
||||||
dge_results.index = dge_results.index.map(lambda x: ids_dic.get(x))
|
df = df.sort_values('FDR', ascending=True)
|
||||||
dge_results = dge_results.sort_values('FDR', ascending=True)
|
df.index.name = 'gene_id'
|
||||||
dge_results.index.name = 'Transcript'
|
DataTable.from_pandas(df, use_index=True)
|
||||||
DataTable.from_pandas(dge_results, use_index=True)
|
|
||||||
|
|
||||||
h5("Results of the edgeR Analysis.")
|
h5("Results of the edgeR Analysis.")
|
||||||
|
|
||||||
@ -150,15 +154,13 @@ def dge_section(dge_file, ids_dic, pval_thresh):
|
|||||||
(False discovery rate - Benjamini-Hochberg) p-value thresholds
|
(False discovery rate - Benjamini-Hochberg) p-value thresholds
|
||||||
defined are shaded as 'Up-' or 'Down-' regulated.
|
defined are shaded as 'Up-' or 'Down-' regulated.
|
||||||
""")
|
""")
|
||||||
|
df['sig'] = None
|
||||||
dge = pd.read_csv(dge_file, sep="\t")
|
df.loc[(df["logFC"] > 0) & (df['PValue'] < pval_thresh), 'sig'] = 'up'
|
||||||
dge['sig'] = None
|
df.loc[(df["logFC"] <= 0) & (df['PValue'] < pval_thresh), 'sig'] = 'down'
|
||||||
dge.loc[(dge["logFC"] > 0) & (dge['PValue'] < pval_thresh), 'sig'] = 'up'
|
df.loc[(df["PValue"] >= pval_thresh), 'sig'] = 'not_sig'
|
||||||
dge.loc[(dge["logFC"] <= 0) & (dge['PValue'] < pval_thresh), 'sig'] = 'down'
|
|
||||||
dge.loc[(dge["PValue"] >= pval_thresh), 'sig'] = 'not_sig'
|
|
||||||
|
|
||||||
plot = scatterplot(
|
plot = scatterplot(
|
||||||
data=dge, x='logCPM', y='logFC', hue='sig',
|
data=df, x='logCPM', y='logFC', hue='sig',
|
||||||
palette=['#E32636', '#7E8896', '#0A22DE'],
|
palette=['#E32636', '#7E8896', '#0A22DE'],
|
||||||
hue_order=['up', 'not_sig', 'down'], marker='circle')
|
hue_order=['up', 'not_sig', 'down'], marker='circle')
|
||||||
plot._fig.x_range.start = 10
|
plot._fig.x_range.start = 10
|
||||||
@ -188,53 +190,58 @@ def salmon_table(salmon_counts):
|
|||||||
|
|
||||||
|
|
||||||
def get_translations(gtf):
|
def get_translations(gtf):
|
||||||
"""Create dict with gene_name and gene_references."""
|
"""Create gene_and transcript id mappings.
|
||||||
|
|
||||||
|
Annotation can be stringtie-generated (GTF) or from the input
|
||||||
|
reference annotation (GTF or GFF3) and the various attributes can differ
|
||||||
|
"""
|
||||||
with open(gtf) as fh:
|
with open(gtf) as fh:
|
||||||
gene_txid = {}
|
txid_to_gene_name = {}
|
||||||
gene_geid = {}
|
gid_to_gene_name = {}
|
||||||
geid_gname = {}
|
tx_id_to_gene_id = {}
|
||||||
|
|
||||||
def get_feature(row, feature):
|
def get_feature(row, feature):
|
||||||
return row.split(feature)[1].split(
|
return row.split(feature)[1].split(
|
||||||
";")[0].replace('=', '').replace("\"", "").strip()
|
";")[0].replace('=', '').replace("\"", "").strip()
|
||||||
|
|
||||||
for i in fh:
|
for gff_entry in fh:
|
||||||
if i.startswith("#"):
|
# Process transcripts features only
|
||||||
|
if gff_entry.startswith("#") or gff_entry.split('\t')[2] != 'transcript':
|
||||||
continue
|
continue
|
||||||
# Different gtf/gff formats contain different attributes
|
# Different gtf/gff formats contain different attributes
|
||||||
# and different formating (eg. gene_name="xyz" or gene_name "xyz")
|
# and different formating (eg. gene_name="xyz" or gene_name "xyz")
|
||||||
gene_name = None
|
gene_name = gene_id = transcript_id = 'unknown'
|
||||||
for var_name in ["gene_name", "gene_id", "gene"]:
|
|
||||||
if var_name in i:
|
|
||||||
gene_name = get_feature(i, var_name)
|
|
||||||
break
|
|
||||||
|
|
||||||
if 'ref_gene_id' in i:
|
if 'ref_gene_id' in gff_entry:
|
||||||
gene_reference = get_feature(i, 'ref_gene_id')
|
# Favour ref_gene_id over gene_id. The latter can be multi-locus merged
|
||||||
elif 'gene_id' in i:
|
# genes from stringtie
|
||||||
gene_reference = get_feature(i, 'gene_id')
|
gene_id = get_feature(gff_entry, 'ref_gene_id')
|
||||||
|
elif 'gene_id' in gff_entry:
|
||||||
|
gene_id = get_feature(gff_entry, 'gene_id')
|
||||||
else:
|
else:
|
||||||
gene_reference = gene_name
|
gene_id = get_feature(gff_entry, 'gene')
|
||||||
if 'transcript_id' in i:
|
|
||||||
transcript_id = get_feature(i, 'transcript_id')
|
if 'transcript_id' in gff_entry:
|
||||||
|
transcript_id = get_feature(gff_entry, 'transcript_id')
|
||||||
|
|
||||||
|
if 'gene_name' in gff_entry:
|
||||||
|
gene_name = get_feature(gff_entry, 'gene_name')
|
||||||
else:
|
else:
|
||||||
transcript_id = "unknown"
|
# Fallback to gene_id if gene_name is not present
|
||||||
if 'gene_id' in i:
|
gene_name = gene_id
|
||||||
gene_id = get_feature(i, 'gene_id')
|
|
||||||
else:
|
txid_to_gene_name[transcript_id] = gene_name
|
||||||
gene_id = gene_name
|
tx_id_to_gene_id[transcript_id] = gene_id
|
||||||
gene_txid[transcript_id] = gene_name
|
gid_to_gene_name[gene_id] = gene_name
|
||||||
gene_geid[gene_id] = gene_reference
|
return txid_to_gene_name, tx_id_to_gene_id, gid_to_gene_name
|
||||||
geid_gname[gene_reference] = gene_name
|
|
||||||
return gene_txid, gene_geid, geid_gname
|
|
||||||
|
|
||||||
|
|
||||||
def de_section(
|
def de_section(
|
||||||
stringtie, dge, dexseq, dtu,
|
annotation, dge, dexseq, dtu,
|
||||||
tpm, report, filtered, unfiltered,
|
tpm, report, filtered, unfiltered,
|
||||||
gene_counts, aln_stats_dir, pval_threshold=0.01):
|
gene_counts, aln_stats_dir, pval_threshold=0.01):
|
||||||
"""Differential expression sections."""
|
"""Differential expression sections."""
|
||||||
with report.add_section("Differential expression", "DE"):
|
with (report.add_section("Differential expression", "DE")):
|
||||||
|
|
||||||
p("""This section shows differential gene expression
|
p("""This section shows differential gene expression
|
||||||
and differential isoform usage. Salmon was used to
|
and differential isoform usage. Salmon was used to
|
||||||
@ -256,39 +263,45 @@ def de_section(
|
|||||||
DataTable.from_pandas(alignment_summary_df, use_index=True)
|
DataTable.from_pandas(alignment_summary_df, use_index=True)
|
||||||
|
|
||||||
salmon_table(tpm)
|
salmon_table(tpm)
|
||||||
gene_txid, gene_name, geid_gname = get_translations(stringtie)
|
|
||||||
|
# Get translations for adding gene names to tables
|
||||||
|
(
|
||||||
|
txid_to_gene_name, txid_to_gene_id, gid_to_gene_name
|
||||||
|
) = get_translations(annotation)
|
||||||
|
|
||||||
# Add gene names columns to counts files and write out
|
# Add gene names columns to counts files and write out
|
||||||
# for publishing to user dir.
|
# for publishing to user dir.
|
||||||
df_dge = pd.read_csv(dge, sep='\t')
|
df_dge = pd.read_csv(dge, sep='\t')
|
||||||
df_dge.insert(0, 'gene_name', df_dge.index.map(lambda x: geid_gname.get(x)))
|
df_dge.insert(0, 'gene_name', df_dge.index.map(
|
||||||
|
lambda x: gid_to_gene_name.get(x)))
|
||||||
df_dge.to_csv('results_dge.tsv', index=True, index_label="gene_id", sep="\t")
|
df_dge.to_csv('results_dge.tsv', index=True, index_label="gene_id", sep="\t")
|
||||||
|
|
||||||
# write_dge(gene_counts, geid_gname, "all_gene_counts.tsv")
|
# write_dge(gene_counts, gid_to_gene_name, "all_gene_counts.tsv")
|
||||||
df_gene_counts = pd.read_csv(gene_counts, sep='\t')
|
df_gene_counts = pd.read_csv(gene_counts, sep='\t')
|
||||||
df_gene_counts.insert(
|
df_gene_counts.insert(
|
||||||
0, 'gene_name', df_gene_counts.index.map(lambda x: geid_gname.get(x)))
|
0, 'gene_name', df_gene_counts.index.map(
|
||||||
|
lambda x: gid_to_gene_name.get(x)))
|
||||||
df_gene_counts.to_csv(
|
df_gene_counts.to_csv(
|
||||||
'results_dge.tsv', index=True, index_label="gene_id", sep="\t")
|
'results_dge.tsv', index=True, index_label="gene_id", sep="\t")
|
||||||
|
|
||||||
df_filtered = pd.read_csv(filtered, sep='\t')
|
df_filtered = pd.read_csv(filtered, sep='\t')
|
||||||
df_filtered.insert(1, "gene_name", df_filtered.gene_id.map(
|
df_filtered.insert(1, "gene_name", df_filtered.gene_id.map(
|
||||||
lambda x: geid_gname.get(x)))
|
lambda x: gid_to_gene_name.get(x)))
|
||||||
df_filtered.to_csv(
|
df_filtered.to_csv(
|
||||||
'filtered_transcript_counts_with_genes.tsv', index=False, sep='\t')
|
'filtered_transcript_counts_with_genes.tsv', index=False, sep='\t')
|
||||||
|
|
||||||
df_unfiltered = pd.read_csv(unfiltered, sep='\t')
|
df_unfiltered = pd.read_csv(unfiltered, sep='\t')
|
||||||
df_unfiltered.insert(1, "gene_name", df_unfiltered.gene_id.map(
|
df_unfiltered.insert(1, "gene_name", df_unfiltered.gene_id.map(
|
||||||
lambda x: geid_gname.get(x)))
|
lambda x: gid_to_gene_name.get(x)))
|
||||||
df_unfiltered.to_csv(
|
df_unfiltered.to_csv(
|
||||||
'unfiltered_transcript_counts_with_genes.tsv', index=False, sep='\t')
|
'unfiltered_transcript_counts_with_genes.tsv', index=False, sep='\t')
|
||||||
|
|
||||||
df_tpm = pd.read_csv(tpm, sep='\t')
|
df_tpm = pd.read_csv(tpm, sep='\t')
|
||||||
df_tpm.insert(1, "gene_name", df_tpm.Reference.map(
|
df_tpm.insert(1, "gene_name", df_tpm.Reference.map(
|
||||||
lambda x: gene_txid.get(x)))
|
lambda x: txid_to_gene_name.get(x)))
|
||||||
df_tpm.to_csv("unfiltered_tpm_transcript_counts.tsv", index=False, sep='\t')
|
df_tpm.to_csv("unfiltered_tpm_transcript_counts.tsv", index=False, sep='\t')
|
||||||
|
|
||||||
# Add tables to report
|
# Add tables to report
|
||||||
dge_section(dge, gene_name, pval_threshold)
|
dge_section(df_dge, pval_threshold)
|
||||||
dexseq_section(dexseq, gene_name, pval_threshold)
|
dexseq_section(dexseq, txid_to_gene_name, txid_to_gene_id, pval_threshold)
|
||||||
dtu_section(dtu, gene_txid, gene_name)
|
dtu_section(dtu, txid_to_gene_name)
|
||||||
|
|||||||
@ -319,15 +319,16 @@ def de_section(report, de_report_dir, de_aln_stats_dir, pval_threshold):
|
|||||||
dexseq = de_report_dir / "results_dexseq.tsv"
|
dexseq = de_report_dir / "results_dexseq.tsv"
|
||||||
dge = de_report_dir / "results_dge.tsv"
|
dge = de_report_dir / "results_dge.tsv"
|
||||||
dtu = de_report_dir / "results_dtu_stageR.tsv"
|
dtu = de_report_dir / "results_dtu_stageR.tsv"
|
||||||
# GFF file can have gtf or gff extension
|
# GFF file can have gtf or gff extension.
|
||||||
stringtie = next(de_report_dir.glob("*.g*f*"))
|
# Will be the original (transcriptome_source=precomputed) or wf-assembled annotation
|
||||||
|
annotation = next(de_report_dir.glob("*.g*f*"))
|
||||||
tpm = de_report_dir / "unfiltered_tpm_transcript_counts.tsv"
|
tpm = de_report_dir / "unfiltered_tpm_transcript_counts.tsv"
|
||||||
filtered = de_report_dir / "filtered_transcript_counts_with_genes.tsv"
|
filtered = de_report_dir / "filtered_transcript_counts_with_genes.tsv"
|
||||||
unfiltered = de_report_dir / "unfiltered_transcript_counts_with_genes.tsv"
|
unfiltered = de_report_dir / "unfiltered_transcript_counts_with_genes.tsv"
|
||||||
gene_counts = de_report_dir / "all_gene_counts.tsv"
|
gene_counts = de_report_dir / "all_gene_counts.tsv"
|
||||||
# This will also add a gene name column to the above counts tsv files
|
# This will also add a gene name column to the above counts tsv files
|
||||||
de_plots.de_section(
|
de_plots.de_section(
|
||||||
stringtie=stringtie,
|
annotation=annotation,
|
||||||
dexseq=dexseq,
|
dexseq=dexseq,
|
||||||
dge=dge,
|
dge=dge,
|
||||||
dtu=dtu,
|
dtu=dtu,
|
||||||
|
|||||||
80
bin/workflow_glue/tests/test_de_plots.py
Normal file
80
bin/workflow_glue/tests/test_de_plots.py
Normal file
@ -0,0 +1,80 @@
|
|||||||
|
"""Test assign_barcodes."""
|
||||||
|
from pathlib import Path
|
||||||
|
|
||||||
|
import pytest
|
||||||
|
from workflow_glue.de_plots import get_translations
|
||||||
|
|
||||||
|
|
||||||
|
@pytest.fixture
|
||||||
|
def test_data(request):
|
||||||
|
"""Define data location fixture."""
|
||||||
|
return Path(request.config.getoption("--test_data")) / "workflow_glue"
|
||||||
|
|
||||||
|
|
||||||
|
@pytest.mark.parametrize(
|
||||||
|
'annotation_file,expected',
|
||||||
|
[
|
||||||
|
[
|
||||||
|
'MSTRG.11088.gtf',
|
||||||
|
dict(gid_to_gene_name={
|
||||||
|
'ENSG00000236051.7': 'MYCBP2-AS1',
|
||||||
|
'ENSG00000283208.2': 'ENSG00000283208',
|
||||||
|
'ENSG00000102805.16': 'CLN5',
|
||||||
|
'MSTRG.11088': 'MSTRG.11088'
|
||||||
|
},
|
||||||
|
txid_to_gene_name={
|
||||||
|
'ENST00000636183.2': 'CLN5',
|
||||||
|
'ENST00000636780.2': 'CLN5',
|
||||||
|
'ENST00000638147.2': 'ENSG00000283208',
|
||||||
|
'ENST00000637192.1': 'ENSG00000283208',
|
||||||
|
'ENST00000636737.1': 'MYCBP2-AS1',
|
||||||
|
'ENST00000450627.6': 'MYCBP2-AS1',
|
||||||
|
'MSTRG.11088.2': 'MSTRG.11088'
|
||||||
|
},
|
||||||
|
txid_to_gene_id={
|
||||||
|
'ENST00000636183.2': 'ENSG00000102805.16',
|
||||||
|
'MSTRG.11088.2': 'MSTRG.11088',
|
||||||
|
'ENST00000636780.2': 'ENSG00000102805.16',
|
||||||
|
'ENST00000638147.2': 'ENSG00000283208.2',
|
||||||
|
'ENST00000637192.1': 'ENSG00000283208.2',
|
||||||
|
'ENST00000636737.1': 'ENSG00000236051.7',
|
||||||
|
'ENST00000450627.6': 'ENSG00000236051.7'
|
||||||
|
})
|
||||||
|
|
||||||
|
],
|
||||||
|
# Small test to check that GFF3 works
|
||||||
|
[
|
||||||
|
'MSTRG.11088.gff3',
|
||||||
|
dict(gid_to_gene_name={
|
||||||
|
"ENSG00000290825.1": "DDX11L2",
|
||||||
|
"ENSG00000236397.3": "DDX11L2"
|
||||||
|
},
|
||||||
|
txid_to_gene_name={
|
||||||
|
"ENST00000456328.2": "DDX11L2",
|
||||||
|
"ENST00000437401.1": "DDX11L2"
|
||||||
|
},
|
||||||
|
txid_to_gene_id={
|
||||||
|
'ENST00000437401.1': 'ENSG00000236397.3',
|
||||||
|
'ENST00000456328.2': 'ENSG00000290825.1'
|
||||||
|
})
|
||||||
|
]
|
||||||
|
]
|
||||||
|
)
|
||||||
|
def test_get_translations(test_data, annotation_file, expected):
|
||||||
|
"""Test that correct feature identifiers are extracted from the annotation.
|
||||||
|
|
||||||
|
`stringtie --merge` can sometimes generate gene models that may span multiple
|
||||||
|
reference genes. Possibly related issue:
|
||||||
|
https://github.com/gpertea/stringtie/issues/217
|
||||||
|
This can lead to the original genes and transcripts being assigned to that
|
||||||
|
incorrectly-merged gene model. The test data contains such a gene model generated
|
||||||
|
from `stringtie --merge` but actually consists of multiple different genes.
|
||||||
|
|
||||||
|
|
||||||
|
"""
|
||||||
|
input_gtf = test_data / annotation_file
|
||||||
|
txid_to_gene_name, txid_to_gene_id, gid_to_gene_name = get_translations(input_gtf)
|
||||||
|
|
||||||
|
assert expected['gid_to_gene_name'] == gid_to_gene_name
|
||||||
|
assert expected['txid_to_gene_name'] == txid_to_gene_name
|
||||||
|
assert expected['txid_to_gene_id'] == txid_to_gene_id
|
||||||
@ -95,7 +95,7 @@ params {
|
|||||||
]
|
]
|
||||||
agent = null
|
agent = null
|
||||||
container_sha = "shad8671ea3a8ed52f2c0f40355e8eb5c6f00d2cbda"
|
container_sha = "shad8671ea3a8ed52f2c0f40355e8eb5c6f00d2cbda"
|
||||||
common_sha="shaf15f9d80aba72c20e3e71f84869619873a56b8af"
|
common_sha="shabadd33adae761be6f2d59c6ecfb44b19cf472cfc"
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|||||||
@ -162,8 +162,9 @@ workflow differential_expression {
|
|||||||
analysis = deAnalysis(sample_sheet, merged, ref_annotation)
|
analysis = deAnalysis(sample_sheet, merged, ref_annotation)
|
||||||
plotResults(analysis.flt_counts, analysis.stageR, sample_sheet)
|
plotResults(analysis.flt_counts, analysis.stageR, sample_sheet)
|
||||||
// Concat files required for making the report
|
// Concat files required for making the report
|
||||||
de_report = analysis.flt_counts.concat(analysis.gene_counts, analysis.dge, analysis.dexseq,
|
de_report = analysis.flt_counts.concat(
|
||||||
analysis.stageR, sample_sheet, merged, ref_annotation, merged_TPM, analysis.unflt_counts).collect()
|
analysis.gene_counts, analysis.dge, analysis.dexseq,
|
||||||
|
analysis.stageR, sample_sheet, merged, ref_annotation, merged_TPM, analysis.unflt_counts).collect()
|
||||||
// Concat files required to be output to user without any changes
|
// Concat files required to be output to user without any changes
|
||||||
de_outputs_concat = analysis.cpm.concat(plotResults.out.dtu_plots, analysis.dge_pdf, analysis.dge_tsv, analysis.dexseq,
|
de_outputs_concat = analysis.cpm.concat(plotResults.out.dtu_plots, analysis.dge_pdf, analysis.dge_tsv, analysis.dexseq,
|
||||||
analysis.dtu_gene, analysis.dtu_transcript, analysis.dtu_stageR, analysis.dtu_pdf, analysis.flt_counts, analysis.gene_counts, merged_TPM).collect()
|
analysis.dtu_gene, analysis.dtu_transcript, analysis.dtu_stageR, analysis.dtu_pdf, analysis.flt_counts, analysis.gene_counts, merged_TPM).collect()
|
||||||
|
|||||||
2
test_data/workflow_glue/MSTRG.11088.gff3
Normal file
2
test_data/workflow_glue/MSTRG.11088.gff3
Normal file
@ -0,0 +1,2 @@
|
|||||||
|
chr1 HAVANA transcript 11869 14409 . + . ID=ENST00000456328.2;Parent=ENSG00000290825.1;gene_id=ENSG00000290825.1;transcript_id=ENST00000456328.2;gene_type=lncRNA;gene_name=DDX11L2;transcript_type=lncRNA;transcript_name=DDX11L2-202;level=2;transcript_support_level=1;tag=basic,Ensembl_canonical;havana_transcript=OTTHUMT00000362751.1
|
||||||
|
chr2 HAVANA transcript 113599036 113601261 . - . ID=ENST00000437401.1;Parent=ENSG00000236397.3;gene_id=ENSG00000236397.3;transcript_id=ENST00000437401.1;gene_type=unprocessed_pseudogene;gene_name=DDX11L2;transcript_type=unprocessed_pseudogene;transcript_name=DDX11L2-201;level=2;transcript_support_level=NA;hgnc_id=HGNC:37103;ont=PGO:0000005;tag=basic,Ensembl_canonical;havana_gene=OTTHUMG00000047823.1;havana_transcript=OTTHUMT00000109036.1
|
||||||
34
test_data/workflow_glue/MSTRG.11088.gtf
Normal file
34
test_data/workflow_glue/MSTRG.11088.gtf
Normal file
@ -0,0 +1,34 @@
|
|||||||
|
chr13 StringTie transcript 76990660 77005117 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636183.2"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie exon 76990660 76992271 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636183.2"; exon_number "1"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie exon 76995063 76995228 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636183.2"; exon_number "2"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie exon 76995902 76996127 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636183.2"; exon_number "3"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie exon 77000458 77005117 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636183.2"; exon_number "4"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie transcript 76991729 77005117 1000 + . gene_id "MSTRG.11088"; transcript_id "MSTRG.11088.2";
|
||||||
|
chr13 StringTie exon 76991729 76991832 1000 + . gene_id "MSTRG.11088"; transcript_id "MSTRG.11088.2"; exon_number "1";
|
||||||
|
chr13 StringTie exon 76995063 76995228 1000 + . gene_id "MSTRG.11088"; transcript_id "MSTRG.11088.2"; exon_number "2";
|
||||||
|
chr13 StringTie exon 76995902 76996127 1000 + . gene_id "MSTRG.11088"; transcript_id "MSTRG.11088.2"; exon_number "3";
|
||||||
|
chr13 StringTie exon 77000458 77005117 1000 + . gene_id "MSTRG.11088"; transcript_id "MSTRG.11088.2"; exon_number "4";
|
||||||
|
chr13 StringTie transcript 76992044 77005117 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636780.2"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie exon 76992044 76992271 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636780.2"; exon_number "1"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie exon 76995063 76995228 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636780.2"; exon_number "2"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie exon 76995902 76996127 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636780.2"; exon_number "3"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie exon 76998043 76998085 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636780.2"; exon_number "4"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie exon 77000458 77005117 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636780.2"; exon_number "5"; gene_name "CLN5"; ref_gene_id "ENSG00000102805.16";
|
||||||
|
chr13 StringTie transcript 76992078 77078025 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000638147.2"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie exon 76992078 76992271 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000638147.2"; exon_number "1"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie exon 76995063 76995228 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000638147.2"; exon_number "2"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie exon 76995902 76996127 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000638147.2"; exon_number "3"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie exon 77075518 77075584 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000638147.2"; exon_number "4"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie exon 77076816 77078025 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000638147.2"; exon_number "5"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie transcript 76995915 77129717 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000637192.1"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie exon 76995915 76996127 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000637192.1"; exon_number "1"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie exon 77109648 77110102 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000637192.1"; exon_number "2"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie exon 77129147 77129717 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000637192.1"; exon_number "3"; gene_name "ENSG00000283208"; ref_gene_id "ENSG00000283208.2";
|
||||||
|
chr13 StringTie transcript 77026767 77078025 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636737.1"; gene_name "MYCBP2-AS1"; ref_gene_id "ENSG00000236051.7";
|
||||||
|
chr13 StringTie exon 77026767 77027122 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636737.1"; exon_number "1"; gene_name "MYCBP2-AS1"; ref_gene_id "ENSG00000236051.7";
|
||||||
|
chr13 StringTie exon 77075518 77075584 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636737.1"; exon_number "2"; gene_name "MYCBP2-AS1"; ref_gene_id "ENSG00000236051.7";
|
||||||
|
chr13 StringTie exon 77076816 77078025 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000636737.1"; exon_number "3"; gene_name "MYCBP2-AS1"; ref_gene_id "ENSG00000236051.7";
|
||||||
|
chr13 StringTie transcript 77075514 77087778 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000450627.6"; gene_name "MYCBP2-AS1"; ref_gene_id "ENSG00000236051.7";
|
||||||
|
chr13 StringTie exon 77075514 77075584 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000450627.6"; exon_number "1"; gene_name "MYCBP2-AS1"; ref_gene_id "ENSG00000236051.7";
|
||||||
|
chr13 StringTie exon 77076816 77076866 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000450627.6"; exon_number "2"; gene_name "MYCBP2-AS1"; ref_gene_id "ENSG00000236051.7";
|
||||||
|
chr13 StringTie exon 77087552 77087778 1000 + . gene_id "MSTRG.11088"; transcript_id "ENST00000450627.6"; exon_number "3"; gene_name "MYCBP2-AS1"; ref_gene_id "ENSG00000236051.7";
|
||||||
Loading…
Reference in New Issue
Block a user