获取multi-point tracking results
一般来说，用ImageJ中的TrackMate插件
就可以实现比较好的Tracking效果
通过设置detector、filter、tracker等一系列参数，然后一路next
到达下面这个步骤时，就可以把每条track上的spots信息导出了
点击Analysis，就会弹出3个statistics results窗口
选择Spots in track statistics窗口
里面就包含了每条track上面，spot随时间的所有位置
点击File/Save As将这个Results保存为csv文件
通过这些数据，就可也计算出MSD（Mean Square Displacement）
MSD是细胞内粒子tracking的一个重要指标
用RStudio处理csv文件
目标非常简单，我希望能够用R计算出每条track对应的MSD(t)
然后利用传说中很好看的ggplot作图
那么代码就得包括这么几个部分
一是文件读取
二是MSD计算和一般统计
三是作图
下面直接上代码
# multi-point auto-tracking MSD analysis and plotting
# 读取文件为 csv 格式
# 2018-5-26
# 设定当前工作路径为数据文件所在路径
filepath <- 'F:/DATA/tracking'
setwd(filepath)
# 读取数据
dat <- read.csv('TrackMate20180526.csv')
# 定义MSD计算公式，参考维基百科
msd <- function(x,y){
n_points <- length(x)
SD <- numeric(length = n_points)
MSD <- numeric(length = n_points)
x1 <- x[1]
y1 <- y[1]
for(i in 1:n_points){
xi <- x[i]
yi <- y[i]
SD <- (xi-x1)^2 + (yi-y1)^2
MSD[i] <- sum(SD)/i
}
return(MSD)
}
# 缩小数据集
tracks <- data.frame(
dat['TRACK_ID'],
dat['POSITION_X'],
dat['POSITION_Y'],
dat['POSITION_T']
)
#按照 TrackID 进行分组
tracks_list <- split(tracks,tracks['TRACK_ID'])
#统计 Track 的总条数
n_tracks <- length(tracks_list)
#构造 data.frame空容器
cdata <- data.frame(TrackID=0,MSD=0,T=0)
cdata <- cdata[-1,]
#批量计算每条 track 的 MSD(t)
for(i in 1:n_tracks){
t <- tracks_list[[i]]
# 从 list 中取值就是要用[[]]
TrackID <- t[['TRACK_ID']]
x <- t[['POSITION_X']]
y <- t[['POSITION_Y']]
T <- t[['POSITION_T']]
MSD <- msd(x,y)
#计算得到 msd 值,存入容器中
c <- data.frame(TrackID,MSD,T)
cdata <- rbind(cdata,c)
}以上代码还只是能完成计算MSD的作用
但是有几个点，想mark一下
第一个就是R的数据类型
list是一个能装任何东西的大杂烩
访问其中内容要用[["item"]]
而dataframe是一种特殊的list
ggplot2要求data得是dataframe类型
第二个就是R的built-in function
比如split(), rbind()
都大大减小了我的工作量
所以一些很常用的功能，
千万不要想着自己去写个function
而是要多查找一下documents
才是最快捷的
用ggplot2包作漂亮的图
## 利用 ggplot 作图
library(ggplot2)
# 所有track的MSD(t)
p1<-ggplot(cdata,aes(T,MSD,color=factor(TrackID)))
p1 + geom_line(alpha=0.2) + guides(color=FALSE)接着上面那段代码，导入ggplot的包
两行简单的命令就可以做出所有track的
MSD随时间变化的曲线图（有点乱）
每条曲线加上透明效果(alpha = 0.2)
可以看到绝大部分粒子处于低速运动状态
这个图花里胡哨不能直观获得关键信息，所以还要进一步处理
接下来咱想看下每个时刻所有例子MSD的均值（含error bar）变化情况
因为是直接从cookbook for R 【链接】里面找的，所以加上一段简单统计的function
叫做summarySE()，可以对dataframe中的数据做统计，得到mean、se、sd、ci
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
## data: a data frame.
## measurevar: the name of a column that contains the variable to be summariezed
## groupvars: a vector containing names of columns that contain grouping variables
## na.rm: a boolean that indicates whether to ignore NA's
## conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
conf.interval=.95, .drop=TRUE) {
library(plyr)
# New version of length which can handle NA's: if na.rm==T, don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else length(x)
}
# This does the summary. For each group's data frame, return a vector with
# N, mean, and sd
datac <- ddply(data, groupvars, .drop=.drop,
.fun = function(xx, col) {
c(N = length2(xx[[col]], na.rm=na.rm),
mean = mean (xx[[col]], na.rm=na.rm),
sd = sd (xx[[col]], na.rm=na.rm)
)
},
measurevar
)
# Rename the "mean" column
datac <- rename(datac, c("mean" = measurevar))
datac$se <- datac$sd / sqrt(datac$N) # Calculate standard error of the mean
# Confidence interval multiplier for standard error
# Calculate t-statistic for confidence interval:
# e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
ciMult <- qt(conf.interval/2 + .5, datac$N-1)
datac$ci <- datac$se * ciMult
return(datac)
}
对数据进行简单统计处理之后，再用ggplot进行绘制
ccdata<-summarySE(cdata,measurevar = "MSD",groupvars = c("T"))
ggplot(ccdata,aes(T,MSD,colour=N,width=1))+
geom_errorbar(aes(ymin=MSD-se,ymax=MSD+se),width=.1) +
geom_point()+geom_line()+
xlab("Lag Time (s)")+ylab("MSD (μm²)")+
scale_x_log10()+ylim(0,3)+
xlim(0,20)效果如下
总体感觉的话
要想用好ggplot2和R
关键还是要学好统计学
脚本语言本身还是相当简单的