fix HVF bug

epigen · Feb 15, 2024 · 3ce524b · 3ce524b
1 parent e47124e
commit 3ce524b
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 5 deletions.
diff --git a/config/config.yaml b/config/config.yaml
@@ -1,11 +1,11 @@
 
-# process output of: https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/feature-bc-analysis, will be loaded with Seurat::Read10X
+# processed output of: https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/feature-bc-analysis, will be loaded with Seurat::Read10X
 
 # alwayse use absolute paths
 
 ##### RESOURCES #####
 mem: '32000'
-threads: 2
+threads: 1
 partition: 'shortq'
 
 ##### GENERAL #####
@@ -61,7 +61,6 @@ extra_metadata: ""
 split_by:
     batch: ['sampleA','sampleB']
 
-
 ##### FILTER #####
 # logical expression for filtering, all variables need to be in metadata (leave empty "" if not applicable)
 # example: "pass_QC == 'True'"
@@ -73,7 +72,7 @@ filter_expression: "mitoPercent <= 20"
 pseudobulk:
     by: ['patient', 'cellType', 'treatment'] # have to be metadata columns
     method: "sum" # options: "sum", "mean", or "median"
-    cell_count_th: 20 # pseudomqbulked samples with less cells than cell_count_th are removed
+    cell_count_th: 20 # pseudobulked samples with less cells than cell_count_th are removed
 
 ##### NORMALIZATION #####
 # using SCTransform defaults https://satijalab.org/seurat/reference/sctransform

diff --git a/workflow/scripts/seurat_plots.R b/workflow/scripts/seurat_plots.R
@@ -109,7 +109,7 @@ if (step=="CORRECTED"){
         features <- snakemake@config[["vis_features"]][[feature_list]]
     }else{
         if (feature_list=="HVG"){
-            HVG_df <- HVFInfo(object = data_object, assay = "SCT")
+            HVG_df <- HVFInfo(object = data_object[["SCT"]], selection.method = "sct")
             features <- rownames(HVG_df[order(-HVG_df$residual_variance),])[1:100]
         }else{
             # load plotting gene list
@@ -154,6 +154,15 @@ if (plot_type=="Heatmap"){
     # Heatmap specs
     height <- max(5, height_row*length(features) + 1)
 
+    ### downsample for efficient hierarchical clustering and plotting of the same number of cells per group (ident)
+    if (ncol(data_object)>30000){
+        # captures count of cells for the ident with the fewest (but at least 100)
+        maxcells  <- min(100,table(Idents(data_object)))
+        # downsample
+        set.seed(42)
+        data_object <- subset(data_object, downsample = maxcells)
+    }
+
     # cluster features (rows) using hclust
     tmp_data <- as.data.frame(GetAssayData(object = data_object, slot = slot, assay = assay))
     # subset by features