Processing Data_Ute_RNAseq_FINAL/Data_RNA-Seq_MKL-1+WaGa

Leave a reply

Here is the optimized version. Key improvements: (1) it verifies that the two Geneid/gene_name columns from the pasted files are identical and keeps only one copy; (2) it never blindly overwrites column names — R’s automatic header sanitization (check.names=TRUE) produces exactly the names you were assigning by hand (e.g. 042_MKL-1_wt_EV → X042_MKL.1_wt_EV, MKL-1_EV-RNA_87 → MKL.1_EV.RNA_87), so we use those; (3) the final friendly rename is done via an explicit name→name map instead of position-based colnames(x) <- c(...), so a sample can never be silently mislabeled — any mismatch throws an error.

# ---------------------------------------------------------------
# 0) Read the merged count table
#    - Do NOT set row.names=1 yet: we first need to compare the two
#      Geneid / gene_name columns that come from the two pasted files.
#    - check.names=TRUE (default) sanitizes the header deterministically:
#        042_MKL-1_wt_EV        -> X042_MKL.1_wt_EV
#        MKL-1_EV-RNA_87        -> MKL.1_EV.RNA_87
#        2nd Geneid / gene_name -> Geneid.1 / gene_name.1
#    (i.e. exactly the names you previously assigned by hand)
# ---------------------------------------------------------------
d.full <- read.delim2("merged_gene_counts_40samples.txt", sep="\t", header=TRUE)

stopifnot(ncol(d.full) == 44)   # 2 x (Geneid + gene_name) + 40 samples
colnames(d.full)                # eyeball-check the auto-generated names

# ---------------------------------------------------------------
# 1) Check that the two Geneid / gene_name columns are identical
#    If this fails, the two files were pasted in different gene order
#    -> stop and fix (diagnosis: which(d.full$Geneid != d.full$Geneid.1))
# ---------------------------------------------------------------
stopifnot(identical(d.full$Geneid,    d.full$Geneid.1))
stopifnot(identical(d.full$gene_name, d.full$gene_name.1))

# identical -> keep only one copy of each
d.full$Geneid.1    <- NULL
d.full$gene_name.1 <- NULL

# Ensembl IDs as row names; gene_name was also dropped before DESeq2 previously
stopifnot(!anyDuplicated(d.full$Geneid))
rownames(d.full) <- d.full$Geneid
d.full$Geneid    <- NULL
d.full$gene_name <- NULL

stopifnot(ncol(d.full) == 40)   # only sample count columns remain

# ---------------------------------------------------------------
# 2) Reorder columns by biological group, using the AUTO-GENERATED names
#    (order matches the condition/donor/batch vectors of the old script)
# ---------------------------------------------------------------
col_order <- c(
  # MKL-1 RNA
  "MKL.1_RNA","MKL.1_RNA_118","MKL.1_RNA_147",
  # MKL-1 wt EV
  "MKL.1_EV.RNA","MKL.1_EV.RNA_2","MKL.1_EV.RNA_118","MKL.1_EV.RNA_87","MKL.1_EV.RNA_27",
  "X042_MKL.1_wt_EV",
  # MKL-1 EV DMSO / Dox
  "X042_MKL.1_sT_DMSO","X0505_MKL.1_sT_DMSO_EV",
  "X042_MKL.1_scr_DMSO_EV","X0505_MKL.1_scr_DMSO_EV",
  "X042_MKL.1_sT_Dox","X0505_MKL.1_sT_Dox_EV",
  "X042_MKL.1_scr_Dox_EV","X0505_MKL.1_scr_Dox_EV",
  # WaGa RNA
  "WaGa_RNA","WaGa_RNA_118","WaGa_RNA_147",
  # WaGa wt EV
  "WaGa_EV.RNA","WaGa_EV.RNA_2","WaGa_EV.RNA_118","WaGa_EV.RNA_147","WaGa_EV.RNA_226",
  "X1107_WaGa_wt_EV","X1605_WaGa_wt_EV","X2706_WaGa_wt_EV",
  # WaGa EV DMSO / Dox
  "X1107_WaGa_sT_DMSO_EV","X1605_WaGa_sT_DMSO_EV","X2706_WaGa_sT_DMSO_EV",
  "X1107_WaGa_scr_DMSO_EV","X1605_WaGa_scr_DMSO_EV","X2706_WaGa_scr_DMSO_EV",
  "X1107_WaGa_sT_Dox_EV","X1605_WaGa_sT_Dox_EV","X2706_WaGa_sT_Dox_EV",
  "X1107_WaGa_scr_Dox_EV","X1605_WaGa_scr_Dox_EV","X2706_WaGa_scr_Dox_EV")

stopifnot(length(col_order) == 40,
          !anyDuplicated(col_order),
          all(col_order %in% colnames(d.full)))   # errors instead of mislabeling

reordered.raw <- d.full[, col_order]

# ---------------------------------------------------------------
# 3) Friendly sample names via an EXPLICIT map (position-independent).
#    Kept because all downstream code refers to e.g. "MKL-1 EV sT DMSO 042".
# ---------------------------------------------------------------
name_map <- c(
  "MKL.1_RNA"               = "MKL-1 RNA",
  "MKL.1_RNA_118"           = "MKL-1 RNA 118",
  "MKL.1_RNA_147"           = "MKL-1 RNA 147",
  "MKL.1_EV.RNA"            = "MKL-1 EV",
  "MKL.1_EV.RNA_2"          = "MKL-1 EV 2",
  "MKL.1_EV.RNA_118"        = "MKL-1 EV 118",
  "MKL.1_EV.RNA_87"         = "MKL-1 EV 87",
  "MKL.1_EV.RNA_27"         = "MKL-1 EV 27",
  "X042_MKL.1_wt_EV"        = "MKL-1 EV 042",
  "X042_MKL.1_sT_DMSO"      = "MKL-1 EV sT DMSO 042",
  "X0505_MKL.1_sT_DMSO_EV"  = "MKL-1 EV sT DMSO 0505",
  "X042_MKL.1_scr_DMSO_EV"  = "MKL-1 EV scr DMSO 042",
  "X0505_MKL.1_scr_DMSO_EV" = "MKL-1 EV scr DMSO 0505",
  "X042_MKL.1_sT_Dox"       = "MKL-1 EV sT Dox 042",
  "X0505_MKL.1_sT_Dox_EV"   = "MKL-1 EV sT Dox 0505",
  "X042_MKL.1_scr_Dox_EV"   = "MKL-1 EV scr Dox 042",
  "X0505_MKL.1_scr_Dox_EV"  = "MKL-1 EV scr Dox 0505",
  "WaGa_RNA"                = "WaGa RNA",
  "WaGa_RNA_118"            = "WaGa RNA 118",
  "WaGa_RNA_147"            = "WaGa RNA 147",
  "WaGa_EV.RNA"             = "WaGa EV",
  "WaGa_EV.RNA_2"           = "WaGa EV 2",
  "WaGa_EV.RNA_118"         = "WaGa EV 118",
  "WaGa_EV.RNA_147"         = "WaGa EV 147",
  "WaGa_EV.RNA_226"         = "WaGa EV 226",
  "X1107_WaGa_wt_EV"        = "WaGa EV 1107",
  "X1605_WaGa_wt_EV"        = "WaGa EV 1605",
  "X2706_WaGa_wt_EV"        = "WaGa EV 2706",
  "X1107_WaGa_sT_DMSO_EV"   = "WaGa EV sT DMSO 1107",
  "X1605_WaGa_sT_DMSO_EV"   = "WaGa EV sT DMSO 1605",
  "X2706_WaGa_sT_DMSO_EV"   = "WaGa EV sT DMSO 2706",
  "X1107_WaGa_scr_DMSO_EV"  = "WaGa EV scr DMSO 1107",
  "X1605_WaGa_scr_DMSO_EV"  = "WaGa EV scr DMSO 1605",
  "X2706_WaGa_scr_DMSO_EV"  = "WaGa EV scr DMSO 2706",
  "X1107_WaGa_sT_Dox_EV"    = "WaGa EV sT Dox 1107",
  "X1605_WaGa_sT_Dox_EV"    = "WaGa EV sT Dox 1605",
  "X2706_WaGa_sT_Dox_EV"    = "WaGa EV sT Dox 2706",
  "X1107_WaGa_scr_Dox_EV"   = "WaGa EV scr Dox 1107",
  "X1605_WaGa_scr_Dox_EV"   = "WaGa EV scr Dox 1605",
  "X2706_WaGa_scr_Dox_EV"   = "WaGa EV scr Dox 2706")

stopifnot(all(colnames(reordered.raw) %in% names(name_map)),
          !anyDuplicated(unname(name_map[colnames(reordered.raw)])))
colnames(reordered.raw) <- unname(name_map[colnames(reordered.raw)])

# ---------------------------------------------------------------
# 4) Write out and filter (same as before)
# ---------------------------------------------------------------
write.csv(reordered.raw, file="counts.txt")

# IMPORTANT: filter low-count genes at this step!
d <- reordered.raw[rowSums(reordered.raw > 3) > 2, ]

What changed and why it is safer

Old code New code Why
row.names=1 at import Import first, compare, then set row names Allows checking the two Geneid/gene_name columns before one is discarded
No check of duplicated id columns identical() + stopifnot Catches a mis-pasted file (different gene order) immediately
Manual colnames(d.raw) <- c(...43 names...) R’s automatic sanitization (which yields the same names) No forced overwrite; names always come from the file header
Position-based colnames(reordered.raw) <- c(...) Explicit name_map keyed by actual column name A sample can never be mislabeled; any unknown/duplicate name throws an error
Filename "merged_gene_counts_40samples.txt " (trailing space) Corrected The stray space would make read.delim2 fail

Note: if stopifnot(identical(...)) fails, the two pasted files do not have the same gene order — in that case do not proceed; re-merge them properly (e.g. with a join on Geneid rather than paste/cbind).

Yes — I checked all six metadata vectors against the reordered column order (the 17 MKL-1 + 23 WaGa samples from your col_order/rename step). All vectors have exactly 40 entries, and every entry maps to the correct sample. The code is correct.

Verification summary

Vector MKL-1 block (17) WaGa block (23) Match?
condition_for_pca RNA×3, EV×6, sT.DMSO×2, scr.DMSO×2, sT.Dox×2, scr.Dox×2 RNA×3, EV×8, sT.DMSO×3, scr.DMSO×3, sT.Dox×3, scr.Dox×3 ✓
condition MKL1.RNA×3, MKL1.EV×6, then 2+2+2+2 WaGa.RNA×3, WaGa.EV×8, then 3+3+3+3 ✓
donor 1,118,147 / 1,2,118,87,27,042 / 042,0505 pairs 1,118,147 / 1,2,118,147,226,1107,1605,2706 / 1107,1605,2706 ×4 ✓
batch 2021.08, 2021.09×2 / 2021.08×2, 2021.09×3 / 2022.08×9 2021.08, 2021.09×2 / 2021.08×2, 2021.09×3 / 2022.11×15 ✓
cell.line MKL-1 ×17 WaGa ×23 ✓
ids unique per sample, order matches unique per sample, order matches ✓

Spot checks that confirm the alignment:

  • Position 9: MKL-1 EV 042 → MKL1.EV, donor 042, batch 2022.08 ✓
  • Positions 26–28: WaGa EV 1107/1605/2706 → WaGa.EV, donors 1107/1605/2706, batch 2022.11 ✓
  • The skipped sample 0505_MKL-1_wt_EV is correctly absent (only one MKL-1 wt EV: 042) ✓

Run these assertions to be safe

stopifnot(ncol(d) == 40,
          length(condition) == 40, length(condition_for_pca) == 40,
          length(donor) == 40, length(batch) == 40,
          length(cell.line) == 40, length(ids) == 40)

# metadata consistent with sample names
stopifnot(all(grepl("^MKL-1", colnames(d)[cell.line == "MKL-1"])),
          all(grepl("^WaGa",  colnames(d)[cell.line == "WaGa"])))
stopifnot(all(paste0(ifelse(cell.line == "MKL-1", "MKL1", "WaGa"), ".") ==
              substr(as.character(condition), 1, 5)))

# visual check — eyeball once before DESeq2
data.frame(sample = colnames(d), condition, donor, batch, cell.line)

Three notes (not errors, but worth knowing)

  1. Design ~batch+condition is full rank here — good. The DMSO/Dox conditions exist only in batch 2022.08 (MKL-1) / 2022.11 (WaGa), but the wt-EV bridge samples (MKL-1 EV 042, WaGa EV 1107/1605/2706) connect those batches to MKL1.EV/WaGa.EV in the earlier batches, so all contrasts remain estimable. DESeq2 will stop with an error if this weren’t the case.
  2. donor is recorded but not used in the design. The sT/scr × DMSO/Dox samples are paired by donor (042/0505 and 1107/1605/2706), so a paired design (~batch+donor+condition) could add power — but given your note that you’re unsure whether “rep1” is the same donor across the 2021 and 2022 experiments, keeping ~batch+condition (as in the original analysis) is the safer, consistent choice.
  3. condition_for_pca is currently unused (it’s only referenced in a commented-out design line). Harmless, but you can drop it if you want to tidy up.

Everything downstream (vst(dds), estimateSizeFactors, the results(dds, name=...) contrast names like MKL1.sT.DMSO_vs_MKL1.scr.DMSO) will work exactly as in your established workflow.

Leave a Reply

Your email address will not be published. Required fields are marked *