Temporal Cox model
Hening Cui
5/9/2023
library(Seurat)
library(spacexr)
library(tidyverse)
library(spatstat)
library(nlme)
library(lgcp)
library(CAinterprTools)loading dataset
load(file= "/Users/cuihening/Desktop/shikun/seurat/robjs/chicken_visium.4.Robj")TM <- FetchData(object=chicken_visium,vars=c('TMSB4X','ident', 'orig.ident'))
tmsb4x=cbind(cell = rownames(TM), TM) %>%
dplyr::rename(time = orig.ident)
tmsb4x$is2=as.factor(ifelse(tmsb4x$ident==2,1,0))The function to get coordinate
getcoor=function(slice, time){
coords <- GetTissueCoordinates(chicken_visium, image= slice)
colnames(coords) <- c("x", "y")
coords <- cbind(cell = rownames(coords),coords)
coords$time = time
return(coords)
}coords_d4 = getcoor('slice1','D4')
coords_d10 = getcoor('slice1_D10-C1','D10')
coords_d7 = getcoor('slice1_D7-B1','D7')
coords_d14 = getcoor('slice1_D14-D1','D14')Function for constuct ppp
buildppp=function(dataset,markers){
X <- dataset$x
Y <- dataset$y
xrange <- range(X, na.rm=T)
print(xrange)
yrange <- range(Y, na.rm=T)
print(yrange)
ppppro = ppp(X,Y,xrange,yrange,marks=markers)
return(ppppro)
}EDA
D4
D4 = tmsb4x %>%
filter(time=='D4') %>%
merge(coords_d4, by="cell") %>%
filter(x >= 385 & y>=267 ) %>%
rename(Y=x, X=y) %>%
rename(x=X, y=Y) %>%
mutate(x=x-138, y=y-208) %>%
mutate(time=4)
summary(D4)## cell TMSB4X ident time.x is2
## Length:141 Min. :0.04300 6 :66 Length:141 0:119
## Class :character 1st Qu.:0.05200 4 :27 Class :character 1: 22
## Mode :character Median :0.05554 2 :22 Mode :character
## Mean :0.05688 8 :21
## 3rd Qu.:0.05953 5 : 5
## Max. :0.07820 0 : 0
## (Other): 0
## y x time.y time
## Min. :182.6 Min. :142.5 Length:141 Min. :4
## 1st Qu.:202.2 1st Qu.:180.0 Class :character 1st Qu.:4
## Median :221.8 Median :206.3 Mode :character Median :4
## Mean :222.8 Mean :201.6 Mean :4
## 3rd Qu.:241.3 3rd Qu.:225.1 3rd Qu.:4
## Max. :274.0 Max. :247.5 Max. :4
## D4_samppp = ppp(D4$x,D4$y, c(0, 450),c(0,550), marks=D4$is2)D4_ty2=D4 %>%
filter(is2==1)
D4_ty2p=ppp(D4_ty2$x,D4_ty2$y, c(0, 450),c(0,550))par(mfrow=c(1,2))
plot(D4_samppp, pch = 20,show.window=FALSE,cols=c("blue", "green"), legend=FALSE, main="D4 point",markscale = 0.1)
plot(density(D4_ty2p), main="D4 density")
D7
D7 = tmsb4x %>%
filter(time=='D7') %>%
merge(coords_d7, by="cell") %>%
filter(x >= 255 & y>=267 ) %>%
rename(Y=x, X=y) %>%
rename(x=X, y=Y) %>%
mutate(x=x-150, y=y-152) %>%
mutate(time=7)
summary(D7)## cell TMSB4X ident time.x is2
## Length:505 Min. :0.00000 1 :166 Length:505 0:438
## Class :character 1st Qu.:0.05050 5 : 83 Class :character 1: 67
## Mode :character Median :0.05403 2 : 67 Mode :character
## Mean :0.05413 4 : 53
## 3rd Qu.:0.05735 8 : 48
## Max. :0.07924 7 : 35
## (Other): 53
## y x time.y time
## Min. :107.8 Min. :131.1 Length:505 Min. :7
## 1st Qu.:179.5 1st Qu.:172.2 Class :character 1st Qu.:7
## Median :225.1 Median :198.6 Mode :character Median :7
## Mean :220.9 Mean :200.3 Mean :7
## 3rd Qu.:264.3 3rd Qu.:228.6 3rd Qu.:7
## Max. :323.1 Max. :281.1 Max. :7
## D7_samppp = ppp(D7$x,D7$y, c(0, 450),c(0,550), marks=D7$is2)D7_ty2=D7 %>%
filter(is2==1)
D7_ty2p=ppp(D7_ty2$x,D7_ty2$y, c(0, 450),c(0,550))par(mfrow=c(1,2))
plot(D7_samppp, pch = 20,show.window=FALSE, cols=c("blue", "green"), legend=FALSE, main="D7 point")
plot(density(D7_ty2p), main="D7 density")
D10
D10 = tmsb4x %>%
filter(time=='D10') %>%
merge(coords_d10, by="cell") %>%
filter(x >= 193 & y>=320) %>%
rename(Y=x, X=y) %>%
rename(x=X, y=Y) %>%
mutate(x=x-217, y=y-128) %>%
mutate(time=10)
summary(D10)## cell TMSB4X ident time.x is2
## Length:888 Min. :0.03259 3 :163 Length:888 0:760
## Class :character 1st Qu.:0.05015 2 :128 Class :character 1:128
## Mode :character Median :0.05414 1 :123 Mode :character
## Mean :0.05418 4 :118
## 3rd Qu.:0.05806 0 :109
## Max. :0.07277 7 :108
## (Other):139
## y x time.y time
## Min. : 72.12 Min. :103.4 Length:888 Min. :10
## 1st Qu.:169.94 1st Qu.:156.0 Class :character 1st Qu.:10
## Median :215.79 Median :197.4 Mode :character Median :10
## Mean :220.37 Mean :200.8 Mean :10
## 3rd Qu.:267.97 3rd Qu.:238.5 3rd Qu.:10
## Max. :378.87 Max. :321.0 Max. :10
## D10_samppp =ppp(D10$x,D10$y, c(0, 450),c(0,550), marks=D10$is2)D10_ty2=D10 %>%
filter(is2==1)
D10_ty2p=ppp(D10_ty2$x,D10_ty2$y, c(0, 450),c(0,550))par(mfrow=c(1,2))
plot(D10_samppp, pch = 20,show.window=FALSE, cols=c("blue", "green"), legend=FALSE, main="D10 point")
plot(density(D10_ty2p), main="D10 density")
D14 = tmsb4x %>%
filter(time=='D14') %>%
merge(coords_d14, by="cell") %>%
rename(Y=x, X=y) %>%
rename(x=X, y=Y) %>%
mutate(x=x-144, y=y-78) %>%
mutate(time=14)
summary(D14)## cell TMSB4X ident time.x is2
## Length:1967 Min. :0.03193 0 :854 Length:1967 0:1643
## Class :character 1st Qu.:0.05184 2 :324 Class :character 1: 324
## Mode :character Median :0.05530 3 :175 Mode :character
## Mean :0.05509 4 :159
## 3rd Qu.:0.05871 8 :138
## Max. :0.07352 7 :128
## (Other):189
## y x time.y time
## Min. : 9.428 Min. : 58.58 Length:1967 Min. :14
## 1st Qu.:139.719 1st Qu.:137.43 Class :character 1st Qu.:14
## Median :224.323 Median :200.63 Mode :character Median :14
## Mean :220.551 Mean :200.92 Mean :14
## 3rd Qu.:302.755 3rd Qu.:261.00 3rd Qu.:14
## Max. :400.476 Max. :361.98 Max. :14
## D14_samppp = ppp(D14$x,D14$y, c(0, 450),c(0,550), marks=D14$is2)D14_ty2=D14 %>%
filter(is2==1)
D14_ty2p=ppp(D14_ty2$x,D14_ty2$y, c(0, 450),c(0,550))par(mfrow=c(1,2))
plot(D14_samppp, pch = 20,show.window=FALSE, cols=c("blue", "green"), legend=FALSE, main="D14 point")
plot(density(D14_ty2p), main="D14 density")
Cox Spatial Temporal Model
prepare data
temp= rbind(D4, D7, D10, D14) %>%
filter(ident==2) %>%
select(x,y,time)tlim=c(0,14)xyt <- stppp(list(data = temp, tlim = tlim, window = owin(c(0,450),c(0,550))))xyt <- integerise(xyt)
xyt## Space-time point pattern
## Planar point pattern: 541 points
## window: rectangle = [0, 450] x [0, 550] units
## Time Window : [ 0 , 14 ]den <- lambdaEst(xyt, axes = TRUE)
plot(den)sar <- spatialAtRisk(den)
sar## SpatialAtRisk object
## X range : [1.7578125,448.2421875]
## Y range : [2.1484375,547.8515625]
## dim : 128 x 128plot(sar)
mut <- muEst(xyt)
mut## temporalAtRisk object
## function (t)
## {
## if (length(t) > 1) {
## stop("Function only works for scalar numeric objects, ie a vector of length 1.")
## }
## if (!any(as.integer(t) == tvec)) {
## return(NA)
## }
## return(obj[which(as.integer(t) == tvec)] * scale)
## }
## <bytecode: 0x7f95652b9cc8>
## <environment: 0x7f95652b8ad8>
## attr(,"tlim")
## [1] 0 14
## attr(,"class")
## [1] "temporalAtRisk" "function"
## Time Window : [ 0 , 14 ]plot(mut)
gin <- ginhomAverage(xyt, spatial.intensity = sar, temporal.intensity = mut)##
|
| | 0%
|
|===================== | 30%
|
|========================================== | 60%
|
|======================================================================| 100%
## Returning an average of 4 curvesplot(gin)
kin <- KinhomAverage(xyt, spatial.intensity = sar, temporal.intensity = mut)##
|
| | 0%
|
|===================== | 30%
|
|========================================== | 60%
|
|======================================================================| 100%
## Returning an average of 4 curves