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

Fix and simplify column name checking

Browse Source
This commit is contained in:
Neil Horner 2025-08-01 09:14:20 +00:00
parent 7a94b2c515
commit 4fe583d4ea
2 changed files with 18 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
bin/de_analysis.R
View File

@ -21,8 +21,8 @@ suppressMessages(library("edgeR"))
# Some functions, including dmDSdata, converts '.' in sample IDs to '-'. # Some functions, including dmDSdata, converts '.' in sample IDs to '-'.
# Make the output DF match the sample IDs in the sample sheet # Make the output DF match the sample IDs in the sample sheet
# start_col is the index of the first sample column in the data frame # start_col is the index of the first sample column in the data frame
rename_sample_columns <- function(df, sample_ids, start_col) { rename_sample_columns <- function(df, alias, start_col) {
colnames(df)[start_col:ncol(df)] <- sample_ids colnames(df)[start_col:ncol(df)] <- alias
return(df) return(df)
} }
@ -37,13 +37,14 @@ if (!dir.exists(argv$merged_out_dir)){
cat("Loading counts, conditions and parameters.\n") cat("Loading counts, conditions and parameters.\n")
cts <- as.matrix(read.csv(argv$all_counts, sep="\t", row.names="Reference", stringsAsFactors=FALSE, check.names=FALSE)) cts <- as.matrix(read.csv(argv$all_counts, sep="\t", row.names="Reference", stringsAsFactors=FALSE, check.names=FALSE))
# Set up sample data frame: # Set up sample data frame:
#changed this to sample_id
coldata <- read.csv(argv$sample_sheet, row.names="alias", sep=",", stringsAsFactors=TRUE, check.names=FALSE) coldata <- read.csv(argv$sample_sheet, row.names="alias", sep=",", stringsAsFactors=TRUE, check.names=FALSE)
coldata$alias <- rownames(coldata)
#dmDSdata looks up sample_id in coldata, we need this to be sample alias
coldata$sample_id <- coldata$alias
coldata$sample_id <- rownames(coldata)
# Reorder the input counts columns to match the sample sheet order. # Reorder the input counts columns to match the sample sheet order.
# This ensures we don't missassign column names, when renaming output DF columns. # This ensures we don't missassign column names, when renaming output DF columns.
cts <- cts[, coldata$sample_id, drop = FALSE] cts <- cts[, coldata$alias, drop = FALSE]
# check if control condition exists, sets as reference # check if control condition exists, sets as reference
if(!"control" %in% coldata$condition) if(!"control" %in% coldata$condition)
stop("sample_sheet.csv does not contain 'control' stop("sample_sheet.csv does not contain 'control'
@ -101,10 +102,11 @@ rownames(txdf) <- NULL
# Create counts data frame: # Create counts data frame:
counts<-data.frame(gene_id=txdf$GENEID, feature_id=txdf$TXNAME, cts) counts<-data.frame(gene_id=txdf$GENEID, feature_id=txdf$TXNAME, cts)
counts <- rename_sample_columns(counts, coldata$sample_id, 3) counts <- rename_sample_columns(counts, coldata$alias, 3)
# output unfiltered version of the counts table now we have paired transcripts with gene ids # output unfiltered version of the counts table now we have paired transcripts with gene ids
write.table(counts, file=file.path(argv$de_out_dir, "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") cat("Filtering counts using DRIMSeq.\n")
d <- dmDSdata(counts=counts, samples=coldata) d <- dmDSdata(counts=counts, samples=coldata)
@ -122,13 +124,13 @@ suppressMessages(library("dplyr"))
# Sum transcript counts into gene counts: # Sum transcript counts into gene counts:
cat("Sum transcript counts into gene counts.\n") cat("Sum transcript counts into gene counts.\n")
trs_cts <- counts(d) trs_cts <- counts(d)
trs_cts <- rename_sample_columns(trs_cts, coldata$sample_id, 3) trs_cts <- rename_sample_columns(trs_cts, coldata$alias, 3)
write.table(trs_cts, file=file.path(argv$merged_out_dir, "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() 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 rownames(gene_cts) <- gene_cts$gene_id
gene_cts$gene_id <- NULL gene_cts$gene_id <- NULL
gene_cts <- rename_sample_columns(gene_cts, coldata$sample_id, 1) gene_cts <- rename_sample_columns(gene_cts, coldata$alias, 1)
write.table(gene_cts, file=file.path(argv$merged_out_dir, "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 # Output count per million of the gene counts using edgeR CPM
@ -137,7 +139,7 @@ cpm_gene_counts <- cpm(gene_cts)
cpm_gene_counts <- cbind(var_name = rownames(cpm_gene_counts), cpm_gene_counts) cpm_gene_counts <- cbind(var_name = rownames(cpm_gene_counts), cpm_gene_counts)
rownames(cpm_gene_counts) <- NULL rownames(cpm_gene_counts) <- NULL
colnames(cpm_gene_counts)[1] <- "gene_id" colnames(cpm_gene_counts)[1] <- "gene_id"
cpm_gene_counts <- rename_sample_columns(cpm_gene_counts, coldata$sample_id, 2) cpm_gene_counts <- rename_sample_columns(cpm_gene_counts, coldata$alias, 2)
write.table(cpm_gene_counts, file=file.path(argv$de_out_dir, "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: # Differential gene expression using edgeR:

11
subworkflows/differential_expression.nf
View File

@ -78,6 +78,7 @@ process deAnalysis {
path "de_analysis/results_dtu_stageR.tsv", emit: dtu_stageR path "de_analysis/results_dtu_stageR.tsv", emit: dtu_stageR
path "de_analysis/results_dtu.pdf", emit: dtu_pdf path "de_analysis/results_dtu.pdf", emit: dtu_pdf
path "de_analysis/cpm_gene_counts.tsv", emit: cpm path "de_analysis/cpm_gene_counts.tsv", emit: cpm
script:
""" """
de_analysis.R \ de_analysis.R \
--annotation annotation.gtf \ --annotation annotation.gtf \
@ -91,11 +92,13 @@ process deAnalysis {
--merged_out_dir merged --merged_out_dir merged
# Check that the original aliases in the input TSV have not been mangled by R's read.csv or other functions # Check that the original aliases in the input TSV have not been mangled by R's read.csv or other functions
head -1 all_counts.tsv | cut -f2- | tr '\t' '\n' > expected_colnames # Sample column order should be in the same order as the input sample sheet
alias_col=\$(awk -v RS=',' '/alias/{print NR; exit}' "sample_sheet.csv")
cut -d',' -f3 sample_sheet.csv | tail -n +2 | paste -sd '\t' - | sed 's/\t*\$//' > expected_colnames
head -1 de_analysis/cpm_gene_counts.tsv | cut -f2- | tr '\t' '\n' > cpm_gene_counts_colnames head -1 de_analysis/cpm_gene_counts.tsv | cut -f2- > cpm_gene_counts_colnames
head -1 merged/all_gene_counts.tsv | tr '\t' '\n' > merged_counts_colnames head -1 merged/all_gene_counts.tsv > merged_counts_colnames
head -1 merged/filtered_transcript_counts_with_genes.tsv | cut -f3- | tr '\t' '\n' > merged_filtered_colnames head -1 merged/filtered_transcript_counts_with_genes.tsv | cut -f3- > merged_filtered_colnames
# Check for mismatches in sample column names # Check for mismatches in sample column names
for file in cpm_gene_counts_colnames merged_counts_colnames merged_filtered_colnames; do for file in cpm_gene_counts_colnames merged_counts_colnames merged_filtered_colnames; do