Merge pull request #27 from mari-ga/main

Correction for old merge
theislab · Jan 8, 2024 · c4736af · c4736af
2 parents 94067fd + 1a243c0
commit c4736af
Show file tree

Hide file tree

Showing 23 changed files with 676 additions and 972 deletions.
diff --git a/.gitignore b/.gitignore
@@ -6,4 +6,5 @@ result*/
 testing/
 testing*
 *.pyc
-docs/build/
+docs/build/
+nextflow_internal.config
diff --git a/bin/bff.R b/bin/bff.R
@@ -74,8 +74,8 @@ if(as.logical(args$preprocess)){
   # Step 3: Create a vector from the barcodesl
   vector <- unlist(words)
   print("Preprocessing")
-  #counts <- Read10X(args$fileHto) 
   counts <- ProcessCountMatrix(rawCountData = args$fileHto, barcodeBlacklist = vector)
+  print("Preprocessing done")
 }else{
   print("No preprocessing")
   counts <- Read10X(args$fileHto) 
@@ -89,19 +89,6 @@ if (!is.null(args$methodsForConsensus)) {
 
 perCell_args <- args$perCellSaturation
 perCell <- ifelse(perCell_args == "null" || perCell_args == "Null", NULL, perCell_args)
-<<<<<<< HEAD
-<<<<<<< HEAD
-print("---------------------")
-print(perCell)
-print("---------------------")
-=======
-print(perCell)
-
->>>>>>> c781241 (bff re-added problematic parameter)
-=======
-print("---------------------")
-print(perCell)
-print("---------------------")
 
 if(args$methodsForConsensus=="bff_raw" || args$methodsForConsensus=="bff_cluster" || args$methodsForConsensus=="bff_raw,bff_cluster" || is.null(args$methodsForConsensus)  )
   #Only Bff in its different variations is available

diff --git a/bin/demuxem.py b/bin/demuxem.py
@@ -33,6 +33,7 @@
     # load input rna data
     rna_data = sc.read_10x_mtx(args.rna_matrix_dir)
     hashing_data = sc.read_10x_mtx(args.hto_matrix_dir,gex_only=False)
+    rna = args.rna_matrix_dir
     filter = ""
     if args.filter_demuxem.lower() in ['true', 't', 'yes', 'y', '1']:
         filter = True

diff --git a/bin/demuxmix.R b/bin/demuxmix.R
diff --git a/bin/dropletUtils.R b/bin/dropletUtils.R
@@ -27,9 +27,11 @@ parser$add_argument("--objectOutEmptyDrops", default = "emptyDroplets",
                     help = "Prefix name for the emptyDrops RDS file")
 parser$add_argument("--assignmentOutEmptyDrops", default = "emptyDroplets",
                     help = "prefex name for emptyDrops assignment CSV file")
+parser$add_argument("--runEmptyDrops", action="store_false",
+                    help = "Executes emptyDrops function only when desired, recomended only for raw data")
 
 #for hashedDrops
-parser$add_argument("--ambient", action = "store_true",
+parser$add_argument("--ambient", action = "store_false",
                     help = "Whether to use the relative abundance of each HTO in the ambient solution from emtpyDrops, set TRUE only when test_ambient is TRUE.")
 parser$add_argument("--minProp", default = 0.05, type = "double",
                     help = "Numeric scalar to be used to infer the ambient profile when ambient=NULL,")
@@ -60,36 +62,40 @@ parser$add_argument("--gene_col", help = "Specify which column of genes.tsv or f
 args <- parser$parse_args()
 
 hto <- Read10X(data.dir = args$raw_hto_matrix_dir, gene.column = args$gene_col)
-
-emptyDrops_out <- emptyDrops(hto, lower = args$lower, niters = args$niters,
-                             test.ambient = args$testAmbient,
-                             ignore = args$ignore,
-                             alpha = args$alpha, round = args$round,
-                             by.rank = args$byRank)
-
-print("------------------- emptyDrops finished ---------------------------------")
+combinations_transformed <- ifelse(tolower(args$combinations) == "null", NULL, args$combinations)
+
+if (args$runEmptyDrops == TRUE) {
+    rna <- Read10X(data.dir = args$raw_rna_matrix_dir,gene.column = args$gene_col)
+    print("------------------- executing emptyDrops ---------------------------------")
+    ignore_transformed <- ifelse(tolower(args$ignore) == "null", NULL, args$ignore)
+    emptyDrops_out <- emptyDrops(rna, lower = args$lower, niters = args$niters,
+                                test.ambient = args$testAmbient,
+                                ignore = NULL,
+                                alpha = args$alpha, round = args$round,
+                                by.rank = args$byRank)
+
+
+    write.csv(emptyDrops_out, paste0(args$outputdir, "/", args$assignmentOutEmptyDrops, ".csv"))
+    saveRDS(emptyDrops_out, file=paste0(args$outputdir, "/", args$objectOutEmptyDrops, ".rds"))
+
+    print("------------------- filtering empty droplets ----------------------------")
+    is.cell <- emptyDrops_out$FDR <= args$isCellFDR
+    colors <- ifelse(is.cell, "red", "black")
+    png(paste0(args$outputdir, "/", "plot_emptyDrops.png"))
+    plot(emptyDrops_out$Total, -emptyDrops_out$LogProb, col=colors, xlab="Total UMI count", ylab="-Log Probability")
+    dev.off()
 
 
-print("-------- Following Files are saved in folder hashedDrops_out ------------")
-print(paste0(args$objectOutEmptyDrops, ".rds"))
-print(paste0(args$assignmentOutEmptyDrops, ".csv"))
-write.csv(emptyDrops_out, paste0(args$outputdir, "/", args$assignmentOutEmptyDrops, ".csv"))
-saveRDS(emptyDrops_out, file=paste0(args$outputdir, "/", args$objectOutEmptyDrops, ".rds"))
-
-print("------------------- filtering empty droplets ----------------------------")
-is.cell <- emptyDrops_out$FDR <= args$isCellFDR
-colors <- ifelse(is.cell, "red", "black")
-png(paste0(args$outputdir, "/", "plot_emptyDrops.png"))
-plot(emptyDrops_out$Total, -emptyDrops_out$LogProb, col=colors, xlab="Total UMI count", ylab="-Log Probability")
-dev.off()
-
-if (args$ambient == TRUE) {
-    hashedDrops_out <- hashedDrops(hto[,which(is.cell)], min.prop = args$minProp, ambient = metadata(emptyDrops_out)$ambient, pseudo.count = args$pseudoCount, constant.ambient = args$constantAmbient, doublet.nmads = args$doubletNmads, doublet.min = args$doubletMin, doublet.mixture = args$doubletMixture, confident.nmads = args$confidentNmads, confident.min = args$confidenMin, combinations = args$combinations)
+    if (args$ambient == TRUE) {
+        hashedDrops_out <- hashedDrops(hto[,which(is.cell)], min.prop = args$minProp, ambient = metadata(emptyDrops_out)$ambient, pseudo.count = args$pseudoCount, constant.ambient = args$constantAmbient, doublet.nmads = args$doubletNmads, doublet.min = args$doubletMin, doublet.mixture = args$doubletMixture, confident.nmads = args$confidentNmads, confident.min = args$confidenMin, combinations = combinations_transformed)
+    } else {
+        hashedDrops_out <- hashedDrops(hto[,which(is.cell)], min.prop = args$minProp, pseudo.count = args$pseudoCount, constant.ambient = args$constantAmbient, doublet.nmads = args$doubletNmads, doublet.min = args$doubletMin, doublet.mixture = args$doubletMixture, confident.nmads = args$confidentNmads, confident.min = args$confidenMin, combinations = combinations_transformed)
+    }
 } else {
-    hashedDrops_out <- hashedDrops(hto[,which(is.cell)], min.prop = args$minProp, pseudo.count = args$pseudoCount, constant.ambient = args$constantAmbient, doublet.nmads = args$doubletNmads, doublet.min = args$doubletMin, doublet.mixture = args$doubletMixture, confident.nmads = args$confidentNmads, confident.min = args$confidenMin, combinations = args$combinations)
+    hashedDrops_out <- hashedDrops(hto,min.prop = args$minProp,pseudo.count = args$pseudoCount, constant.ambient = args$constantAmbient,doublet.nmads = args$doubletNmads, doublet.min = args$doubletMin,confident.nmads = args$confidentNmads,confident.min = args$confidenMin)
+
 }
 
-print("------------------- hashedDrops finished ---------------------------------")
 
 ignore <- args$ignore
 if (is.null(ignore)) {

diff --git a/bin/generate_data.py b/bin/generate_data.py
@@ -26,5 +26,17 @@
         adata.write("adata_with_donor_matching.h5ad")
 
     if args.generate_mudata:
-        # write mudata_with_donor_matching.h5mu data
-        pass
+        rna_data = sc.read_10x_mtx(args.read_rna_mtx)
+        hto_data = sc.read_10x_mtx(args.read_hto_mtx, gex_only=False)
+        assignment_dir = os.path.join(args.assignment, 
+                                    [filename for filename in os.listdir(args.assignment) if filename == "all_assignment_after_match.csv"][0])
+
+        assignment = pd.read_csv(assignment_dir, index_col = 0)
+        mudata = MuData({"rna": rna_data, "hto": hto_data })
+
+        mudata['rna'].obs = mudata['rna'].obs.merge(args.assignment, left_index=True, right_index=True, how='left')
+        mudata['rna'].obs.rename(columns={mudata['rna'].obs.columns[0]: 'donor'}, inplace=True)
+        mudata['rna'].obs.donor = mudata['rna'].obs.donor.fillna("negative")
+        mudata['rna'].obs.donor = mudata['rna'].obs.donor.astype(str)
+        mudata.update()
+        mudata.write("mudata_with_donor_matching.h5mu")