R语言之于流式结果分析

需要用到的R包有:flowCore,flowStats

这里以Guava Mini2, InCyte 3.1的数据为例

library(flowCore)  
library(ggplot2)  #我比较习惯于用ggplot2作图
library(MASS)   #用 fitdistr 来进行正态分布曲线拟合
library(plyr)

# 可以先设一下locale，不然会报错，当然报错估计也影响不大
Sys.setlocale(category = "LC_ALL", locale = "English_United States.1252")

dat <- read.FCS("filename", alter.names=TRUE)   # filename改成相应的文件路径, alter.names会把-等改成点.

str(dat) # 看看数据结构，

head(dat@exprs)      # 看一下数据有几列，列名叫什么

dat@parameters@varMetadata

dat@parameters@data        # 看一下exprs各列都是什么含义

dat@parameters@dimLabels

# 为了进一步分析，我们先把数据转成数据框
Dat <- as.data.frame(dat@exprs)

# plot(Dat)  # 可以用plot(Dat)看看各列之间的关系，当然，这个图太大，不好看 

plot(Dat[ , grepl("FSC|SSC", names(Dat))], col="#00000022")   # 看看FSC和SSC的关系

很显然能够看到*.HLin和*.HLog之间似乎没有什么本质的区别，因此直接用*.HLin即可。而且也能看到信号分为两群呢，左下角的一群应该是细胞碎片。

能够看到，其它通道中RED.B.*的信号是最强的，然后我们看看

plot(Dat[ , grepl("RED.B", names(Dat))], col="#00000022")   # 看看FSC和SSC的关系

同样可以看出来，*.Lin和*.Log是一致的。 接下来用*.Lin分析就好了。这样我们主要使用的FSC.HLin、SSC.HLin、RED.B.ALin、RED.B.HLin就好了

DAT <- Dat[ , c("FSC.HLin", "SSC.HLin", "RED.B.ALin", "RED.B.HLin", "RED.B.W")]
names(DAT) <- c("FSC", "SSC", "PIA", "PIH", "PIW")  # 改一下名称以便后边好用, 假装是用PI染色的

接下来就要进行数据的过滤，先通过FSC对大小进行过滤，太小的可能是碎片，太大的可能是细胞团。可能还有别的办法，我这里提供一种我认为可行的。

d <- dist(DAT[ , "FSC"], method="euclidean")   # FSC代表细胞的大小，通过聚类区分
h <- hclust(d, method="ward.D2")
C <- cutree(h, k=3)                          # 把细胞分为三群
DAT <- data.frame(DAT[ , 1:5], Cluster = as.character(C))

ggplot(DAT)+geom_point(aes(x=FSC, y=SSC, colour=Cluster), alpha=0.25) +theme_bw()

看起来还可以，我们需要保留的是第2群细胞。然后用PIA来看细胞周期。 A是area的意思，代表DNA的总量。 先看看DNA总量的分布情况。
应该还要有FSC-A和FSC-W的值，通过FSC-A和H 或者FSC-A和W来对细胞的二聚体或者多聚体进行排除效果估计会更好。

DAT <- DAT[DAT$Cluster == "2", ]

ggplot(DAT) + 
    geom_histogram(aes(x=PIA)， binwidth=5) +
    scale_x_continuous(breaks=seq(0, 1000, 50))+
    theme_bw()

看起来在做流式的时候没有设置得太好，一般来说，在做流失的时候我们把G1期那个峰对应的信号值调到200。

再看看PIA和PIH之间的关系，PIH小而PIA大的可能不是两个或者多个细胞黏在一起了。

ggplot(DAT) + 
    geom_point(aes(x = PIW, y = PIA)， alpha=0.25) +
    scale_x_continuous(breaks=seq(0, 1000, 50))+
    scale_y_continuous(breaks=seq(0, 1000, 50))+
    theme_bw()

看看PIA和PIH、PIW的关系，H代表信号的最大值，W代表信号的宽度。 想象一下一个球形的细胞通过一束光，刚进入的时候开始有信号，球心进入的时候就达到了最大值，这个值就是H，整个过程信号强度就能画出一个曲线，有宽度W，而曲线的下边的面积就是A。

有少部分的细胞偏离，不过界限看起来也不是很明显，根据自己的实际情况决定是否要把偏离的gate掉吧。 一般来说PIW的变化范围应该要比PIA的小，G2/M期的细胞PIW不会相比于G0/G1期的细胞成倍地增加。 但是本例的这个图很难说清楚为什么是这个样子的。

从上边的histogram能够看出各种DNA含量的细胞的分布。最高的峰对应G1期细胞群体的中点。 我们先来预测一下这个中点。

DATA <- ddply(DAT[DAT$PIA >= 200 & DAT$PIA <= 350, ], .(PIA), summarize, Count=length(PIA))
DATA <- merge(data.frame(PIA=min(DATA$PIA):max(DATA$PIA)), DATA, by="PIA", all.x=TRUE)  #有的PIA对应没有细胞
DATA[is.na(DATA)] <- 0   # 把PIA对应没有细胞的细胞数设为0

Fit <- loess(Count~PIA, data=DATA, span=0.75)    # 调节这个span，使得平滑效果合适

DATA$sCount <- predict(Fit, data.frame(PIA=DATA$PIA), se = TRUE)$fit

head(DATA[order(-DATA$sCount), ])

DATA[order(-DATA$sCount), ][1, 1]  # 这里得到的值是 [1] 278


ggplot() +
     geom_histogram(data=DAT, aes(x=PIA), binwidth=1) +
     geom_line(data=DATA, aes(x=PIA, y=sCount), colour="red") +
     scale_y_continuous(limits=c(0, 50), expand=c(0,0))+
     theme_bw()

可以看看拟合的情况，还是比较好的，因此选取278作为G1 DNA含量的均值。 这个时候，我们假设G1细胞的PIA值呈正态分布，那么均值右边的细胞数和左边的细胞数应该是一致的。 而G2期细胞的信号的均值可以看做是G1的两倍。 由于这里G2/M期的峰不明显，没法知道实际的值，也只能按照两倍G1期的来算了。

# G1的细胞数：
nG1 <- 2*nrow(DAT[DAT$PIA < 278 & DAT$PIA > 150, ]) + nrow(DAT[DAT$PIA == 278, ]) #200 这里是比较武断的选择的。

nG2 <- 2*nrow(DAT[DAT$PIA > 556 & DAT$PIA < 750, ]) + nrow(DAT[DAT$PIA == 556, ]) #750 这里是比较武断的选择的。

nS <- nrow(DAT[DAT$PIA > 150 & DAT$PIA < 750, ])  - nG1 - nG2

这样我们得到G1、S、G2/M的细胞数依次为2128、1101、1640，百分比就依次为43.71%、22.61%、33.68%。很显然，本例中还有很多PIA比较大的，有可能是多倍体细胞。仅从FSC-H进行了过滤，不确定是不是把细胞二聚体或者多聚体都去掉了。

<html>

<div style=“font-size:100%;color:red;”> <b>关于G1、S、G2/M各时期的细胞数的计算只是提出一种办法，其可靠性还需要读者自己判断。 特别是G2期的计算，像本例这样峰不明显的情况这么处理是否合适还有待商榷。 如果G2的峰也比较明显，可以采用G1类似的办法。 另外，上边这个例子还有很多多倍体或者细胞团的干扰，因此G2估计可能更加不可靠了，然后S期的比例相应地也不那么可靠了。 </b> </div>

</html>

做流式的时候，把G1期的峰都调到200，这样就便于不同样本之间一样的程序分析了，不然就需要先归一化一下数据。 可以自己研究一下flowStats包。

flowStats看起来不错。待续。

---

自2019年9月30日以来，本网页已被访问9184次。