Added real bulk with no ground truth

Author: zgr2788

.gitignore

@@ -6,7 +6,6 @@ Output
 Scratch
 .snakemake
 *pbulks*
-*props*
 notes
 log.txt
 C_1*

Modules/ADAPTS/proportionsInAdmixture_run.R

@@ -38,7 +38,7 @@ res[is.na(res)] <- 0
 res = apply(res,2,function(x) ifelse(x < 0, 0, x)) #explicit non-negativity constraint
 res = apply(res,2,function(x) x/sum(x)) #explicit STO constraint
 res <- res[-nrow(res),]
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/BisqueRNA/Bisque_run.R

@@ -38,7 +38,7 @@ C0$SubjectName <- C0$sampleID #BisqueRNA uses hard-coded colname
 
 #Get results and reorder the matrices for correspondence
 res <- BisqueRNA::ReferenceBasedDecomposition(T, C0, markers=NULL, use.overlap=FALSE)$bulk.props
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/CDSeq/CDSeq_run.R

@@ -31,7 +31,7 @@ T  <- T[common,]
 #Get results and reorder the matrices for correspondence
 res <- CDSeq(bulk_data = T, reference_gep = C1, cell_type_number = ncol(C1), mcmc_iterations = 1000, cpu_number = cores, block_number = 6, gene_subset_size=15)
 res <- res$estProp
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/CIBERSORT/CIBERSORT_run.R

@@ -33,7 +33,7 @@ T <- data.frame(T)
 #Run deconv, params may be exported to config later
 res <- CIBERSORT(sig_matrix = C1, mixture_file = T, QN = FALSE, absolute = FALSE, perm = 0)
 res <- t(res[,1:(ncol(res)-3)])
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/CPM/CPM_run.R

@@ -40,7 +40,7 @@ cellTypes <- as.character(phenData$cellType)
 res <- t(CPM(C0, cellTypes, T, p, quantifyTypes = TRUE, no_cores = cores)$cellTypePredictions)
 res <- apply(res,2,function(x) ifelse(x < 0, 0, x)) #explicit non-negativity constraint
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/CellMix/DSA_run.R

@@ -34,7 +34,7 @@ [email protected] <- tapply(as.character(C2$geneName),as.character(C2$cellType),list)
 res <- CellMix::ged(as.matrix(T), ML, method = "DSA", log = FALSE)@fit@H
 res <- apply(res,2,function(x) ifelse(x < 0, 0, x)) #explicit non-negativity constraint
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/CellMix/deconf_run.R

@@ -30,7 +30,7 @@ [email protected] <- tapply(as.character(C2$geneName),as.character(C2$cellType),list)
 #Get res and reorder the matrices for correspondence
 res <- CellMix::ged(T, x=length(unique(as.character(C2$cellType))), method = "deconf", maxIter = 500, verbose= TRUE)
 res <- CellMix::match.nmf(res, ML)@fit@H
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/CellMix/ssFrobenius_run.R

@@ -29,7 +29,7 @@ [email protected] <- tapply(as.character(C2$geneName),as.character(C2$cellType),list)
 
 #Get res and reorder the matrices for correspondence
 res <- CellMix::ged(as.matrix(T), ML, method = "ssFrobenius", sscale = TRUE, maxIter = 500, log = FALSE)@fit@H
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/CellMix/ssKL_run.R

@@ -31,7 +31,7 @@ [email protected] <- tapply(as.character(C2$geneName),as.character(C2$cellType),list)
 #Get res and reorder the matrices for correspondence
 res <- CellMix::ged(as.matrix(T), x=length(unique(as.character(C2$cellType))), method = "ssKL",markers= "semi", sscale = FALSE, maxIter=500, log = FALSE, verbose= TRUE)
 res <- CellMix::match.nmf(res, ML)@fit@H
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/DCQ/DCQ_run.R

@@ -32,7 +32,7 @@ T  <- T[common,]
 res <- t(dcq(reference_data = C1, mix_data = T, marker_set = as.data.frame(row.names(C1)) , alpha_used = 0.05, lambda_min = 0.2, number_of_repeats = 10)$average)
 res <- apply(res,2,function(x) ifelse(x < 0, 0, x)) #explicit non-negativity constraint
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/DWLS/DWLS_run.R

@@ -47,7 +47,7 @@ res <- future_apply(T,2, function(x){
 }, future.seed = TRUE)
 
 rownames(res) <- as.character(unique(phenData$cellType))
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/DeconRNASeq/DeconRNASeq_run.R

@@ -29,7 +29,7 @@ C1 <- data.frame(C1)
 #Run DeconRNASeq and reorder rows for correspondence
 res <- DeconRNASeq(datasets = T, signatures = C1, proportions = NULL, checksig = FALSE, known.prop = FALSE, use.scale = TRUE, fig = FALSE)
 res <- t(res$out.all)
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/EPIC/EPIC_run.R

@@ -49,7 +49,7 @@ C_EPIC[["refProfiles.var"]] <- refVar[markers, cellTypes]
 res <- t(EPIC(bulk=as.matrix(T), reference=C_EPIC, withOtherCells=TRUE, scaleExprs=FALSE)$cellFractions) #scaleExprs=TRUE by default: only keep genes in common between matrices
 res <- res[!rownames(res) %in% "otherCells",]
 res[is.na(res)] <- 0
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/EpiDISH/EpiDISH_run.R

@@ -30,7 +30,7 @@ C1 <- as.matrix(C1) #Explicit typecasting needed
 res <- t(EpiDISH::epidish(beta.m = T, ref.m = C1, method = "RPC")$estF)
 res <- apply(res,2,function(x) ifelse(x < 0, 0, x)) #explicit non-negativity constraint
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/FARDEEP/FARDEEP_run.R

@@ -31,7 +31,7 @@ message('Started running')
 res <- fardeep(C1, T, nn = TRUE, intercept = TRUE, permn = 10, QN = FALSE)
 res <- t(res$abs.beta)
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/MuSiC/MuSiC_run.R

@@ -52,7 +52,7 @@ C0 <- ExpressionSet(C0, phenoData = as(phenData, "AnnotatedDataFrame"))
 
 #Get results and reorder the matrices for correspondence
 res <- t(music_prop(bulk.eset = T, sc.eset = C0, clusters = "cellType", samples = "sampleID", select.ct = levels(C0@phenoData$cellType), markers = NULL, normalize = FALSE, verbose = FALSE)$Est.prop.weighted)
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/NNLS/NNLS_run.R

@@ -26,7 +26,7 @@ T  <- T[common,]
 res <- do.call(cbind.data.frame,lapply(apply(T,2,function(x) nnls::nnls(as.matrix(C1),x)), function(y) y$x))
 rownames(res) <- colnames(C1)
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/OLS/OLS_run.R

@@ -27,7 +27,7 @@ rownames(res) <- unlist(lapply(strsplit(rownames(res),")"),function(x) x[2])) #C
 res[is.na(res)] <- 0       #Convert NA's to zeros prior to applying constraints
 res = apply(res,2,function(x) ifelse(x < 0, 0, x)) #Explicit non-negativity constraint
 res = apply(res,2,function(x) x/sum(x)) #Explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/Psuedobulk/dummy_props.rds

@@ -28,7 +28,7 @@ res <- do.call(cbind.data.frame,lapply(apply(T,2,function(x) MASS::rlm(x ~ as.ma
 rownames(res) <- unlist(lapply(strsplit(rownames(res),")"),function(x) x[2]))
 res <- apply(res,2,function(x) ifelse(x < 0, 0, x)) #explicit non-negativity constraint
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 
 #Save and exit

Modules/RLR/RLR_run.R

@@ -23,6 +23,9 @@ filenames <- list.files("Output") %>% lapply(., function(x) { x <- paste0("Outpu
 files <- lapply(filenames, readRDS)
 names(files) <- lapply(filenames, function(x) gsub("Output/.*res_(.*).rds", "\\1", x))
 
+#Reorder files in case this was without proportions
+files <- lapply(files, function(x) x[order(match(rownames(x), rownames(files[[1]]))),])
+
 #Generate combs to iterate over
 combs <- combn(files, 2)
 colnames(combs) <- combn(names(files), 2, paste0, collapse="@")

Modules/Res_explore/Heatmap_agr.R

@@ -12,8 +12,8 @@ if (!require("L1pack", quietly = TRUE)){
 
 suppressMessages(library(remotes))
 suppressMessages(library(devtools))
-remotes::install_github("renozao/xbioc", auth_token = "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")
 options(timeout=400) # to avoid problems when connection slow
+remotes::install_github("renozao/xbioc", auth_token = "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")
 devtools::install_github("meichendong/SCDC", auth_token = "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")
 suppressMessages(library(SCDC))
 suppressMessages(library(Biobase))
@@ -47,7 +47,7 @@ C0 <- ExpressionSet(C0, phenoData = as(phenData, "AnnotatedDataFrame"))
 
 #Get results and reorder the matrices for correspondence
 res <- t(SCDC::SCDC_prop(bulk.eset = T, sc.eset = C0, ct.varname = "cellType", sample = "sampleID", ct.sub = levels(phenData$cellType), iter.max = 200)$prop.est.mvw)
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/SCDC/SCDC_run.R

@@ -38,7 +38,7 @@ names(ref_anno) <- phenData$cellType
 
 #Get results and reorder the matrices for correspondence
 res <-t(TIMER_deconv(mix = T, ref = C0, curated.cell.types = ref_anno, sig = rownames(T)))
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/TIMER/TIMER_run.R

@@ -57,7 +57,7 @@ B <- bseqsc_basis(C0, markers, clusters = 'cellType', samples = 'sampleID', ct.s
 message("Running CIBERSORT with B...")
 res <- CIBERSORT(B, T, QN = FALSE, absolute = FALSE, perm = 0)
 res <- t(res[,1:(ncol(res)-3)])
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/bseqsc/bseqsc_run.R

@@ -5,11 +5,11 @@
 ##
 ## Description:
 ## This file includes all rules associated with deconvolution using debCAM
-## package. The 3 different modes are
+## package
+##
+##      Semi-supervised with C2 reference
+##
 ##
-##      1- Unsupervised with k coefficient
-##      2- Semi-supervised with C1 reference
-##      3- Semi-supervised with C2 reference
 ##
 ##
 ##
@@ -18,33 +18,6 @@
 ##
 ##
 ##
-
-rule debCAM_unsupervised:
-    input:
-        getBulks('{sample}')
-
-    output:
-        "Output/{sample}_res_debCAM_unsupervised.rds"
-
-    params:
-        cellTypes_n = config['unsupervisedParams']['expectedTypes']
-
-    conda:
-        "env.yaml"
-
-    threads: config["cores"]["debCAM_unsupervised"]
-
-    priority: -50
-
-    benchmark:
-        "Benchmarks/{sample}_debCAM_unsupervised_benchmark.txt"
-
-    resources:
-        mem_mb=getMB(config['mem_gb']['minorOps'])
-
-    shell:
-        "Rscript Modules/debCAM/debCAM_u.R {input} {params} {output} {threads}"
-
 
 rule debCAM_marker:
     input:
@@ -67,24 +40,3 @@ rule debCAM_marker:
 
     shell:
         "Rscript Modules/debCAM/debCAM_marker.R {input} {params} {output}"
-
-
-
-rule debCAM_C1:
-    input:
-        getC1(getBulks('{sample}'), '{sample}')
-
-    output:
-        "Output/{sample}_res_debCAM_C1.rds"
-
-    conda:
-        "env.yaml"
-
-    benchmark:
-        "Benchmarks/{sample}_debCAM_C1_benchmark.txt"
-
-    resources:
-        mem_mb=getMB(config['mem_gb']['minorOps'])
-
-    shell:
-        "Rscript Modules/debCAM/debCAM_C1.R {input} {output}"

Modules/debCAM/Snakefile

@@ -25,7 +25,7 @@ P <- readRDS(filename_P)
 
 #Switch format to debCAM MGlist
 MGlist <- C2 %>% split(., .[, 'cellType']) %>% lapply(., function(x) x[order(-abs(x$log2FC)),]) %>% lapply(., function(x) rownames(x))
-MGlist <- MGlist[order(match(names(MGlist), rownames(P)))] #Makes sure list order is same as rownames(P) for future metrics
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') MGlist <- MGlist[order(match(names(MGlist), rownames(P)))] #Makes sure list order is same as rownames(P) for future metrics
 
 #Run deconv with marker list generated by Seurat
 res <- redoASest(T, MGlist, maxIter = 50) #Adjust maxIter later

Modules/debCAM/debCAM_marker.R

@@ -41,7 +41,7 @@ MD <- lapply(MD,function(x) sapply(x, function(y) which(y==rownames(C1))))
 
 #Get results and reorder the matrices for correspondence
 res <- t(dtangle::dtangle(Y=mixtures, reference=refs, markers = MD)$estimates)
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/dtangle/dtangle_run.R

@@ -28,7 +28,7 @@ T  <- T[common,]
 res <- apply(T, 2, function(z) coef(glmnet::glmnet(x = as.matrix(C1), y = z, alpha = 0.2, standardize = TRUE, lambda = glmnet::cv.glmnet(as.matrix(C1), z)$lambda.1se))[1:ncol(C1)+1,])
 res <- apply(res,2,function(x) ifelse(x < 0, 0, x)) #explicit non-negativity constraint
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/elasticNET/elasticNET_run.R

@@ -30,7 +30,7 @@ res <- apply(T, 2, function(z) coef(glmnet::glmnet(x = as.matrix(C1), y = z, alp
 res <- apply(res,2,function(x) ifelse(x < 0, 0, x)) #explicit non-negativity constraint
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
 res[is.na(res)] <- 0 #Needed for models where glmnet drops all terms of a model and fit an intercept-only model (very unlikely but possible).
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/lasso/lasso_run.R

@@ -31,7 +31,7 @@ cellTypes <- phenData$cellType
 model <- omnideconv::build_model(C0, cellTypes, bulk_gene_expression = T, method = "autogenes")
 res <- omnideconv::deconvolute(T, signature = model, single_cell_object = C0, method = "autogenes")
 res <- t(res/rowSums(res))
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/omnideconv/Autogenes_run.R

@@ -30,7 +30,7 @@ rm(common)
 cellTypes <- phenData$cellType
 model <- omnideconv::build_model(C0, cellTypes, bulk_gene_expression = T, method = "momf")
 res <- t(omnideconv::deconvolute(T, signature = model, single_cell_object = C0, method = "momf"))
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/omnideconv/MOMF_run.R

@@ -29,7 +29,7 @@ T  <- T[common,]
 res <- apply(T, 2, function(z) coef(glmnet::glmnet(x = as.matrix(C1), y = z, alpha = 0, standardize = TRUE, lambda = glmnet::cv.glmnet(as.matrix(C1), z)$lambda.1se))[1:ncol(C1)+1,])
 res <- apply(res,2,function(x) ifelse(x < 0, 0, x)) #explicit non-negativity constraint
 res <- apply(res,2,function(x) x/sum(x)) #explicit STO constraint
-res <- res[order(match(rownames(res), rownames(P))),]
+if (filename_P != 'Modules/Psuedobulk/dummy_props.rds') res <- res[order(match(rownames(res), rownames(P))),]
 
 #Save and exit
 saveRDS(res, file=filename_O)

Modules/ridge/ridge_run.R

@@ -26,6 +26,8 @@ cat: bug
 </pre>
 
 # Running the pipeline
+**IMPORTANT**: If running for the first time, use one sample only as the environments are
+installed throughout the workflow
 ## Self-reference
 ### Description
 Uses **one** single-cell reference to generate the pseudobulks and references for deconvolution benchmarking. Important assumptions:
@@ -57,9 +59,9 @@ Uses **one** single-cell reference to generate the pseudobulks and references fo
 		mamba create -n snakemake snakemake
 
  4. Build the pipeline.
- 		
+
 		make
-		
+
  5.  Place the input file in the `Input` directory.
  6. Adjust settings in `config.yaml`.
  7. **(Optional)** Run `getDag.sh` to generate the updated DAG after adjusting config.
@@ -114,7 +116,7 @@ Same as self-reference, except after the **3rd** step, note the following direct
 
 	Input/Psuedobulks
 
-This is where the bulks should go. Also, **enable realBulk in config.** If you wish to use all the data for the reference, go to: 
+This is where the bulks should go. Also, **enable realBulk in config.** If you wish to use all the data for the reference, go to:
 
 	Input/Cell_splits