yaml

This package, which builds on a seurat object, provides functions for the analysis and visualization of single-cell spatial perturb-seq data. It includes functionality for annotating barcodes (representing perturbations) which facilites downstream analyses like differential expression, analysis of barcode spatial distribution patterns, and calling cellular neighbours and neighbourhoods for the exploration of cell-cell interactions.

Adds an image component to a Seurat object

Adapted from: https://github.com/STOmics/STCellbin

add_image_slice(obj)
## An object of class Seurat 
## 30280 features across 3757 samples within 2 assays 
## Active assay: Spatial (26280 features, 2000 variable features)
##  3 layers present: counts, data, scale.data
##  1 other assay present: BANKSY
##  2 dimensional reductions calculated: spatial, pca
##  1 image present: slice1

Creates a metadata column based on expression of a user-given list of barcodes, with values “Pos” or “Neg”

obj <- annotate_pos_neg(obj, gRNA_list)
## 
##  Neg  Pos 
## 3624  133

Creates barplot showing the % of barcode-positive cells in each identity class. Ident can be seurat cluster, or annotated cell type, etc

barplot_pos(obj)
## # A tibble: 2 × 2
##   seurat_clusters Percent_Pos
##   <fct>                 <dbl>
## 1 3                      3.28
## 2 12                     4.43

Calculates DEGs for Pos vs Neg cell populations and shows top upregulated genes for each population

Spatial plot of negative and positive cells

spatial_pos_neg(obj)

Frequency histogram of nearest neighbor distances for Pos points in a spatial region, highlighting distribution of distances

A left-skewed graph suggests clustering of Pos points while a right-skewed graph suggests a more dispersed pattern. A uniform graph could indicate random distribution without clustering or dispersion, while multimodal peaks could indicate multiple clusters.

nndist_pos_neg(obj)

Creates a metadata column in a Seurat object based on expression of a list of features, labeling each cell with its barcode, “Neg” or “Multiple”

obj <- annotate_gRNA(obj, gRNA_list)
## 
##      Multiple           Neg sgrna-C9orf72 sgrna-Cfap410     sgrna-clu 
##            10          3623             1             6            21 
##    sgrna-dpp5    sgrna-fasn    sgrna-flcn    sgrna-gfap   sgrna-lrrk2 
##             2            18             3             6             2 
##   sgrna-msafe  sgrna-ndufaf  sgrna-oligo2   sgrna-rbfox   sgrna-rraga 
##             4             4            13             4             8 
##  sgrna-sh3gl2     sgrna-srf   sgrna-stk39    sgrna-tbk1   sgrna-trem2 
##             4             2            19             2             5

Creates a piechart showing the proportion of cells having 1, 2, 3 etc unique features

In most cases, most cells ideally have only 1 unique barcode

piechart_barcode_uniqueness(obj, gRNA_list)
##   gRNA_count numCells
## 1          1      124
## 2          2        8
## 3          3        1

Bar plot showing the number of cells for each perturbation

barplot_barcode_count(obj, Multiple = TRUE, Neg = FALSE)

Creates 2 barplots: number of DEGs and average fold change of DEGs, for each perturbation

plot_perturbation_effect(obj, perturbations = c("sgrna-oligo2", "sgrna-rraga", "sgrna-fasn", "sgrna-Cfap410", "sgrna-clu", "sgrna-gfap", "sgrna-stk39", "sgrna-trem2"))
##               numDEG   Avg_FC
## sgrna-oligo2      23 3.348977
## sgrna-rraga       51 4.094054
## sgrna-fasn         5 5.160619
## sgrna-Cfap410     28 4.612094
## sgrna-clu          7 4.051291
## sgrna-gfap        46 3.835338
## sgrna-stk39       13 3.365290
## sgrna-trem2      110 5.649040

Dotplot of upregulated genes for each perturbation

barcodes_to_include <- c("sgrna-Cfap410", "sgrna-clu", "sgrna-fasn", "sgrna-flcn", "sgrna-gfap", 
                         "sgrna-msafe", "sgrna-ndufaf", "sgrna-oligo2", "sgrna-rbfox", 
                         "sgrna-rraga", "sgrna-sh3gl2", "sgrna-stk39", "sgrna-trem2")

dotplot_upregulated(obj, barcodes_to_include = barcodes_to_include)
## Calculating cluster sgrna-stk39
## Calculating cluster sgrna-sh3gl2
## Calculating cluster sgrna-msafe
## Calculating cluster sgrna-trem2
## Calculating cluster sgrna-rbfox
## Calculating cluster sgrna-rraga
## Calculating cluster sgrna-clu
## Calculating cluster sgrna-ndufaf
## Calculating cluster sgrna-fasn
## Calculating cluster sgrna-Cfap410
## Calculating cluster sgrna-gfap
## Calculating cluster sgrna-flcn
## Calculating cluster sgrna-oligo2

dotplot_upregulated(obj, barcodes_to_include = barcodes_to_include, include_sgrna = TRUE)
## Calculating cluster sgrna-stk39
## Calculating cluster sgrna-sh3gl2
## Calculating cluster sgrna-msafe
## Calculating cluster sgrna-trem2
## Calculating cluster sgrna-rbfox
## Calculating cluster sgrna-rraga
## Calculating cluster sgrna-clu
## Calculating cluster sgrna-ndufaf
## Calculating cluster sgrna-fasn
## Calculating cluster sgrna-Cfap410
## Calculating cluster sgrna-gfap
## Calculating cluster sgrna-flcn
## Calculating cluster sgrna-oligo2

Scatterplot showing the spatial location of different perturbations in the tissue

cells <- subset(x = obj, subset = x < 19500 & y > 5000 & x > 15000 & y < 8000)
Idents(cells) <- cells$gRNA
plot_barcode_spread(cells)

Given a list of cellIDs and the Seurat object, returns the cellIDs of its nearest neighbours (euclidean distance)

Adapted from: https://cran.rstudio.com/web/packages/RcppHNSW/index.html

cell_list <- c("69243462750838", "82231443856126", "22166326236902", "67598490276766", "81999515620420", "87024627356747", "76965813951447", "32336808792990", "21689584866416")
get_euclidean_neighbours(obj, cell_list = cell_list)
##   [1] "70403103920949" "67740224197433" "69136088569111" "71111773525255"
##   [5] "66249870546059" "65085934408629" "69836168238731" "67748814132850"
##   [9] "73985106645788" "71008694310656" "70579197581101" "69900592748368"
##  [13] "70257075033928" "71897752540863" "70871255357225" "82106889804524"
##  [17] "82029580393236" "82291573398327" "82201379084991" "82154134444867"
##  [21] "82055350196928" "82544976468712" "81930796145445" "82502026795832"
##  [25] "82433307319111" "81823421963028" "81905026341568" "82613695945522"
##  [29] "82205674052457" "82390357645974" "22286585321211" "22016002381543"
##  [33] "21921513100986" "22449794078397" "22376779634337" "22286585321270"
##  [37] "22522808522538" "22664542443238" "22484153816909" "21582210684620"
##  [41] "22741851854610" "22621592770377" "22793391462124" "21569325782683"
##  [45] "22836341135106" "69243462750838" "66189741004463" "68633577396133"
##  [49] "68002217203654" "68947110008756" "68053756811232" "67306432501739"
##  [53] "81960860914947" "80925773796399" "82617990911432" "79998060860335"
##  [57] "81973745816102" "83558588749347" "79929341384122" "85203561222922"
##  [61] "78297253811485" "83554293782423" "85358180046566" "86053964748288"
##  [65] "85869281154682" "77618648978177" "77902116820601" "88085484279420"
##  [69] "89803471197210" "86560770889577" "91057601647897" "89275190220662"
##  [73] "91126321124607" "18013092859954" "21431886828416" "22024592315291"
##  [77] "77068893166570" "76918569311263" "77129022708628" "77283641531398"
##  [81] "77137612643368" "76944339115065" "77283641531286" "77395310681050"
##  [85] "76519137352671" "77163382447163" "76991583755349" "76811195128918"
##  [89] "76394583301101" "77257871727692" "76394583301126" "32216549708686"
##  [93] "32461362844552" "32263794348918" "32422708138857" "32083405722535"
##  [97] "32693291078585" "32723355849691" "31825707684790" "31885837226976"
## [101] "32701881013265" "31898722128901" "31748398273491" "31649614025683"
## [105] "31585189516230" "33114197873641" "21547850945603" "21921513100370"
## [109] "21431886828673" "21689584866490" "21638045258788" "21324512646236"
## [113] "21388937155643" "22067541988433" "21659520095444" "21242908267636"
## [117] "21762599310345" "21517786174466" "21333102580931" "22144851399730"
## [121] "21187073692744"

Given a list of CellIDs, return CellIDs of K nearest neighbours

Adapted from: https://github.com/prabhakarlab/Banksy

cell_list <- c("69243462750838", "82231443856126", "22166326236902", "67598490276766", "81999515620420", "87024627356747", "76965813951447", "32336808792990", "21689584866416")
get_KNN(obj, cell_list)
## Computing neighbors...
## Spatial mode is kNN_r
## Parameters: k_geom=15
## Done
## [1] "Number of unique neighbours: 122"
##   [1] "70403103920949" "67740224197433" "81565723923563" "89008902247666"
##   [5] "82149839477514" "77068893166570" "32345398727559" "21814138917999"
##   [9] "22290880288476" "66361539695565" "81432579937347" "86053964748288"
##  [13] "82291573398292" "77154792512528" "32422708138928" "21522081141861"
##  [17] "22286585321211" "67598490276766" "66314295055388" "81960860914947"
##  [21] "84993107825966" "82334523071214" "77266461662150" "32220844676014"
##  [25] "21547850945603" "22016002381543" "69136088569111" "67160403612942"
##  [29] "80925773796399" "85980950304285" "82240033790689" "76918569311263"
##  [33] "32216549708686" "21921513100370" "22046067152578" "66249870546059"
##  [37] "80406082753475" "85203561222922" "82106889804524" "77129022708628"
##  [41] "32461362844552" "21431886828673" "22144851400371" "71111773525255"
##  [45] "65786014077839" "83421149795616" "89558658061357" "82407537515263"
##  [49] "77283641531398" "32263794348918" "21754009375793" "22222160811699"
##  [53] "69243462750838" "82617990911432" "88085484279420" "82162724379428"
##  [57] "77137612643368" "32422708138857" "21689584866490" "22415434340064"
##  [61] "70497593202037" "80212809225537" "85869281154682" "82029580393236"
##  [65] "76828374998066" "32169305068403" "21638045258788" "22282290353837"
##  [69] "65085934408629" "79998060860335" "89803471197210" "82291573398327"
##  [73] "76944339115065" "32323923891044" "21792664081594" "21921513100986"
##  [77] "81973745816102" "85298050504302" "82471962024730" "77283641531286"
##  [81] "32083405722535" "21324512646236" "22436909176585" "67748814132850"
##  [85] "80315888440776" "88124138985182" "82201379084991" "77395310681050"
##  [89] "32633161536392" "21388937155643" "22449794078397" "69917772617213"
##  [93] "68676527068814" "83558588749347" "85358180046566" "82033875360557"
##  [97] "77038828395579" "32023276180360" "21393232123048" "22376779634337"
## [101] "73199127630572" "79929341384122" "86560770889577" "82343113005761"
## [105] "76519137352671" "32693291078585" "21466246567097" "22286585321270"
## [109] "67757404067540" "83532818945685" "82154134444867" "77279346564130"
## [113] "32723355849691" "21328807613590" "22531398457058" "69836168238731"
## [117] "85392539785072" "82055350196928" "76583561862043" "31860067423172"
## [121] "22067541988433" "21981642643095"

Given a list of CellIDs, return CellIDs of neighbours within a certain radius

Adapted from: https://cran.rstudio.com/web/packages/RcppHNSW/index.html

cell_list <- c("69243462750838", "82231443856126", "22166326236902", "67598490276766", "81999515620420", "87024627356747", "76965813951447", "32336808792990", "21689584866416")
get_radial_neighbours(obj, cell_list)
## $`69243462750838`
## character(0)
## 
## $`82231443856126`
##  [1] "82390357645974" "82368882809458" "82334523071214" "82257213659788"
##  [5] "82240033790556" "82201379084991" "82106889804524" "82055350196928"
##  [9] "81905026341568" "81879256537716" "81776177322702" "81703162878719"
## [13] "81655918238381" "81505594383020" "81393925233403" "82875688950429"
## [17] "82665235552958" "82544976468712" "82433307319111" "82381767711641"
## [21] "82291573398327" "82205674052457" "82154134444867" "82111184771968"
## [25] "82029580393236" "81930796145445" "81823421963028" "81806242093873"
## [29] "81720342747962" "81621558500134" "81561428958021" "81389630266144"
## [33] "82948703394668" "82691005356890" "82613695945522" "82502026795832"
## 
## $`22166326236902`
##  [1] "21968757741125" "21921513100986" "21779779180111" "21582210684620"
##  [5] "21569325782683" "22827751200403" "22793391462124" "22664542443238"
##  [9] "22449794078397" "22376779634337" "22286585321211" "22256520549957"
## [13] "22097606759981" "22016002381543" "22887880742672" "22836341135106"
## [17] "22780506560358" "22741851854610" "22621592770377" "22604412901225"
## [21] "22522808522538" "22484153816909" "22286585321270"
## 
## $`67598490276766`
## [1] "67740224197433"
## 
## $`81999515620420`
## [1] "81565723923563" "81960860914947"
## 
## $`87024627356747`
## character(0)
## 
## $`76965813951447`
##  [1] "76519137352671" "76394583301101" "76291504085940" "76261439314806"
##  [5] "76192719838115" "76124000361453" "77799037606884" "77730318130099"
##  [9] "77661598653427" "77575699307355" "77395310681050" "77283641531286"
## [13] "77129022708628" "77068893166570" "76944339115065" "76918569311263"
## [17] "76811195128918" "76690936044673" "76596446764140" "76527727287356"
## [21] "76459007810669" "76416058137686" "76394583301126" "76373108464699"
## [25] "77661598653503" "77378130811999" "77283641531398" "77257871727692"
## [29] "77163382447163" "77137612643368" "76991583755349"
## 
## $`32336808792990`
##  [1] "31825707684790" "31748398273491" "31649614025683" "31585189516230"
##  [5] "32723355849691" "32693291078585" "32461362844552" "32422708138857"
##  [9] "32263794348918" "32216549708686" "32083405722535" "31885837226976"
## [13] "33114197873641" "32701881013265" "31898722128901"
## 
## $`21689584866416`
##  [1] "21835613754349" "21582210683870" "21380347220939" "22076131922929"
##  [5] "21998822511549" "20882131014719" "21981642642654" "21947282904314"
##  [9] "21921513100370" "21814138917999" "21762599310345" "21689584866490"
## [13] "21659520095444" "21638045258788" "21547850945603" "21517786174466"
## [17] "21444771730680" "21431886828673" "21388937155643" "21333102580931"
## [21] "21324512646236" "21251498202329" "21242908267636" "21187073692744"
## [25] "21118354216184" "21049634739391" "21006685066321" "20920785720512"
## [29] "20920785720428" "22497038718024" "22424024273953" "22393959502932"
## [33] "22144851399730" "22067541988433" "21590800618790" "21483426436389"
## [37] "21294447875351" "21191368660235" "22084721857793"

Frequency plot tabulating the number of positive neighbours each positive cell has

The more the number of neighbours (k_geom), the more likely a cell will have positive neighbours

Idents(obj) <- obj$gRNA_exp
cell_list <- WhichCells(obj, ident = "Pos")

positive_neighbor_histogram(obj, cell_list, k_geom = 10)

positive_neighbor_histogram(obj, cell_list, k_geom = 20)

VlnPlot showing module score of each perturbation

Interrogates the LIANA database to identify LR interactions, visualizating cell-cell communication between 2 groups of cells (source group and target group)

results <- liana_table(obj)

results %>%
  filter(source == "sgrna-srf" | source == "sgrna-tbk1" | source == "sgrna-rraga") %>%
  liana_dotplot(ntop = 20)

To compare at cell-cell communication between Lrrk2-KO cells and mSafe-KO cells and their respective WT neighbours.

This takes into account spatial information.

Idents(obj) <- obj$gRNA
cell_list <- WhichCells(object = obj, ident = "sgrna-lrrk2")
neighbours <- get_KNN(obj, cell_list = cell_list)
## [1] "Number of unique neighbours: 30"
cell_list <- c(cell_list, neighbours)
obj_selected <- subset(obj, cells = cell_list)
results1 <- liana_table(obj_selected)

cell_list <- WhichCells(object = obj, ident = "sgrna-msafe")
neighbours <- get_KNN(obj, cell_list = cell_list)
## [1] "Number of unique neighbours: 60"
cell_list <- c(cell_list, neighbours)
obj_selected <- subset(obj, cells = cell_list)
results2 <- liana_table(obj_selected)

results <- rbind(results1, results2)

results %>%
  arrange(desc(sca.LRscore)) %>%
  filter(source == "sgrna-lrrk2" | source == "sgrna-msafe") %>%
  liana_dotplot(ntop = 20)

Plot of each cell and their neighbours, with each cell coloured by expression level of a specified ligand/receptor pair

cell_list_source = c("70596377451715", "69020124453852") # top and bottom
cell_list_target = get_KNN(obj, cell_list_source)
## [1] "Number of unique neighbours: 30"
plot_ligand_receptor(obj, cell_list_source, cell_list_target, ligand = "Calm1", receptor = "Ptpra")

Estimating power based on a gamma fit model in a scRNA-seq dataset to determine the ability to detect true DE, using scPower.

power_estimation(obj, perturbation = "sgrna-clu", control = "sgrna-msafe", nSamples = 10)
## [1] "The power is estimated to be 0.651017159753292. A power of 0.8 is generally accepted to be sufficient."
power_estimation(obj, perturbation = "sgrna-clu", control = "sgrna-msafe", nSamples = 20)
## [1] "The power is estimated to be 0.861902591224339. A power of 0.8 is generally accepted to be sufficient."
end_time <- Sys.time()
print(paste0("Vignette runtime (min): ", as.numeric(difftime(end_time, start_time, units = "mins"))))
## [1] "Vignette runtime (min): 3.205151450634"