@prefix rdf: . @prefix CommandLineBinding: . @prefix xsd: . @prefix DockerRequirement: . @prefix Workflow: . @prefix cwl: . @prefix sld: . @prefix rdfs: . @prefix ResourceRequirement: . @prefix CommandOutputBinding: . @prefix ns2: . @prefix ns1: . cwl:source . rdfs:comment "One or more TSV files belonging to set B" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--files_setB" ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . rdfs:comment "Tool in comparison" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--tool" ] ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--tool" ] ; sld:default "annoFuse" ; sld:type xsd:string , sld:null . rdfs:comment "One or more TSV files belonging to set B" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--file-b" ] ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--file-b" ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . a cwl:Workflow ; rdfs:comment "# GENCODE Benchmark\n\n## Introduction\n\nThe GENCODE benchmark workflow generates a series of plots and other files used to compare two sets of data. While made specifically for GENCODE comparisons,\nthe workflow will compare any two sets of KF RNAseq data. In total the workflow generates four plots; two are generated from the annoFuse inputs and two are\ngenerated from the RSEM gene results of the Kids First RNAseq pipeline. The first of these plots is a raw count comparison of the annoFuse fusions between \nthe two sets. For each sample, it shows the number of fusions that are shared or unique to each sample. The second of these plots is the first plot but\nrather than showing counts it shows percentage of total fusions. The third plot is a heatmap showing the correlation of RSEM TMP values among samples. The \nfourth and final plot is an RSEM TPM tSNE plot for all the samples. In additon to the plot the workflow also returns the a CSV containing the raw correlation\nvalues used to generate the heatmap as well as a series of files that detail the fusions that are unique to each Sample and Set.\n\nThese plots and other files are for evaluating the effect of going from one GENCODE version to another. The general idea is that greater differences between\nversions will result in:\n- more pronounced version-based tSNE clustering\n- lower same-sample heatmap correlation\n- greater fusion count differences\n\n## Usage\n\n### Inputs\n\n- `baseline_annofuses`: *.annoFuse_filter.tsv files from baseline samples\n- `baseline_rsems`: *.rsem.genes.results.gz files from baseline samples\n- `baseline_label`: Description of baseline samples\n- `call_annofuses`: *.annoFuse_filter.tsv files from call samples\n- `call_rsems`: *.rsem.genes.results.gz files from call samples\n- `call_label`: Description of samples\n- `target_genes`: TSV containing target ENSGs for heatmap and tsne plots. Recommendation: stable housekeeping gene list\n- `output_basename`: Basename for output files\n\n### Outputs\n\n- `heatmap`: Heatmap plot comparing baseline and call samples by RSEM gene TPM values for ENSGs in the target_genes file\n- `corr`: Heatmap CSV comparing baseline and call samples by RSEM gene TPM values for ENSGs in the target_genes file\n- `tsne`: tSNE plot clustering baseline and call samples by RSEM gene TMP values for ENSGs in the target_genes files\n- `fusion_comp_raw`: Stacked bar chart showing shared and unique fusions for each sample\n- `fusion_comp_pct`: Stacked bar chart showing percente of shared and unique fusions for each sample\n- `fusion_comp_uniqs`: Files with unique fusions for each sample x label\n" ; rdfs:label "GENCODE Bencharmking" ; ns2:original_cwlVersion "v1.2" ; ns1:sevenbridges.comcategories "CORR" , "GENCODE" , "FUSION" , "HEATMAP" , "RSEM" , "GENE" , "ANNOFUSE" , "RNASEQ" , "TSNE" ; ns1:sevenbridges.comlicense "Apache License 2.0" ; ns1:sevenbridges.compublisher "KFDRC" ; Workflow:steps , , , ; cwl:cwlVersion cwl:v1.2 ; cwl:inputs , , , , , , , ; cwl:outputs , , , , , ; cwl:requirements [ a cwl:StepInputExpressionRequirement ] ; cwl:requirements [ a cwl:InlineJavascriptRequirement ] . cwl:outputBinding [ CommandOutputBinding:glob "*.tsne_plot.png" ] ; sld:type cwl:File . rdfs:comment "TSV containing target ENSGs for heatmap and tsne plots. Recommendation: stable housekeeping gene list" ; sld:type cwl:File . cwl:source . cwl:in , , , , , , ; cwl:out , ; cwl:run . cwl:outputSource ; sld:type cwl:File . rdfs:comment "*.annoFuse_filter.tsv files from Set B" ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . rdfs:comment "CPUs to allocate to this task" ; sld:default 1 ; sld:type xsd:int , sld:null . cwl:outputSource ; sld:type cwl:File . rdfs:comment "Basename for output files" ; sld:type xsd:string . rdfs:comment "TSV containing ENSG names on which the TSNE will be performed" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--gene_list_path" ] ; sld:type cwl:File . cwl:source . rdfs:comment "Label to use for set A in the legend and filenames" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--setA_name" ] ; sld:type xsd:string , sld:null . cwl:valueFrom "annoFuse" . rdfs:comment "Description of call files" ; sld:type xsd:string . rdfs:comment "Label to use for set A in the legend and filenames" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--label-a" ] ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--label-a" ] ; sld:type xsd:string , sld:null . cwl:outputSource ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . cwl:source . cwl:outputSource ; sld:type cwl:File . cwl:source . cwl:source . rdfs:comment "GB of RAM to allocate to this task" ; sld:default 8 ; sld:type xsd:int , sld:null . rdfs:comment "One or more TSV files belonging to set A" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--files_setA" ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . cwl:source ; cwl:valueFrom "$(self).annofuse.percent.png" . cwl:source . rdfs:comment "Label to use for set B in the legend and filenames" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--setB_name" ] ; sld:type xsd:string , sld:null . rdfs:comment "One or more TSV files belonging to set A" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--file-a" ] ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--file-a" ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . rdfs:comment "Basename for output files (no extension). If omitted, defaults to v." ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--output_basename" ] ; sld:type xsd:string , sld:null . cwl:source . cwl:source . cwl:outputBinding [ CommandOutputBinding:glob "unique/*" ] ; cwl:outputBinding [ CommandOutputBinding:glob "unique/*" ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . rdfs:comment "Description of baseline files" ; sld:type xsd:string . rdfs:comment "Basename for output files (no extension). If omitted, defaults to v." ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--output_basename" ] ; sld:type xsd:string , sld:null . rdfs:comment "CPUs to allocate to this task" ; sld:default 1 ; sld:type xsd:int , sld:null . cwl:source . cwl:source . a cwl:CommandLineTool ; ns2:original_cwlVersion "v1.2" ; cwl:arguments ( [ CommandLineBinding:position 1 ; CommandLineBinding:shellQuote false ; cwl:valueFrom "Rscript corr_heatmap.R" ] ) ; cwl:baseCommand () ; cwl:cwlVersion cwl:v1.2 ; cwl:inputs , , , , , , , ; cwl:outputs , ; cwl:requirements [ a cwl:InitialWorkDirRequirement ; cwl:listing [ cwl:entry "#!/usr/bin/env Rscript\n#'\n#' Sample–sample correlation analysis using RSEM TPM expression\n#'\n#' This script computes and visualizes sample–sample correlations from\n#' gene-level TPM expression values produced by RSEM. Expression values\n#' are log2-transformed, filtered to retain expressed genes, and used to\n#' calculate Pearson correlations between samples.\n#'\n#' The resulting correlation matrix is clustered using correlation\n#' distance (1 − Pearson correlation) and visualized as a heatmap.\n#' The script is intended for exploratory quality control and comparison\n#' of two sets of RNA-seq samples.\n#'\n#' Main steps:\n#' 1. Read RSEM gene-level TPM files (TSV or TSV.GZ).\n#' 2. Normalize gene identifiers by removing version suffixes.\n#' 3. Subset to a provided gene list.\n#' 4. Filter low-expression genes.\n#' 5. Log2-transform TPM values.\n#' 6. Compute sample–sample Pearson correlations.\n#' 7. Cluster samples using correlation distance and generate a heatmap.\n#'\n#' Notes:\n#' - Correlation-based similarity is robust to global expression shifts\n#' and is appropriate for cross-sample comparisons.\n#' - Clustering uses 1 − correlation as a distance metric to ensure\n#' consistency between similarity and hierarchical clustering.\n#'\n#' Usage:\n#' sample_correlation_heatmap.R \\\\\n#' --gene_list_path genes.tsv \\\\\n#' --setA_name SetA --files_setA A1.tsv.gz A2.tsv.gz \\\\\n#' --setB_name SetB --files_setB B1.tsv.gz B2.tsv.gz \\\\\n#' --output_basename SetA_vs_SetB\n#'\n#' Outputs:\n#' - .sample_correlation_heatmap.csv\n#' - .sample_correlation_heatmap.pdf\n#' - .removed_genes_tpm.csv\n#'\n#' Requirements:\n#' R >= 4.0\n#' Packages: readr, dplyr, tidyr, purrr, pheatmap, argparse\n#'\n\nsuppressPackageStartupMessages({\n library(readr)\n library(dplyr)\n library(tidyr)\n library(tibble)\n library(purrr)\n library(pheatmap)\n library(argparse)\n})\n\n#' Read RSEM TPM file (optionally gzipped)\n#'\n#' @param path Path to a TSV or TSV.GZ file with columns `gene_id` and `TPM`\n#' @return A tibble with columns `gene_id` and `TPM`\nread_rsem_tpm_gz <- function(path) {\n read_tsv(path, show_col_types = FALSE, progress = FALSE) %>%\n select(gene_id, TPM)\n}\n\n# ---- Args ----\np <- ArgumentParser()\np$add_argument(\"--gene_list_path\", required = TRUE)\n\np$add_argument(\"--setA_name\", default = \"setA\")\np$add_argument(\"--files_setA\", nargs = \"+\", required = TRUE)\n\np$add_argument(\"--setB_name\", default = \"setB\")\np$add_argument(\"--files_setB\", nargs = \"+\", required = TRUE)\n\np$add_argument(\"--output_basename\", default = NULL,\n help = \"Basename for output files (no extension).\")\n\nargs <- p$parse_args()\n\ngene_list_path <- args$gene_list_path\nsetA_name <- args$setA_name\nsetB_name <- args$setB_name\nout_basename <- args$output_basename\nif (is.null(out_basename) || !nzchar(out_basename)) {\n out_basename <- paste0(setA_name, \"v\", setB_name)\n}\n\nfiles_setA <- setNames(args$files_setA, paste0(\"A\", seq_along(args$files_setA)))\nfiles_setB <- setNames(args$files_setB, paste0(\"B\", seq_along(args$files_setB)))\n\nfiles <- c(files_setA, files_setB)\n# ---- END Args ----\n\n\n# Long -> wide using a versionless gene key\ntpm_long <- imap_dfr(files, ~ read_rsem_tpm_gz(.x) %>% mutate(sample = .y)) %>%\n mutate(\n gene_ensg = sub(\"_.*$\", \"\", gene_id), # ENSG...\n gene_ensg = sub(\"\\\\..*$\", \"\", gene_ensg), # ENSG... (no version)\n ) %>%\n select(gene_ensg, sample, TPM)\n\n# If the same gene_ensg appears more than once per sample, collapse it (usually not needed, but safe)\ntpm_long <- tpm_long %>%\n group_by(gene_ensg, sample) %>%\n summarise(TPM = sum(TPM, na.rm = TRUE), .groups = \"drop\")\n\ntpm_mat <- tpm_long %>%\n pivot_wider(names_from = sample, values_from = TPM, values_fill = 0) %>%\n column_to_rownames(\"gene_ensg\") %>%\n as.matrix()\n\ngenes_tbl <- readr::read_tsv(\n gene_list_path,\n col_names = c(\"gene_id\", \"symbol\"),\n show_col_types = FALSE\n)\n\n# Keep list as Ensembl-without-version (robust)\ngenes_keep <- unique(sub(\"\\\\..*$\", \"\", genes_tbl$gene_id))\n\n# Extract Ensembl-without-version from tpm_mat rownames like \"ENSG....15_TSPAN6\"\ntpm_ensg <- sub(\"_.*$\", \"\", rownames(tpm_mat)) # \"ENSG....15\"\ntpm_ensg <- sub(\"\\\\..*$\", \"\", tpm_ensg) # \"ENSG....\"\n\n# Subset\ntpm_mat <- tpm_mat[tpm_ensg %in% genes_keep, , drop = FALSE]\n\n# Optional: filter low-expression genes (recommended for correlations)\nkeep <- rowSums(tpm_mat >= 1) >= 2\n\nremoved_rows <- tpm_mat[!keep, , drop = FALSE]\nwrite.csv(\n removed_rows,\n paste0(out_basename, \".removed_genes_tpm.csv\"),\n quote = FALSE\n)\n\ntpm_mat <- tpm_mat[keep, , drop = FALSE]\n\n# Samples x genes, log-transform\nx <- t(log2(tpm_mat + 1))\ndim(x)\n\n# Sample–sample correlation\nsystem.time(cors <- cor(t(x), method = \"pearson\", use = \"pairwise.complete.obs\"))\nwrite.csv(\n cors,\n paste0(out_basename, \".sample_correlation_heatmap.csv\"),\n quote = FALSE\n)\n\n# Optional annotation (edit the rule to match your sample names)\nann <- data.frame(Set = ifelse(grepl(\"^A\", rownames(cors)), setA_name, setB_name))\nrownames(ann) <- rownames(cors)\n\n# ---- CORRELATION DISTANCE FOR CLUSTERING ----\ndist_corr <- as.dist(1 - cors)\n\n# Hierarchical clustering can vary slightly due to ties. Set seed for reproducibility\nset.seed(1)\n# Pearson correlation on log2(TPM+1); clustering uses 1 - correlation distance\npheatmap(\n cors,\n clustering_distance_rows = dist_corr,\n clustering_distance_cols = dist_corr,\n annotation_col = ann,\n annotation_row = ann,\n main = \"Sample correlation (Pearson) on log2(TPM+1)\",\n filename = paste0(out_basename, \".sample_correlation_heatmap.pdf\"),\n width = 8,\n height = 8\n)\n\n" ; cwl:entryname "corr_heatmap.R" ] ] ; cwl:requirements [ a cwl:ResourceRequirement ; ResourceRequirement:coresMin "$(inputs.cpu)" ; ResourceRequirement:ramMin "$(inputs.ram * 1000)" ] ; cwl:requirements [ a cwl:DockerRequirement ; DockerRequirement:dockerPull "pgc-images.sbgenomics.com/danmiller/plot_gencode_rnaseq_r:4.5.0" ] ; cwl:requirements [ a cwl:InlineJavascriptRequirement ] ; cwl:requirements [ a cwl:ShellCommandRequirement ] . cwl:source . cwl:valueFrom "$(true)" . cwl:source . cwl:outputSource ; sld:type cwl:File . cwl:outputBinding [ CommandOutputBinding:glob "*.sample_correlation_heatmap.csv" ] ; sld:type cwl:File . rdfs:comment "*.rsem.genes.results.gz files from Set B" ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . cwl:in , , , , , , ; cwl:out , ; cwl:run . cwl:source . rdfs:comment "GB of RAM to allocate to this task" ; sld:default 8 ; sld:type xsd:int , sld:null . rdfs:comment "Label to use for set A in the legend and filenames" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--setA_name" ] ; sld:type xsd:string , sld:null . cwl:source . cwl:source . cwl:outputBinding [ CommandOutputBinding:glob "*.sample_correlation_heatmap.pdf" ] ; sld:type cwl:File . rdfs:comment "Output figure path (e.g., overlap.png, overlap.pdf, overlap.svg)" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--out" ] ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--out" ] ; sld:type xsd:string . rdfs:comment "One or more TSV files belonging to set A" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--files_setA" ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . rdfs:comment "*.annoFuse_filter.tsv files from baseline" ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . cwl:source . rdfs:comment "CPUs to allocate to this task" ; sld:default 1 ; sld:type xsd:int , sld:null . rdfs:comment "*.rsem.genes.results.gz files from baseline" ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . cwl:source . a cwl:CommandLineTool ; ns2:original_cwlVersion "v1.2" ; cwl:arguments ( [ CommandLineBinding:position 1 ; CommandLineBinding:shellQuote false ; cwl:valueFrom "Rscript tsne_plot.R" ] ) ; cwl:baseCommand () ; cwl:cwlVersion cwl:v1.2 ; cwl:inputs , , , , , , , ; cwl:outputs ; cwl:requirements [ a cwl:DockerRequirement ; DockerRequirement:dockerPull "pgc-images.sbgenomics.com/danmiller/plot_gencode_rnaseq_r:4.5.0" ] ; cwl:requirements [ a cwl:InlineJavascriptRequirement ] ; cwl:requirements [ a cwl:InitialWorkDirRequirement ; cwl:listing [ cwl:entry "#!/usr/bin/env Rscript\n#'\n#' t-SNE visualization of RNA-seq samples using RSEM TPM values\n#'\n#' This script performs an exploratory t-SNE analysis of RNA-seq samples\n#' based on gene-level TPM expression values produced by RSEM. Expression\n#' values are log2-transformed, filtered to retain expressed genes, scaled\n#' per gene, and embedded using t-SNE with Euclidean distance.\n#'\n#' The script is designed for sample-level quality control and comparison\n#' between two sets of samples. It supports gzip-compressed input files,\n#' robust gene ID normalization, and produces a publication-ready t-SNE\n#' scatter plot.\n#'\n#' Main steps:\n#' 1. Read RSEM gene-level TPM files (TSV or TSV.GZ).\n#' 2. Normalize gene identifiers by removing version suffixes.\n#' 3. Filter low-expression and zero-variance genes.\n#' 4. Log2-transform and scale expression values per gene.\n#' 5. Compute a t-SNE embedding of samples.\n#' 6. Generate a labeled t-SNE plot colored by sample set.\n#'\n#' Notes:\n#' - t-SNE is an exploratory visualization; distances and clusters should\n#' not be over-interpreted, especially for small sample sizes.\n#' - Euclidean distance is used on scaled log2(TPM+1) values.\n#'\n#' Usage:\n#' tsne_samples.R \\\\\n#' --gene_list_path genes.tsv \\\\\n#' --setA_name SetA --files_setA A1.tsv.gz A2.tsv.gz \\\\\n#' --setB_name SetB --files_setB B1.tsv.gz B2.tsv.gz \\\\\n#' --output_basename SetA_vs_SetB\n#'\n#' Outputs:\n#' - .tsne_plot.png\n#'\n#' Requirements:\n#' R >= 4.0, packages: readr, dplyr, tidyr, purrr, Rtsne, ggplot2, ggrepel\n#'\n\nsuppressPackageStartupMessages({\n library(readr)\n library(dplyr)\n library(tidyr)\n library(tibble)\n library(purrr)\n library(Rtsne)\n library(ggplot2)\n library(ggrepel)\n library(argparse)\n})\n\n\n#' Read an RSEM TPM file (optionally gzipped)\n#'\n#' Reads an RSEM gene-level expression file containing TPM values.\n#' The file may be plain TSV or gzip-compressed (.gz). Only the\n#' `gene_id` and `TPM` columns are retained.\n#'\n#' @param path Character string giving the path to a TSV or TSV.GZ file\n#' with columns `gene_id` and `TPM`.\n#'\n#' @return A tibble with two columns:\n#' \\describe{\n#' \\item{gene_id}{Gene identifier (e.g., Ensembl ID)}\n#' \\item{TPM}{Transcripts Per Million values}\n#' }\n#'\n#' @examples\n#' \\dontrun{\n#' tpm <- read_rsem_tpm_gz(\"sample.rsem.genes.results.gz\")\n#' head(tpm)\n#' }\nread_rsem_tpm_gz <- function(path) {\n read_tsv(path, show_col_types = FALSE, progress = FALSE) %>%\n select(gene_id, TPM)\n}\n\n# ---- Args ----\np <- ArgumentParser()\np$add_argument(\"--gene_list_path\", required = TRUE)\n\np$add_argument(\"--setA_name\", default = \"setA\")\np$add_argument(\"--files_setA\", nargs = \"+\", required = TRUE)\n\np$add_argument(\"--setB_name\", default = \"setB\")\np$add_argument(\"--files_setB\", nargs = \"+\", required = TRUE)\n\np$add_argument(\"--output_basename\", default = NULL,\n help = \"Basename for output files (no extension).\")\n\nargs <- p$parse_args()\n\ngene_list_path <- args$gene_list_path\nsetA_name <- args$setA_name\nsetB_name <- args$setB_name\nout_basename <- args$output_basename\nif (is.null(out_basename) || !nzchar(out_basename)) {\n out_basename <- paste0(setA_name, \"v\", setB_name)\n}\n\nfiles_setA <- setNames(args$files_setA, paste0(\"A\", seq_along(args$files_setA)))\nfiles_setB <- setNames(args$files_setB, paste0(\"B\", seq_along(args$files_setB)))\n\nfiles <- c(files_setA, files_setB)\n# ---- END Args ----\n\n# Long -> wide using a versionless gene key\ntpm_long <- imap_dfr(files, ~ read_rsem_tpm_gz(.x) %>% mutate(sample = .y)) %>%\n mutate(\n gene_ensg = sub(\"_.*$\", \"\", gene_id), # ENSG...\n gene_ensg = sub(\"\\\\..*$\", \"\", gene_ensg), # ENSG... (no version)\n ) %>%\n select(gene_ensg, sample, TPM)\n\n# If the same gene_ensg appears more than once per sample, collapse it (usually not needed, but safe)\ntpm_long <- tpm_long %>%\n group_by(gene_ensg, sample) %>%\n summarise(TPM = sum(TPM, na.rm = TRUE), .groups = \"drop\")\n\ntpm_mat <- tpm_long %>%\n pivot_wider(names_from = sample, values_from = TPM, values_fill = 0) %>%\n column_to_rownames(\"gene_ensg\") %>%\n as.matrix()\n\ngenes_tbl <- readr::read_tsv(\n gene_list_path,\n col_names = c(\"gene_id\", \"symbol\"),\n show_col_types = FALSE\n)\n\n# Keep list as Ensembl-without-version (robust)\ngenes_keep <- unique(sub(\"\\\\..*$\", \"\", genes_tbl$gene_id))\n\n# Extract Ensembl-without-version from tpm_mat rownames like \"ENSG....15_TSPAN6\"\ntpm_ensg <- sub(\"_.*$\", \"\", rownames(tpm_mat)) # \"ENSG....15\"\ntpm_ensg <- sub(\"\\\\..*$\", \"\", tpm_ensg) # \"ENSG....\"\n\n# Subset\ntpm_mat <- tpm_mat[tpm_ensg %in% genes_keep, , drop = FALSE]\n\n# Filter low-expression genes (tune if desired)\nkeep <- rowSums(tpm_mat >= 1) >= 2\ntpm_mat <- tpm_mat[keep, , drop = FALSE]\n\n# Samples x genes, log-transform\nx <- t(log2(tpm_mat + 1))\n\n# remove zero-variance genes BEFORE scaling (prevents NaN)\nnzv <- apply(x, 2, sd) > 0\nx <- x[, nzv, drop = FALSE]\n\n# scale genes\nx <- scale(x)\n\n# Warn for small groups\nif (nrow(x) < 6) {\n warning(\"Very small sample size for t-SNE; interpretation may be unstable.\")\n}\n\n\n# final safety: remove any remaining non-finite values\nok <- is.finite(x)\nif (any(!ok)) {\n x <- x[, colSums(!ok) == 0, drop = FALSE]\n}\n\n# t-SNE (perplexity must be small for ~10 samples)\nset.seed(1)\nperp <- max(2, min(30, floor((nrow(x) - 1) / 3)))\n\n# t-SNE uses Euclidean distance on scaled log2(TPM+1)\nif (any(duplicated(x))) {\n warning(\"Duplicate samples detected in expression matrix.\")\n}\nts <- Rtsne(x, perplexity = perp, pca = TRUE, check_duplicates = FALSE)\n\nemb <- as.data.frame(ts$Y)\ncolnames(emb) <- c(\"tSNE1\", \"tSNE2\")\nemb$sample <- rownames(x)\nemb$set <- ifelse(emb$sample %in% names(files_setA), setA_name, setB_name)\n\n\nplt <- ggplot(emb, aes(tSNE1, tSNE2, color = set, label = sample)) +\n geom_point(size = 3) +\n geom_text_repel(size = 3, max.overlaps = Inf) +\n theme_classic() +\n labs(title = \"t-SNE of samples using RSEM TPM\",\n subtitle = \"log2(TPM+1), scaled per gene\")\n\nggsave(paste0(out_basename, \".tsne_plot.png\"),\n plot = plt, width = 7, height = 5, dpi = 300)\n" ; cwl:entryname "tsne_plot.R" ] ] ; cwl:requirements [ a cwl:ShellCommandRequirement ] ; cwl:requirements [ a cwl:ResourceRequirement ; ResourceRequirement:coresMin "$(inputs.cpu)" ; ResourceRequirement:ramMin "$(inputs.ram * 1000)" ] . cwl:in , , , , , ; cwl:out ; cwl:run . cwl:in , , , , , ; cwl:out , ; cwl:run . cwl:source ; cwl:valueFrom "$(self).annofuse.raw_count.png" . cwl:outputSource ; sld:type cwl:File . cwl:valueFrom "annoFuse" . cwl:source . rdfs:comment "One or more TSV files belonging to set B" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--files_setB" ] ; sld:type [ sld:items cwl:File ; sld:type sld:array ] . rdfs:comment "GB of RAM to allocate to this task" ; sld:default 8 ; sld:type xsd:int , sld:null . rdfs:comment "Label to use for set B in the legend and filenames" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--setB_name" ] ; sld:type xsd:string , sld:null . rdfs:comment "Plot title" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--title" ] ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--title" ] ; sld:type xsd:string , sld:null . rdfs:comment "Stack the bars by percentage rather than raw count" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--percent" ] ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--percent" ] ; sld:type xsd:boolean , sld:null . rdfs:comment "Label to use for set B in the legend and filenames" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--label-b" ] ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--label-b" ] ; sld:type xsd:string , sld:null . cwl:outputBinding [ CommandOutputBinding:glob "$(inputs.output_filename)" ] ; cwl:outputBinding [ CommandOutputBinding:glob "$(inputs.output_filename)" ] ; sld:type cwl:File . cwl:source . a cwl:CommandLineTool ; ns2:original_cwlVersion "v1.2" ; cwl:arguments ( [ CommandLineBinding:position 1 ; CommandLineBinding:shellQuote false ; cwl:valueFrom "python compare_sets.py --dump-unique-dir unique" ] ) ; cwl:arguments ( [ CommandLineBinding:position 1 ; CommandLineBinding:shellQuote false ; cwl:valueFrom "python compare_sets.py --dump-unique-dir unique" ] ) ; cwl:baseCommand () ; cwl:cwlVersion cwl:v1.2 ; cwl:inputs , , , , , , , , , ; cwl:outputs , ; cwl:requirements [ a cwl:InlineJavascriptRequirement ] ; cwl:requirements [ a cwl:ResourceRequirement ; ResourceRequirement:coresMin "$(inputs.cpu)" ; ResourceRequirement:ramMin "$(inputs.ram * 1000)" ] ; cwl:requirements [ a cwl:ShellCommandRequirement ] ; cwl:requirements [ a cwl:DockerRequirement ; DockerRequirement:dockerPull "pgc-images.sbgenomics.com/danmiller/plot_gencode_rnaseq_python:3.13.0" ] ; cwl:requirements [ a cwl:InitialWorkDirRequirement ; cwl:listing [ cwl:entry "#!/usr/bin/env python3\nimport argparse\nimport re\nfrom pathlib import Path\n\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\nSAMPLE_COL = \"Sample\"\nLEFT_COL = \"LeftBreakpoint\"\nRIGHT_COL = \"RightBreakpoint\"\n\n\ndef safe_name(s):\n \"\"\"Sanitize a string for safe use in filenames.\n\n Replaces characters that are not alphanumeric, underscore, dot,\n or hyphen with underscores.\n\n Args:\n s (Any): Input value to sanitize (will be converted to string).\n\n Returns:\n str: A filename-safe version of the input string.\n \"\"\"\n return re.sub(r'[^A-Za-z0-9_.-]', '_', str(s))\n\n\ndef load_files_with_key(paths):\n \"\"\"Load multiple TSV files and add a breakpoint-pair key column.\n\n Each input file is read as a tab-separated file with all columns\n treated as strings. A new column named ``key`` is created by\n concatenating ``LeftBreakpoint`` and ``RightBreakpoint`` with\n a pipe character (\"|\").\n\n Args:\n paths (Iterable[str | Path]): Paths to TSV files to load.\n\n Returns:\n pandas.DataFrame: Concatenated DataFrame containing all rows\n from the input files, with an added ``key`` column. If no\n files are provided, an empty DataFrame with ``Sample`` and\n ``key`` columns is returned.\n\n Raises:\n ValueError: If any input file is missing the required\n ``LeftBreakpoint`` or ``RightBreakpoint`` columns.\n \"\"\"\n dfs = []\n for p in paths:\n df = pd.read_csv(p, sep=\"\\t\", dtype=str)\n # ensure breakpoint columns exist\n if LEFT_COL not in df.columns or RIGHT_COL not in df.columns:\n raise ValueError(\n f\"File {p} missing required column(s): {LEFT_COL}, {RIGHT_COL}\"\n )\n df = df.copy()\n df[\"key\"] = df[LEFT_COL].str.cat(df[RIGHT_COL], sep=\"|\")\n dfs.append(df)\n\n if not dfs:\n return pd.DataFrame(columns=[SAMPLE_COL, \"key\"])\n\n return pd.concat(dfs, ignore_index=True)\n\n\ndef write_unique_rows(out_dir, sample, set_label, rows_df):\n \"\"\"Write rows unique to a comparison set to a per-sample TSV file.\n\n Output files are named as:\n ``_.unique.tsv``\n\n Filenames are sanitized to ensure filesystem safety. If the\n DataFrame is empty, no file is written.\n\n Args:\n out_dir (str | Path): Output directory for TSV files.\n sample (str): Sample identifier used in the output filename.\n set_label (str): Label indicating which comparison set\n (e.g., SetA or SetB) the rows belong to.\n rows_df (pandas.DataFrame): DataFrame containing rows\n unique to the given set and sample.\n\n Returns:\n None\n \"\"\"\n if rows_df.empty:\n return\n\n out_dir = Path(out_dir)\n out_dir.mkdir(parents=True, exist_ok=True)\n\n fname = f\"{safe_name(sample)}_{safe_name(set_label)}.unique.tsv\"\n out_path = out_dir / fname\n rows_df.to_csv(out_path, sep=\"\\t\", index=False)\n\n\ndef main():\n \"\"\"CLI entry point for comparing two sets of TSV files.\n\n This function:\n 1. Parses command-line arguments.\n 2. Loads TSV files for two comparison sets (A and B).\n 3. Compares breakpoint pairs (LeftBreakpoint|RightBreakpoint)\n per sample.\n 4. Computes counts of:\n - Entries unique to set A\n - Entries shared between A and B\n - Entries unique to set B\n 5. Generates a stacked bar plot (count or percent).\n 6. Optionally writes per-sample TSVs of unique entries.\n\n Returns:\n None\n \"\"\"\n ap = argparse.ArgumentParser(\n description=(\n \"Compare two sets of TSVs (A vs B) by \"\n \"LeftBreakpoint|RightBreakpoint per Sample, \"\n \"plot counts and optionally dump unique entries.\"\n )\n )\n ap.add_argument(\"-a\", \"--file-a\", nargs=\"+\", required=True,\n help=\"One or more TSV files belonging to set A\")\n ap.add_argument(\"-b\", \"--file-b\", nargs=\"+\", required=True,\n help=\"One or more TSV files belonging to set B\")\n\n ap.add_argument(\"--label-a\", default=\"SetA\",\n help=\"Label to use for set A in the legend and filenames\")\n ap.add_argument(\"--label-b\", default=\"SetB\",\n help=\"Label to use for set B in the legend and filenames\")\n\n ap.add_argument(\"--title\",\n help=\"Plot title\")\n ap.add_argument(\"--tool\", default=\"Tool\",\n help=\"Tool in comparison\")\n ap.add_argument(\"--figsize\", nargs=2, type=float, default=[10, 4],\n metavar=(\"W\", \"H\"), help=\"Figure size in inches\")\n\n ap.add_argument(\"-o\", \"--out\", required=True,\n help=\"Output figure path (e.g., overlap.png, overlap.pdf)\")\n\n ap.add_argument(\"--dump-unique-dir\", default=None,\n help=(\n \"Directory where per-sample TSVs unique to A \"\n \"and B will be written (optional)\"\n ))\n\n ap.add_argument(\"--percent\", action=\"store_true\",\n help=\"Stack the bars by percentage rather than raw count\")\n\n args = ap.parse_args()\n\n if not args.title:\n args.title = (\n f\"{args.label_a} vs {args.label_b} Fusion {args.tool} \"\n f\"{'Percent' if args.percent else 'Count'} Comparison\"\n )\n\n A = load_files_with_key(args.file_a)\n B = load_files_with_key(args.file_b)\n\n samples = sorted(\n set(A[SAMPLE_COL].dropna() if SAMPLE_COL in A else []) |\n set(B[SAMPLE_COL].dropna() if SAMPLE_COL in B else [])\n )\n\n col_unique_a = f\"Unique to {args.label_a}\"\n col_unique_b = f\"Unique to {args.label_b}\"\n col_shared = \"Shared\"\n\n rows = []\n for s in samples:\n Arows = A.loc[A[SAMPLE_COL] == s].copy() if not A.empty else pd.DataFrame(columns=A.columns)\n Brows = B.loc[B[SAMPLE_COL] == s].copy() if not B.empty else pd.DataFrame(columns=B.columns)\n\n Akeys = set(Arows[\"key\"].dropna())\n Bkeys = set(Brows[\"key\"].dropna())\n\n uniqueA_keys = Akeys - Bkeys\n uniqueB_keys = Bkeys - Akeys\n shared_keys = Akeys & Bkeys\n\n rows.append({\n \"Sample\": s,\n col_unique_a: len(uniqueA_keys),\n col_shared: len(shared_keys),\n col_unique_b: len(uniqueB_keys),\n })\n\n if args.dump_unique_dir:\n uniqueA_rows = Arows[Arows[\"key\"].isin(uniqueA_keys)].drop(columns=[\"key\"], errors=\"ignore\")\n write_unique_rows(args.dump_unique_dir, s, args.label_a, uniqueA_rows)\n\n uniqueB_rows = Brows[Brows[\"key\"].isin(uniqueB_keys)].drop(columns=[\"key\"], errors=\"ignore\")\n write_unique_rows(args.dump_unique_dir, s, args.label_b, uniqueB_rows)\n\n counts = pd.DataFrame(rows)\n\n pivot = counts.set_index(\"Sample\")[[col_unique_a, col_shared, col_unique_b]]\n\n if args.percent:\n totals = pivot.sum(axis=1)\n pivot = pivot.div(totals.where(totals != 0, 1), axis=0) * 100.0\n ylabel = \"Percent (%)\"\n else:\n ylabel = \"Count (by LeftBreakpoint|RightBreakpoint)\"\n\n samples_order = pivot.index.tolist()\n categories = pivot.columns.tolist()\n\n plt.figure(figsize=args.figsize)\n ax = plt.gca()\n bottom = None\n colors = sns.color_palette(\"tab10\", n_colors=len(categories))\n\n for i, cat in enumerate(categories):\n vals = pivot[cat].values\n if bottom is None:\n ax.bar(samples_order, vals, label=cat, color=colors[i])\n bottom = vals.copy()\n else:\n ax.bar(samples_order, vals, bottom=bottom, label=cat, color=colors[i])\n bottom = bottom + vals\n\n ax.set_xlabel(\"Sample\")\n ax.set_ylabel(ylabel)\n ax.set_title(args.title)\n ax.legend()\n\n if len(samples_order) > 5:\n plt.xticks(rotation=45, ha=\"right\")\n\n plt.tight_layout()\n plt.savefig(args.out, dpi=300, bbox_inches=\"tight\")\n plt.close()\n\n\nif __name__ == \"__main__\":\n main()" ; cwl:entryname "compare_sets.py" ] ] ; cwl:requirements [ a cwl:ShellCommandRequirement ] ; cwl:requirements [ a cwl:InitialWorkDirRequirement ; cwl:listing [ cwl:entry "#!/usr/bin/env python3\nimport argparse\nimport re\nfrom pathlib import Path\n\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\nSAMPLE_COL = \"Sample\"\nLEFT_COL = \"LeftBreakpoint\"\nRIGHT_COL = \"RightBreakpoint\"\n\n\ndef safe_name(s):\n \"\"\"Sanitize a string for safe use in filenames.\n\n Replaces characters that are not alphanumeric, underscore, dot,\n or hyphen with underscores.\n\n Args:\n s (Any): Input value to sanitize (will be converted to string).\n\n Returns:\n str: A filename-safe version of the input string.\n \"\"\"\n return re.sub(r'[^A-Za-z0-9_.-]', '_', str(s))\n\n\ndef load_files_with_key(paths):\n \"\"\"Load multiple TSV files and add a breakpoint-pair key column.\n\n Each input file is read as a tab-separated file with all columns\n treated as strings. A new column named ``key`` is created by\n concatenating ``LeftBreakpoint`` and ``RightBreakpoint`` with\n a pipe character (\"|\").\n\n Args:\n paths (Iterable[str | Path]): Paths to TSV files to load.\n\n Returns:\n pandas.DataFrame: Concatenated DataFrame containing all rows\n from the input files, with an added ``key`` column. If no\n files are provided, an empty DataFrame with ``Sample`` and\n ``key`` columns is returned.\n\n Raises:\n ValueError: If any input file is missing the required\n ``LeftBreakpoint`` or ``RightBreakpoint`` columns.\n \"\"\"\n dfs = []\n for p in paths:\n df = pd.read_csv(p, sep=\"\\t\", dtype=str)\n # ensure breakpoint columns exist\n if LEFT_COL not in df.columns or RIGHT_COL not in df.columns:\n raise ValueError(\n f\"File {p} missing required column(s): {LEFT_COL}, {RIGHT_COL}\"\n )\n df = df.copy()\n df[\"key\"] = df[LEFT_COL].str.cat(df[RIGHT_COL], sep=\"|\")\n dfs.append(df)\n\n if not dfs:\n return pd.DataFrame(columns=[SAMPLE_COL, \"key\"])\n\n return pd.concat(dfs, ignore_index=True)\n\n\ndef write_unique_rows(out_dir, sample, set_label, rows_df):\n \"\"\"Write rows unique to a comparison set to a per-sample TSV file.\n\n Output files are named as:\n ``_.unique.tsv``\n\n Filenames are sanitized to ensure filesystem safety. If the\n DataFrame is empty, no file is written.\n\n Args:\n out_dir (str | Path): Output directory for TSV files.\n sample (str): Sample identifier used in the output filename.\n set_label (str): Label indicating which comparison set\n (e.g., SetA or SetB) the rows belong to.\n rows_df (pandas.DataFrame): DataFrame containing rows\n unique to the given set and sample.\n\n Returns:\n None\n \"\"\"\n if rows_df.empty:\n return\n\n out_dir = Path(out_dir)\n out_dir.mkdir(parents=True, exist_ok=True)\n\n fname = f\"{safe_name(sample)}_{safe_name(set_label)}.unique.tsv\"\n out_path = out_dir / fname\n rows_df.to_csv(out_path, sep=\"\\t\", index=False)\n\n\ndef main():\n \"\"\"CLI entry point for comparing two sets of TSV files.\n\n This function:\n 1. Parses command-line arguments.\n 2. Loads TSV files for two comparison sets (A and B).\n 3. Compares breakpoint pairs (LeftBreakpoint|RightBreakpoint)\n per sample.\n 4. Computes counts of:\n - Entries unique to set A\n - Entries shared between A and B\n - Entries unique to set B\n 5. Generates a stacked bar plot (count or percent).\n 6. Optionally writes per-sample TSVs of unique entries.\n\n Returns:\n None\n \"\"\"\n ap = argparse.ArgumentParser(\n description=(\n \"Compare two sets of TSVs (A vs B) by \"\n \"LeftBreakpoint|RightBreakpoint per Sample, \"\n \"plot counts and optionally dump unique entries.\"\n )\n )\n ap.add_argument(\"-a\", \"--file-a\", nargs=\"+\", required=True,\n help=\"One or more TSV files belonging to set A\")\n ap.add_argument(\"-b\", \"--file-b\", nargs=\"+\", required=True,\n help=\"One or more TSV files belonging to set B\")\n\n ap.add_argument(\"--label-a\", default=\"SetA\",\n help=\"Label to use for set A in the legend and filenames\")\n ap.add_argument(\"--label-b\", default=\"SetB\",\n help=\"Label to use for set B in the legend and filenames\")\n\n ap.add_argument(\"--title\",\n help=\"Plot title\")\n ap.add_argument(\"--tool\", default=\"Tool\",\n help=\"Tool in comparison\")\n ap.add_argument(\"--figsize\", nargs=2, type=float, default=[10, 4],\n metavar=(\"W\", \"H\"), help=\"Figure size in inches\")\n\n ap.add_argument(\"-o\", \"--out\", required=True,\n help=\"Output figure path (e.g., overlap.png, overlap.pdf)\")\n\n ap.add_argument(\"--dump-unique-dir\", default=None,\n help=(\n \"Directory where per-sample TSVs unique to A \"\n \"and B will be written (optional)\"\n ))\n\n ap.add_argument(\"--percent\", action=\"store_true\",\n help=\"Stack the bars by percentage rather than raw count\")\n\n args = ap.parse_args()\n\n if not args.title:\n args.title = (\n f\"{args.label_a} vs {args.label_b} Fusion {args.tool} \"\n f\"{'Percent' if args.percent else 'Count'} Comparison\"\n )\n\n A = load_files_with_key(args.file_a)\n B = load_files_with_key(args.file_b)\n\n samples = sorted(\n set(A[SAMPLE_COL].dropna() if SAMPLE_COL in A else []) |\n set(B[SAMPLE_COL].dropna() if SAMPLE_COL in B else [])\n )\n\n col_unique_a = f\"Unique to {args.label_a}\"\n col_unique_b = f\"Unique to {args.label_b}\"\n col_shared = \"Shared\"\n\n rows = []\n for s in samples:\n Arows = A.loc[A[SAMPLE_COL] == s].copy() if not A.empty else pd.DataFrame(columns=A.columns)\n Brows = B.loc[B[SAMPLE_COL] == s].copy() if not B.empty else pd.DataFrame(columns=B.columns)\n\n Akeys = set(Arows[\"key\"].dropna())\n Bkeys = set(Brows[\"key\"].dropna())\n\n uniqueA_keys = Akeys - Bkeys\n uniqueB_keys = Bkeys - Akeys\n shared_keys = Akeys & Bkeys\n\n rows.append({\n \"Sample\": s,\n col_unique_a: len(uniqueA_keys),\n col_shared: len(shared_keys),\n col_unique_b: len(uniqueB_keys),\n })\n\n if args.dump_unique_dir:\n uniqueA_rows = Arows[Arows[\"key\"].isin(uniqueA_keys)].drop(columns=[\"key\"], errors=\"ignore\")\n write_unique_rows(args.dump_unique_dir, s, args.label_a, uniqueA_rows)\n\n uniqueB_rows = Brows[Brows[\"key\"].isin(uniqueB_keys)].drop(columns=[\"key\"], errors=\"ignore\")\n write_unique_rows(args.dump_unique_dir, s, args.label_b, uniqueB_rows)\n\n counts = pd.DataFrame(rows)\n\n pivot = counts.set_index(\"Sample\")[[col_unique_a, col_shared, col_unique_b]]\n\n if args.percent:\n totals = pivot.sum(axis=1)\n pivot = pivot.div(totals.where(totals != 0, 1), axis=0) * 100.0\n ylabel = \"Percent (%)\"\n else:\n ylabel = \"Count (by LeftBreakpoint|RightBreakpoint)\"\n\n samples_order = pivot.index.tolist()\n categories = pivot.columns.tolist()\n\n plt.figure(figsize=args.figsize)\n ax = plt.gca()\n bottom = None\n colors = sns.color_palette(\"tab10\", n_colors=len(categories))\n\n for i, cat in enumerate(categories):\n vals = pivot[cat].values\n if bottom is None:\n ax.bar(samples_order, vals, label=cat, color=colors[i])\n bottom = vals.copy()\n else:\n ax.bar(samples_order, vals, bottom=bottom, label=cat, color=colors[i])\n bottom = bottom + vals\n\n ax.set_xlabel(\"Sample\")\n ax.set_ylabel(ylabel)\n ax.set_title(args.title)\n ax.legend()\n\n if len(samples_order) > 5:\n plt.xticks(rotation=45, ha=\"right\")\n\n plt.tight_layout()\n plt.savefig(args.out, dpi=300, bbox_inches=\"tight\")\n plt.close()\n\n\nif __name__ == \"__main__\":\n main()" ; cwl:entryname "compare_sets.py" ] ] ; cwl:requirements [ a cwl:InlineJavascriptRequirement ] ; cwl:requirements [ a cwl:DockerRequirement ; DockerRequirement:dockerPull "pgc-images.sbgenomics.com/danmiller/plot_gencode_rnaseq_python:3.13.0" ] ; cwl:requirements [ a cwl:ResourceRequirement ; ResourceRequirement:coresMin "$(inputs.cpu)" ; ResourceRequirement:ramMin "$(inputs.ram * 1000)" ] . rdfs:comment "TSV containing ENSG names on which the TSNE will be performed" ; cwl:inputBinding [ CommandLineBinding:position 2 ; CommandLineBinding:prefix "--gene_list_path" ] ; sld:type cwl:File . cwl:valueFrom "$(false)" .