# RNA Atlas: export a Seurat RDS on YOUR computer, without uploading it.
# RStudio: open this script and click Source. Choose your original RDS file.
# Command line: Rscript seurat-export.R "path/to/object.rds"
# Requires an installed SeuratObject package compatible with your object.
# No filtering, normalization, significance testing, or sample labels are inferred.
# Seurat v5 layer API: https://satijalab.org/seurat/articles/seurat5_essential_commands.html

write_rna_atlas_export <- function(x, metadata, output_dir, umap = NULL) {
  dimensions <- attr(x, "Dim")
  identifiers <- attr(x, "Dimnames")
  p <- attr(x, "p"); i <- attr(x, "i"); v <- attr(x, "x")
  if (length(dimensions) != 2L || is.null(p) || is.null(i) || is.null(v))
    stop("Export needs an in-memory dgCMatrix raw-count layer. Use a native workflow for external BPCells layers.")
  genes <- identifiers[[1]]; cells <- identifiers[[2]]
  stopifnot(length(genes) == dimensions[1], length(cells) == dimensions[2],
            !anyNA(genes), !anyNA(cells), all(nzchar(genes)), all(nzchar(cells)),
            !anyDuplicated(genes), !anyDuplicated(cells),
            length(p) == length(cells) + 1L, p[1] == 0, tail(p, 1) == length(v),
            length(i) == length(v), all(diff(p) >= 0),
            all(i >= 0 & i < length(genes)), all(is.finite(v)),
            all(v >= 0), all(v == floor(v)), all(v <= 1e12))
  if (!all(cells %in% rownames(metadata))) stop("Cell metadata IDs do not match raw-count columns.")
  metadata <- metadata[cells, , drop = FALSE]
  if (any(c("cell", "UMAP_1", "UMAP_2") %in% names(metadata)))
    names(metadata) <- make.unique(ifelse(names(metadata) %in% c("cell", "UMAP_1", "UMAP_2"), paste0("original_", names(metadata)), names(metadata)))
  if (!is.null(umap) && ncol(umap) >= 2L && all(cells %in% rownames(umap))) {
    metadata$UMAP_1 <- umap[cells, 1]
    metadata$UMAP_2 <- umap[cells, 2]
  }
  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  if (any(file.exists(file.path(output_dir, c("matrix.mtx.gz", "features.tsv.gz", "barcodes.tsv.gz", "cell_metadata.csv")))))
    stop("Output files already exist. Choose a new output directory; nothing was overwritten.")
  con <- gzfile(file.path(output_dir, "matrix.mtx.gz"), "wt")
  tryCatch({
    writeLines(c("%%MatrixMarket matrix coordinate integer general",
                 "% RNA Atlas: unchanged raw counts exported locally from Seurat",
                 sprintf("%d %d %.0f", length(genes), length(cells), sum(v != 0))), con)
    for (column in seq_along(cells)) {
      first <- p[column] + 1; last <- p[column + 1]
      if (first <= last) for (start in seq.int(first, last, by = 50000)) {
        at <- seq.int(start, min(last, start + 49999)); at <- at[v[at] != 0]
        if (length(at)) writeLines(sprintf("%d %d %.0f", i[at] + 1L, column, v[at]), con)
      }
    }
  }, finally = close(con))
  write_gz_table <- function(table, filename) {
    connection <- gzfile(file.path(output_dir, filename), "wt")
    tryCatch(write.table(table, connection, sep = "\t", quote = FALSE,
                         row.names = FALSE, col.names = FALSE), finally = close(connection))
  }
  write_gz_table(data.frame(genes, genes, "Gene Expression"), "features.tsv.gz")
  write_gz_table(data.frame(cells), "barcodes.tsv.gz")
  write.csv(data.frame(cell = cells, metadata, check.names = FALSE),
            file.path(output_dir, "cell_metadata.csv"), row.names = FALSE)
  writeLines(c("Raw counts exported without filtering or normalization.",
               paste(length(genes), "genes;", length(cells), "cells;", sum(v != 0), "stored nonzero entries."),
               "Cluster labels/other metadata and matching supplied UMAP coordinates are retained.",
               "Sample IDs must represent real biological samples; technical libraries and pooled mice are not independent replicates.",
               "Open RNA Atlas > Import data > 10x sparse matrix. Select the three .gz files, inspect, then load.",
               "Next choose Add metadata CSV / TSV and select cell_metadata.csv.",
               "For the original UMAP, select PCA & UMAP > Display > Supplied embedding.",
               "A count matrix may still exceed browser analysis memory; the import preflight will check it."),
             file.path(output_dir, "README.txt"))
  writeLines(capture.output(sessionInfo()), file.path(output_dir, "R_sessionInfo.txt"))
  invisible(normalizePath(output_dir))
}

export_seurat_for_rna_atlas <- function(input = NULL, assay_name = "RNA", output_dir = NULL) {
  if (!requireNamespace("SeuratObject", quietly = TRUE))
    stop("Install SeuratObject in your local R environment first: install.packages('SeuratObject').")
  if (is.null(input)) {
    if (!interactive()) stop("Provide an RDS path: Rscript seurat-export.R 'path/to/object.rds'")
    input <- file.choose()
  }
  object <- readRDS(input)
  if (!inherits(object, "Seurat")) stop("The RDS is not a Seurat object.")
  if (!assay_name %in% names(object@assays)) stop("Raw RNA assay not found. Available assays: ", paste(names(object@assays), collapse = ", "))
  assay <- object[[assay_name]]
  if (inherits(assay, "Assay5")) {
    count_layers <- SeuratObject::Layers(assay, search = "^counts($|\\.)")
    if (!length(count_layers)) stop("No raw-count layer found. Normalized/integrated values are not raw counts.")
    if (length(count_layers) > 1L) {
      message("Joining split raw-count layers locally: ", paste(count_layers, collapse = ", "))
      assay <- SeuratObject::JoinLayers(assay, layers = "counts", new = "counts")
      count_layers <- "counts"
    }
    x <- SeuratObject::LayerData(assay, layer = count_layers[[1]])
  } else if (utils::packageVersion("SeuratObject") >= "5.0.0") {
    x <- SeuratObject::GetAssayData(object, assay = assay_name, layer = "counts")
  } else x <- SeuratObject::GetAssayData(object, assay = assay_name, slot = "counts")
  umap <- if ("umap" %in% names(object@reductions)) SeuratObject::Embeddings(object[["umap"]]) else NULL
  if (is.null(output_dir)) output_dir <- file.path(dirname(input), paste0(tools::file_path_sans_ext(basename(input)), "_RNAAtlas_", format(Sys.time(), "%Y%m%d_%H%M%S")))
  result <- write_rna_atlas_export(x, object[[]], output_dir, umap)
  message("Done. Original RDS unchanged. Export folder: ", result)
  invisible(result)
}

if (interactive() || sys.nframe() == 0L) {
  arguments <- commandArgs(trailingOnly = TRUE)
  export_seurat_for_rna_atlas(if (length(arguments)) arguments[[1]] else NULL)
}