kfukuta/wf-transcriptomes-v202
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Use args for R script inputs [CW-3398]

Browse Source
This commit is contained in:
Kiah McIntosh 2025-02-07 08:42:49 +00:00
parent 91a514abd2
commit eab8ef1efe
4 changed files with 71 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
bin/de_analysis.R
View File

@ -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")

25
bin/plot_dtu_results.R
View File

@ -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)
}

2
nextflow.config
View File

@ -94,7 +94,7 @@ params {
"--sample_sheet 'wf-transcriptomes-demo/sample_sheet.csv'",
]
agent = null
container_sha = "shad8671ea3a8ed52f2c0f40355e8eb5c6f00d2cbda"
container_sha = "shac733d952a14257cf3c5c5d5d44c6aed84d5fe5a1"
common_sha = "shaabceef445fb63214073cbf5836fdd33c04be4ac7"
}
}

19
subworkflows/differential_expression.nf
View File

@ -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
"""
}