diff --git a/bin/de_analysis.R b/bin/de_analysis.R index cb31755..83f13b3 100755 --- a/bin/de_analysis.R +++ b/bin/de_analysis.R @@ -1,22 +1,37 @@ #!/usr/bin/env Rscript +suppressMessages(library(argparser)) + +parser <- arg_parser("Run differential expression analysis") +parser <- add_argument(parser, "--annotation", help="Reference annotation.") +parser <- add_argument(parser, "--min_samps_gene_expr", help="Minimum number of samples a gene must be expressed in to be included in differential gene expression.", type="numeric") +parser <- add_argument(parser, "--min_samps_feature_expr", help="Minimum number of samples for differential transcript usage.", type="numeric") +parser <- add_argument(parser, "--min_gene_expr", help="Minimum counts per gene required for differential gene expression.", type="numeric") +parser <- add_argument(parser, "--min_feature_expr", help="Minimum counts per transcript required for differential transcript usage.", type="numeric") +parser <- add_argument(parser, "--sample_sheet", help="Sample sheet.") +parser <- add_argument(parser, "--all_counts", help="All transcript counts CSV file.") +parser <- add_argument(parser, "--de_out_dir", help="Directory where differential expression out files will be saved. Directory will be created if it does not exist", default="de_analysis") +parser <- add_argument(parser, "--merged_out_dir", help="Directory where merged count files will be saved. Directory will be created if it does not exist", default="merged") +argv <- parse_args(parser) + suppressMessages(library("DRIMSeq")) suppressMessages(library("GenomicFeatures")) suppressMessages(library("edgeR")) -args <- commandArgs(trailingOnly=TRUE) -ref_annotation <- args[1] -min_samps_gene_expr <- as.numeric(args[2]) -min_samps_feature_expr <- as.numeric(args[3]) -min_gene_expr <- as.numeric(args[4]) -min_feature_expr <- as.numeric(args[5]) -sample_sheet <- args[6] + +# Create output directories +if (!dir.exists(argv$de_out_dir)){ + dir.create(argv$de_out_dir, recursive=TRUE) +} +if (!dir.exists(argv$merged_out_dir)){ + dir.create(argv$merged_out_dir, recursive=TRUE) +} cat("Loading counts, conditions and parameters.\n") -cts <- as.matrix(read.csv("all_counts.tsv", sep="\t", row.names="Reference", stringsAsFactors=FALSE)) +cts <- as.matrix(read.csv(argv$all_counts, sep="\t", row.names="Reference", stringsAsFactors=FALSE)) # Set up sample data frame: #changed this to sample_id -coldata <- read.csv(sample_sheet, row.names="alias", sep=",", stringsAsFactors=TRUE) +coldata <- read.csv(argv$sample_sheet, row.names="alias", sep=",", stringsAsFactors=TRUE) coldata$sample_id <- rownames(coldata) # check if control condition exists, sets as reference @@ -30,7 +45,7 @@ coldata$condition <- relevel(coldata$condition, ref = "control") # see https://www.ensembl.org/info/website/upload/gff.html # and http://gmod.org/wiki/GFF2#Converting_GFF2_to_GFF3 cat("Checking annotation file type.\n") -lines <- readLines(file(ref_annotation), n=10000) +lines <- readLines(file(argv$annotation), n=10000) # If transcript_id containing '=' (format eg. transcript_id=xxx) # annotation type is gff3 check_file_type <- sum(grepl("transcript_id=", lines)) @@ -49,7 +64,7 @@ if (check_file_type != 0){ # The following handles this. cat("Checking annotation file for presence of transcript_id versions.\n") # Get the first transcript_id from the annotation file by parsing -lines <- readLines(file(ref_annotation), n=100000) +lines <- readLines(file(argv$annotation), n=100000) # Find transcript_ids in first 1000 lines and check if they contain dot (format eg. ENTXXX.1) check_version <- sum(grepl("transcript_id[^;]+\\.", lines)) if (check_version != 0){ @@ -62,7 +77,7 @@ if (check_version != 0){ } cat("Loading annotation database.\n") -txdb <- makeTxDbFromGFF(ref_annotation, format = annotation_type) +txdb <- makeTxDbFromGFF(argv$annotation, format = annotation_type) txdf <- select(txdb, keys(txdb,"GENEID"), "TXNAME", "GENEID") tab <- table(txdf$GENEID) txdf$ntx<- tab[match(txdf$GENEID, names(tab))] @@ -78,15 +93,15 @@ rownames(txdf) <- NULL 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) +write.table(counts, file=file.path(argv$de_out_dir, "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) trs_cts_unfiltered <- counts(d) -d <- dmFilter(d, min_samps_gene_expr = min_samps_gene_expr, min_samps_feature_expr = min_samps_feature_expr, - min_gene_expr = min_gene_expr, min_feature_expr = min_feature_expr) +d <- dmFilter(d, min_samps_gene_expr=argv$min_samps_gene_expr, min_samps_feature_expr=argv$min_samps_feature_expr, + min_gene_expr=argv$min_gene_expr, min_feature_expr=argv$min_feature_expr) cat("Building model matrix.\n") design <- model.matrix(~condition, data=DRIMSeq::samples(d)) @@ -98,12 +113,12 @@ 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/filtered_transcript_counts_with_genes.tsv", sep="\t", row.names = FALSE, quote=FALSE) +write.table(trs_cts, file=file.path(argv$merged_out_dir, "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", quote=FALSE) +write.table(gene_cts, file=file.path(argv$merged_out_dir, "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) @@ -111,7 +126,7 @@ cpm_gene_counts <- cpm(gene_cts) 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) +write.table(cpm_gene_counts, file=file.path(argv$de_out_dir, "cpm_gene_counts.tsv"), sep="\t", quote=FALSE, row.names = FALSE) # Differential gene expression using edgeR: cat("Running differential gene expression analysis using edgeR.\n") @@ -139,7 +154,7 @@ plotMD(qlf, status=status, values=c("up","down","notsig"), hl.col=c("red","blue abline(h=c(-1,1), col="blue") plotQLDisp(fit) -write.table(as.data.frame(edger_res), file="de_analysis/results_dge.tsv", sep="\t") +write.table(as.data.frame(edger_res), file=file.path(argv$de_out_dir, "results_dge.tsv"), sep="\t") # Differential transcript usage using DEXSeq: suppressMessages(library("DEXSeq")) @@ -166,8 +181,8 @@ dxr.g <- dxr.g[order(dxr.g$qval),] dxr_out <- as.data.frame(dxr[,c("featureID", "groupID", "pvalue")]) dxr_out <- dxr_out[order(dxr$pvalue),] -write.table(dxr.g, file="de_analysis/results_dtu_gene.tsv", sep="\t") -write.table(dxr_out, file="de_analysis/results_dtu_transcript.tsv", sep="\t") +write.table(dxr.g, file=file.path(argv$de_out_dir, "results_dtu_gene.tsv"), sep="\t") +write.table(dxr_out, file=file.path(argv$de_out_dir, "results_dtu_transcript.tsv"), sep="\t") # and writing out some of the DEXSeq metrics to accompany EPI2ME Labs tutorial colnames(dxr)[grep("log2fold", colnames(dxr))] <- "log2fold" @@ -175,7 +190,7 @@ MADTUdata <- data.frame(dxr)[order(dxr$padj),c("exonBaseMean", "log2fold", "pval MADTUdata$exonBaseMean <- log2(MADTUdata$exonBaseMean) colnames(MADTUdata)[which(colnames(MADTUdata)=="exonBaseMean")] <- "Log2MeanExon" colnames(MADTUdata)[which(colnames(MADTUdata)=="log2fold")] <- "Log2FC" -write.table(MADTUdata, file="de_analysis/results_dexseq.tsv", sep="\t") +write.table(MADTUdata, file=file.path(argv$de_out_dir, "results_dexseq.tsv"), sep="\t") # stageR analysis of DEXSeq results: cat("stageR analysis\n") @@ -194,4 +209,4 @@ stageRObj <- stageWiseAdjustment(stageRObj, method="dtu", alpha=0.10) suppressWarnings({dex.padj <- getAdjustedPValues(stageRObj, order=FALSE, onlySignificantGenes=FALSE)}) # dex.padj <- dex.padj[,-1] -write.table(dex.padj, file="de_analysis/results_dtu_stageR.tsv", sep="\t") +write.table(dex.padj, file=file.path(argv$de_out_dir, "results_dtu_stageR.tsv"), sep="\t") diff --git a/bin/plot_dtu_results.R b/bin/plot_dtu_results.R index 12850b1..16df20a 100755 --- a/bin/plot_dtu_results.R +++ b/bin/plot_dtu_results.R @@ -1,21 +1,30 @@ #!/usr/bin/env Rscript +suppressMessages(library(argparser)) + +parser <- arg_parser("Plot results") +parser <- add_argument(parser, "--counts", help="Filtered transcript counts with genes.") +parser <- add_argument(parser, "--results_dtu", help="stageR results.") +parser <- add_argument(parser, "--sample_sheet", help="Sample sheet.") +parser <- add_argument(parser, "--pdf_out", help="PDF file name.") +argv <- parse_args(parser) + suppressMessages(library(dplyr)) suppressMessages(library(ggplot2)) suppressMessages(library(tidyr)) # Set up sample data frame: -coldata <- read.csv("sample_sheet.tsv", row.names="alias", sep=",") +coldata <- read.csv(argv$sample_sheet, 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("results_dtu_stageR.tsv", sep="\t") +stageR <- read.csv(argv$results_dtu, sep="\t") names(stageR) <- c("gene_id", "transcript_id", "p_gene", "p_transcript"); # Read filtered counts: -counts <- read.csv("filtered_transcript_counts_with_genes.tsv", sep="\t"); +counts <- read.csv(argv$counts, sep="\t"); names(counts)[2]<-"transcript_id" # Join counts and stageR results: @@ -44,13 +53,13 @@ sig_level <- 0.05 genes <- as.character(tdf[which(tdf$p_gene < sig_level),]$gene_id) genes <- unique(genes) -pdf("dtu_plots.pdf") +pdf(argv$pdf_out) for(gene in genes){ gdf<-tdf[which(tdf$gene_id==gene),] p_gene <- unique(gdf$p_gene) - p <- ggplot(gdf, aes(x=transcript_id, y=norm_count)) + geom_bar(stat="identity", aes(fill=sample), position="dodge") - p <- p + facet_wrap(~ group) + coord_flip() - p <- p + ggtitle(paste(gene," : p_value=",p_gene,sep="")) - print(p) + dtu_plot <- ggplot(gdf, aes(x=transcript_id, y=norm_count)) + geom_bar(stat="identity", aes(fill=sample), position="dodge") + dtu_plot <- dtu_plot + facet_wrap(~ group) + coord_flip() + dtu_plot <- dtu_plot + ggtitle(paste(gene," : p_value=",p_gene,sep="")) + print(dtu_plot) } diff --git a/nextflow.config b/nextflow.config index 6154534..d45d0e5 100644 --- a/nextflow.config +++ b/nextflow.config @@ -94,7 +94,7 @@ params { "--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv'", ] agent = null - container_sha = "shad8671ea3a8ed52f2c0f40355e8eb5c6f00d2cbda" + container_sha = "shac733d952a14257cf3c5c5d5d44c6aed84d5fe5a1" common_sha = "shaabceef445fb63214073cbf5836fdd33c04be4ac7" } } diff --git a/subworkflows/differential_expression.nf b/subworkflows/differential_expression.nf index 29863b4..0483c89 100644 --- a/subworkflows/differential_expression.nf +++ b/subworkflows/differential_expression.nf @@ -81,9 +81,16 @@ process deAnalysis { path "de_analysis/results_dtu.pdf", emit: dtu_pdf path "de_analysis/cpm_gene_counts.tsv", emit: cpm """ - mkdir merged - mkdir de_analysis - de_analysis.R annotation.gtf $params.min_samps_gene_expr $params.min_samps_feature_expr $params.min_gene_expr $params.min_feature_expr "sample_sheet.csv" + de_analysis.R \ + --annotation annotation.gtf \ + --min_samps_gene_expr $params.min_samps_gene_expr \ + --min_samps_feature_expr $params.min_samps_feature_expr \ + --min_gene_expr $params.min_gene_expr \ + --min_feature_expr $params.min_feature_expr \ + --sample_sheet sample_sheet.csv \ + --all_counts all_counts.tsv \ + --de_out_dir de_analysis \ + --merged_out_dir merged """ } @@ -99,7 +106,11 @@ process plotResults { output: path "dtu_plots.pdf", emit: dtu_plots """ - plot_dtu_results.R + plot_dtu_results.R \ + --counts filtered_transcript_counts_with_genes.tsv \ + --results_dtu results_dtu_stageR.tsv \ + --sample_sheet sample_sheet.tsv \ + --pdf_out dtu_plots.pdf """ }