dtw_example.R

# An example to run DTW time delay analysis across modalities.

library(tidyverse)
library(dtw)
library(zoo)

gene = 'EGR1'
bins = 20

# input file contains three columns: latent time, motif accessibility, and TF gene expression
res = read.delim(paste0(gene, '_res.txt'), sep=' ', header = F, stringsAsFactors = F)
colnames(res) = c('t', 'm', 'r')
res = res[order(res$t),]

res$t = res$t * bins
res$window = floor(res$t / 1)
res$window[res$window==bins] = bins-1
res2 = aggregate(res, by=list(res$window), FUN=mean)

rv = res2$r - min(res2$r)
rv = c(rv[1], rv, rv[length(rv)])
rv = rollmean(rv, 3)
rv[1] = 0
rv[length(rv)] = 0
rv = rv / max(rv)
mv = res2$m - min(res2$m)
mv = c(mv[1], mv, mv[length(mv)])
mv = rollmean(mv, 3)
mv[1] = 0
mv[length(mv)] = 0
mv = mv / max(mv)

# compute dtw
dtw.res = dtw(rv, mv, keep.internals = T)
dtwPlotTwoWay(dtw.res, col=c('magenta', '#f7ba41'), lwd=3, xlab='latent time', ylab='accessibility and expression')