KJ optionally use h5seurat instead of rds

R/convert_rds_to_anndata.R

@@ -0,0 +1,29 @@
+## Export data to be used as input for analysis with python
+## packages such as scanpy
+
+# Libraries ----
+
+stopifnot(
+  require(optparse),
+  require(Seurat),
+  require(SeuratDisk),
+  require(tenxutils)
+)
+
+# Options ----
+
+option_list <- list(
+  make_option(c("--seuratobject"), default="begin.Robj",
+              help="A seurat object after dimensionality reduction")
+)
+
+opt <- parse_args(OptionParser(option_list=option_list))
+
+cat("Running with options:\n")
+print(opt)
+
+output_file = gsub(".h5seurat", ".h5ad", opt$seuratobject)
+
+Convert(opt$seuratobject, dest = output_file, overwrite = TRUE, verbose = TRUE) 
+
+print("Completed")

R/export_for_python.R

@@ -53,10 +53,10 @@ exportMetaData <- function(seurat_object, outdir=NULL) {
                 quote=FALSE, sep="\t", row.names = FALSE, col.names = TRUE)
     }
 
-
 # Read RDS seurat object
-message("readRDS")
+message(sprintf("readRDS: %s", opt$seuratobject))
 s <- readRDS(opt$seuratobject)
+
 message("export_for_python running with default assay: ", DefaultAssay(s))
 
 # Write out the cell and feature names

R/plot_group_numbers.R

@@ -179,8 +179,8 @@ if(opt$stat %in% c("total_UMI", "ngenes"))
     {
         require(Seurat)
         # add the statistic from the seurat object
-        s <- readRDS(opt$seuratobject)
-
+        s <- loadSeurat(path=opt$seuratobject)
+      
         ## set the default assay
         message("Setting default assay to: ", opt$seuratassay)
         DefaultAssay(s) <- opt$seuratassay

R/plot_rdims_cluster_genes.R

@@ -117,7 +117,7 @@ print(opt)
 ## read in the assay data
 #data <- read.table(opt$assaydata, sep="\t", header=T, as.is=T, check.names=F)
 message("reading in assay data")
-data <- readRDS(opt$assaydata)
+data <- loadSeurat(opt$assaydata)
 
 message("reading in the coordinates")
 plot_data <- read.table(opt$table,sep="\t",header=TRUE)

R/plot_rdims_gene.R

@@ -119,9 +119,7 @@ if(opt$pointpch != ".") { opt$pointpch <- as.numeric(opt$pointpch) }
 cat("Running with options:\n")
 print(opt)
 
-
-## read in the raw count matrix
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 
 ## set the default assay
 message("Setting default assay to: ", opt$seuratassay)

R/plot_violins.R

@@ -59,7 +59,8 @@ genes <- read.table(opt$genetable,
 
 
 ## read in the raw count matrix
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
+
 Idents(s) <- readRDS(opt$clusterids)
 
 ## set the default assay

R/scanpy_post_process_clusters.R

@@ -34,8 +34,7 @@ opt <- parse_args(OptionParser(option_list=option_list))
 cat("Running with options:\n")
 print(opt)
 
-message(sprintf("readRDS: %s", opt$seuratobject))
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 
 message("seurat_cluster.R running with default assay: ", DefaultAssay(s))
 
@@ -71,16 +70,16 @@ names(cluster_ids) <- clusters$barcode
 unique_cluster_ids <- unique(cluster_ids)
 print(unique_cluster_ids)
 
-message(sprintf("saving a unique list ofe cluster ids"))
+print(sprintf("saving a unique list of cluster ids"))
 write.table(unique_cluster_ids, file=file.path(opt$outdir,"cluster_ids.tsv"),
             quote=FALSE, col.names = FALSE, row.names = FALSE)
 
 cluster_colors <- gg_color_hue(length(unique_cluster_ids))
-message(sprintf("saving the cluster colors (ggplot)"))
+print(sprintf("saving the cluster colors (ggplot)"))
 write.table(cluster_colors, file=file.path(opt$outdir,"cluster_colors.tsv"),
             quote=FALSE, col.names = FALSE, row.names = FALSE)
 
-message(sprintf("saveRDS"))
+print(sprintf("saveRDS of cluster ids"))
 saveRDS(cluster_ids, file=file.path(opt$outdir,"cluster_ids.rds"))
 
 cluster_assignments <- data.frame(barcode=as.character(names(cluster_ids)),
@@ -91,6 +90,6 @@ write.table(cluster_assignments,
             sep="\t", col.names=T, row.names=F, quote=F)
 
 
-message("seurat_cluster.R final default assay: ", DefaultAssay(s))
+print(sprintf("seurat_cluster.R final default assay: %s", DefaultAssay(s)))
 
-message("Completed")
+print("Completed")

R/seurat_FindMarkers.R

@@ -78,27 +78,43 @@ plan("multiprocess",
      workers = opt$numcores)
 
 
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 cluster_ids <- readRDS(opt$clusterids)
 
+if (class(s@misc) == "list"){
+  cat("Convert the s@misc to dataframe")
+  s@misc <- as.data.frame(s@misc, stringsAsFactors=FALSE)
+}
+
 have_gene_ids <- "gene_id" %in% colnames(s@misc)
 
 if(!have_gene_ids)
 {
-    cat("Missing ensembl gene_id column in @misc, reading in annotation.\n")
-
+  cat("Missing ensembl gene_id column in @misc, reading in annotation from other source.\n")
+
+  if(endsWith(opt$seuratobject, ".h5seurat")){
+    message("Using s@assays$RNA@meta.features in Seurat object for gene annotation. This slot contains unique gene names and is produced during conversion from scanpy anndata object.")
+    ann <- s@assays$RNA@meta.features
+    ann$seurat_id = rownames(ann)
+    ann <- unique(ann[,c("gene_ids", "seurat_id")])
+    colnames(ann) <- c("ensembl_id", "gene_name")
+  } else {
+    message("Reading in annotation from annotation.dir")
+    message("Please note that only genes with unique gene names will be kept.")
+    message("The other ones will be lost due to Seurat's process of making gene names unique!")
     ann <- read.table(gzfile(opt$annotation), header=T, sep="\t", as.is=T)
 
     ## subset to columns of interest
     ann <- unique(ann[,c("ensembl_id", "gene_name")])
+  }
 
-    ## the gene names have been made unique so we want only 1:1 mappings.
-    gn_tab <- table(ann$gene_name)
+  ## the gene names have been made unique so we want only 1:1 mappings.
+  gn_tab <- table(ann$gene_name)
 
-    unique_gn <- names(gn_tab)[gn_tab == 1]
-    ann <- ann[ann$gene_name %in% unique_gn,]
+  unique_gn <- names(gn_tab)[gn_tab == 1]
+  ann <- ann[ann$gene_name %in% unique_gn,]
 
-    rownames(ann) <- ann$gene_name
+  rownames(ann) <- ann$gene_name
 }
 
 xx <- names(cluster_ids)

R/seurat_begin.R

@@ -42,7 +42,6 @@ option_list <- list(
         c("--matrixtype"), default="10X",
         help="either 10X or rds"
     ),
-
     make_option(
         c("--project"),
         default="SeuratAnalysis",
@@ -53,6 +52,11 @@ option_list <- list(
         default="seurat.out.dir",
         help="Location for outputs files. Must exist."
         ),
+    make_option(
+      c("--file_type"),
+      default="rds",
+      help="Which type of file to use (rds or h5seurat)"
+    ),
     make_option(
         c("--metadata"),
         default="none",
@@ -674,6 +678,11 @@ print(
 message("seurat_begin.R object final default assay: ", DefaultAssay(s))
 
 # Save the R object
-saveRDS(s, file=file.path(opt$outdir, "begin.rds"))
+if (opt$file_type == "rds") {
+  outfile = file.path(opt$outdir, "begin.rds")
+} else {
+  outfile = file.path(opt$outdir, "begin.h5seurat")
+}
+saveSeurat(path=outfile, format=opt$file_type)
 
 message("Completed")

R/seurat_characteriseClusterDEGenes.R

@@ -187,7 +187,8 @@ tex <- c(tex, getFigureTex(defn, deCaption,plot_dir_var=opt$plotdirvar))
 ## enforce significance
 
 ## read in the seurat object
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
+
 cluster_ids <- readRDS(opt$clusterids)
 Idents(s) <- cluster_ids
 

R/seurat_clusterStats.R

@@ -52,7 +52,7 @@ cat("Running with options:\n")
 print(opt)
 
 
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 cluster_ids <- readRDS(opt$clusterids)
 
 xx <- names(cluster_ids)

R/seurat_compare_clusters.R

@@ -36,8 +36,7 @@ opt <- parse_args(OptionParser(option_list=option_list))
 cat("Running with options:\n")
 print(opt)
 
-message(sprintf("readRDS: %s", opt$seuratobject))
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 
 message("seurat_cluster.R running with default assay: ", DefaultAssay(s))
 

R/seurat_dm.R

@@ -48,9 +48,7 @@ opt <- parse_args(OptionParser(option_list=option_list))
 cat("Running with options:\n")
 print(opt)
 
-
-message("readRDS")
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 cluster_ids <- readRDS(opt$clusterids)
 Idents(s) <- cluster_ids
 

R/seurat_explore_hvg_and_cell_cycle.R

@@ -200,7 +200,7 @@ options(future.globals.maxSize = opt$memory * 1024^2)
 
 # Input data ----
 
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 
 ## ######################################################################### ##
 ## # (i) Initial normalisation, variable gene identification and scaling ### ##

R/seurat_find_neighbors.R

@@ -87,7 +87,7 @@ options(future.globals.maxSize = opt$memory * 1024^2)
 
 # Input data ----
 
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 
 if(opt$usesigcomponents)
 {

R/seurat_get_assay_data.R

@@ -63,8 +63,7 @@ cat("Running with options:\n")
 print(opt)
 
 
-## read in the raw count matrix
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 
 ## set the default assay
 message("Setting default assay to: ", opt$seuratassay)

R/seurat_jackstraw.R

@@ -30,7 +30,7 @@ opt <- parse_args(OptionParser(option_list=option_list))
 cat("Running with options:\n")
 print(opt)
 
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 
 ## In Macosko et al, we implemented a resampling test inspired by the jackStraw procedure.
 ## We randomly permute a subset of the data (1% by default) and rerun PCA,

R/seurat_normalise_and_scale.R

@@ -221,7 +221,7 @@ options(future.globals.maxSize = opt$memory * 1024^2)
 
 # Input data ----
 
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 
 ## ######################################################################### ##
 ## # (i) Initial normalisation, variable gene identification and scaling ## ##
@@ -351,6 +351,6 @@ if(opt$normalizationmethod=="sctransform")
 message("seurat_begin.R object final default assay: ", DefaultAssay(s))
 
 # Save the R object
-saveRDS(s, file=file.path(opt$outdir, "begin.rds"))
+saveSeurat(path=opt$seuratobject)
 
 message("Completed")

R/seurat_pca.R

@@ -115,7 +115,7 @@ options(future.globals.maxSize = opt$memory * 1024^2)
 
 # Input data ----
 
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 
 ## ######################################################################### ##
 ## ################### (vi) Dimension reduction (PCA) ###################### ##
@@ -200,6 +200,6 @@ if(opt$normalizationmethod!="sctransform" & opt$jackstraw)
 message("seurat_begin.R object final default assay: ", DefaultAssay(s))
 
 # Save the R object
-saveRDS(s, file=file.path(opt$outdir, "begin.rds"))
+saveSeurat(path=opt$seuratobject)
 
 message("Completed")

R/seurat_summariseMarkers.R

@@ -40,7 +40,8 @@ if(!is.null(opt$subgroup)) { opt$subgroup <- strsplit(opt$subgroup,",")[[1]]}
 cat("Running with options:\n")
 print(opt)
 
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
+
 cluster_ids <- readRDS(opt$clusterids)
 Idents(s) <- cluster_ids
 

R/seurat_summariseMarkersBetween.R

@@ -48,7 +48,7 @@ print(opt)
 # set the run specs
 run_specs <- paste(opt$numpcs,opt$resolution,opt$algorithm,opt$testuse,sep="_")
 
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
 cluster_ids <- readRDS(opt$clusterids)
 
 message("Setting the default seurat assay to: ", opt$seuratassay)

R/seurat_tsne.R

@@ -45,8 +45,8 @@ cat("Running with options:\n")
 print(opt)
 
 
-message("readRDS")
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
+
 cluster_ids <- readRDS(opt$clusterids)
 Idents(s) <- cluster_ids
 print(head(Idents(s)))

R/seurat_umap.R

@@ -50,9 +50,8 @@ opt <- parse_args(OptionParser(option_list=option_list))
 cat("Running with options:\n")
 print(opt)
 
+s <- loadSeurat(path=opt$seuratobject)
 
-message("readRDS")
-s <- readRDS(opt$seuratobject)
 s@graphs <- readRDS(opt$seuratgraphs)
 
 message("seurat_umap.R running with default assay: ", DefaultAssay(s))

R/singleR_extra_plots.R

@@ -35,8 +35,8 @@ opt <- parse_args(OptionParser(option_list=option_list))
 cat("Running with options:\n")
 print(opt)
 
-message(sprintf("readRDS: %s", opt$seuratobject))
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
+
 a <- ifelse("SCT" %in% names(s), yes = "SCT", no = "RNA")
 
 sce <- as.SingleCellExperiment(s, assay = a) # Seems to be using SCT assay (so does it take the defualt?)

R/singleR_plots.R

@@ -38,8 +38,8 @@ opt <- parse_args(OptionParser(option_list=option_list))
 cat("Running with options:\n")
 print(opt)
 
-message(sprintf("readRDS: %s", opt$seuratobject))
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
+
 a <- ifelse("SCT" %in% names(s), yes = "SCT", no = "RNA")
 
 sce <- as.SingleCellExperiment(s, assay = a) # Seems to be using SCT assay (so does it take the defualt?)

R/singleR_run.R

@@ -35,8 +35,8 @@ opt <- parse_args(OptionParser(option_list=option_list))
 cat("Running with options:\n")
 print(opt)
 
-message(sprintf("readRDS: %s", opt$seuratobject))
-s <- readRDS(opt$seuratobject)
+s <- loadSeurat(path=opt$seuratobject)
+
 a <- ifelse("SCT" %in% names(s), yes = "SCT", no = "RNA")
 
 sce <- as.SingleCellExperiment(s, assay = a) # Seems to be using SCT assay (so does it take the defualt?)