just wrap labelTransfer

wbaopaul · Jun 30, 2021 · b6380d9 · b6380d9
1 parent a08a2b2
commit b6380d9
Show file tree

Hide file tree

Showing 4 changed files with 1,419 additions and 42 deletions.
diff --git a/README.md b/README.md
@@ -51,7 +51,7 @@ Updates
 - Now provide [scATAC-pro tutorial in R](https://scatacpro-in-r.netlify.app/index.html) for access QC metrics and perform downstream analysis
 - Current version: 1.4.0
 - Recent updates
-    * *labelTransfer*: new module added, to do label trasfer (for cell annotation) from scRNA-seq annotation. First construct a gene by cell activity matrix, then use *FindTransferAnchors* and *TransferData* function from Seurat R package to predicted cell type annotation from the cell annotaiton in scRNA-seq data.
+    * *labelTransfer*: new module added, to do label trasfer (for cell annotation) from cell annotation of scRNA-seq data. First construct a gene by cell activity matrix, then use *FindTransferAnchors* and *TransferData* function from Seurat R package to predicted cell type annotation from the cell annotaiton in scRNA-seq data.
     * *rmDoublets*: new module added, to remove potential doublets using [DoubletFinder](https://github.com/chris-mcginnis-ucsf/DoubletFinder) algorithm.
     * *qc_per_barcode*: add tss enrichment score per cell into the QC metrics
     * *call_cell*: enable filtering barcodes with minimal tss enrichment score cutoff (parameter **min_tss_escore** in the updated [configure_user.txt](configure_user.txt) file)

diff --git a/scripts/src/dsAnalysis_utilities.R b/scripts/src/dsAnalysis_utilities.R
@@ -1318,3 +1318,51 @@ generate_gene_cisActivity <- function(gene_ann, mtx, include_body = T,
 }
 
 generate_gene_cisActivity = cmpfun(generate_gene_cisActivity)
+
+
+labelTransfer_R <- function(seurat.atac, seurat.rna, gene_ann,
+                            rna_ann_var = 'Cell_Type', include_genebody = T){
+  atac.mtx = seurat.atac@assays$ATAC@counts
+  rn = rownames(atac.mtx)
+  rownames(atac.mtx) <- sapply(rn, function(x) unlist(strsplit(x, ','))[1])
+  activity.matrix = generate_gene_cisActivity(gene_ann = gene_ann,
+                                              atac.mtx, 
+                                              include_body = include_genebody)
+
+  seurat.atac[["ACTIVITY"]] <- CreateAssayObject(counts = activity.matrix)
+
+  DefaultAssay(seurat.atac) <- "ACTIVITY"
+  seurat.atac <- NormalizeData(seurat.atac)
+  seurat.atac <- ScaleData(seurat.atac, features = rownames(seurat.atac))
+  seurat.atac <- FindVariableFeatures(seurat.atac)
+
+  DefaultAssay(seurat.atac) <- "ATAC"
+  seurat.atac$tech = 'ATAC'
+  seurat.rna$tech = 'RNA'
+
+  ## transfer label 
+  genes4anchors = VariableFeatures(object = seurat.rna)
+  #genes4anchors = NULL
+  transfer.anchors <- FindTransferAnchors(reference = seurat.rna,
+                                          query = seurat.atac,
+                                          features = genes4anchors,
+                                          reference.assay = "RNA",
+                                          query.assay = "ACTIVITY",
+                                          reduction = "cca",
+                                          k.anchor = 5)
+
+
+  celltype.predictions <- TransferData(anchorset = transfer.anchors,
+                                       refdata = seurat.rna[[rna_ann_var]][, 1],
+                                       weight.reduction = seurat.atac[["pca"]],
+                                       dims = 1:ncol(seurat.atac[["pca"]]),
+                                       k.weight = 50)
+  celltype.predictions = subset(celltype.predictions, select = c('predicted.id', 'prediction.score.max'))
+  names(celltype.predictions)[1] = 'Predicted_Cell_Type'
+  seurat.atac <- AddMetaData(seurat.atac, metadata = celltype.predictions)
+  rm(transfer.anchors)
+
+  seurat.atac[["ACTIVITY"]] <- NULL ## don't save activity assay
+  return(seurat.atac)
+}
+labelTransfer_R = cmpfun(labelTransfer_R)
diff --git a/scripts/src/labelTransfer.R b/scripts/src/labelTransfer.R
@@ -124,54 +124,15 @@ if(F){
 seurat.rna = readRDS(inputSeurat_rna)
 seurat.atac = readRDS(inputSeurat_atac)
 
-atac.mtx = seurat.atac@assays$ATAC@counts
-rn = rownames(atac.mtx)
-rownames(atac.mtx) <- sapply(rn, function(x) unlist(strsplit(x, ','))[1])
-activity.matrix = generate_gene_cisActivity(gene_ann = gene_ann,
-                                            atac.mtx, 
-                                            include_body = T)
-
-seurat.atac[["ACTIVITY"]] <- CreateAssayObject(counts = activity.matrix)
-
-DefaultAssay(seurat.atac) <- "ACTIVITY"
-seurat.atac <- NormalizeData(seurat.atac)
-seurat.atac <- ScaleData(seurat.atac, features = rownames(seurat.atac))
-seurat.atac <- FindVariableFeatures(seurat.atac)
-
-DefaultAssay(seurat.atac) <- "ATAC"
-seurat.atac$tech = 'ATAC'
-seurat.rna$tech = 'RNA'
-
-## transfer label 
-genes4anchors = VariableFeatures(object = seurat.rna)
-#genes4anchors = NULL
-transfer.anchors <- FindTransferAnchors(reference = seurat.rna,
-                                        query = seurat.atac,
-                                        features = genes4anchors,
-                                        reference.assay = "RNA",
-                                        query.assay = "ACTIVITY",
-                                        reduction = "cca",
-                                        k.anchor = 5)
-
-
-celltype.predictions <- TransferData(anchorset = transfer.anchors,
-                                     refdata = seurat.rna$Cell_Type,
-                                     weight.reduction = seurat.atac[["pca"]],
-                                     dims = 1:ncol(seurat.atac[["pca"]]),
-                                     k.weight = 50)
-celltype.predictions = subset(celltype.predictions, select = c('predicted.id', 'prediction.score.max'))
-names(celltype.predictions)[1] = 'Predicted_Cell_Type'
-seurat.atac <- AddMetaData(seurat.atac, metadata = celltype.predictions)
-rm(transfer.anchors)
+seurat.atac = labelTransfer_R(seurat.atac, seurat.rna, gene_ann,
+                              rna_ann_var = 'Cell_Type', include_genebody = T)
 
 p1 <- DimPlot(seurat.atac, group.by = "Predicted_Cell_Type",
               label = TRUE, repel = TRUE) + ggtitle("scATAC-seq") 
 
 p2 <- DimPlot(seurat.rna, group.by = "Cell_Type",
               label = TRUE, repel = TRUE) + ggtitle("scRNA-seq") 
 
-seurat.atac[["ACTIVITY"]] <- NULL ## don't save activity assay
-
 outputPath = dirname(inputSeurat_atac)
 outputName = basename(inputSeurat_atac)
 saveRDS(seurat.atac, file = paste0(outputPath, '/updated_', outputName))