R語言實現流式細胞數據分析

  • 2020 年 2 月 11 日
  • 筆記

流式細胞術通過光學檢測系統快速檢測多參數的細胞流。許多因素使得流式細胞術能夠成功和廣泛的應用，比如檢測速度（能夠允許大量的細胞被檢測），高度的準確性和解析度，低成本。此外，流式細胞術還是一種非破壞性技術，可以分選出活細胞用於後續分析。能夠分析和分選單個細胞的能力使流式細胞術在生物學和醫學領域有非常廣泛的應用。在免疫學中，流式細胞術用來鑒定和量化免疫細胞亞群，因此可以監控病人的免疫狀態，通過比較不同的病人組也可以找出生物標誌物。具體的原理：一定波長的雷射束直接照射到高壓驅動的液流，產生的光訊號被多個接收器接受，一個是機關束直線方向上接受的前向角散射光訊號。其他是在雷射束垂直方向上接受的光訊號，包括側向角散射光訊號和熒光訊號，這些光訊號被相應的接受器接受後，根據接收訊號的強弱就能反應出細胞的物理和化學特徵。

其中主要的訊號包括：前向角散射訊號（FSC）反映細胞的大小；側向角散射訊號（SSC）反映細胞顆粒度值越大，細胞內顆粒結構越複雜，品質越大；熒光訊號（FL）反映胞內、胞外染色情況。接下來我們看下在R語言如何進行整個的分析流程。首先我們需要安裝其幾個重要的包：

##數據源：  BiocManager::install("flowWorkspaceData")    ##數據讀取和質控：  BiocManager::install("flowCore")  BiocManager::install("flowClean")    ##批量標準化：  BiocManager::install("flowStats")    ##自動圈門包：  BiocManager::install("openCyto")  BiocManager::install("flowDensity")    ##結果的可視化：  BiocManager::install("ggcyto")  BiocManager::install("flowViz")

以上的包便可以完成流式的整個分析過程，接下來我們就看下如何實現整個分析流程：

1. 載入數據

library(flowWorkspaceData)  dataDir <-system.file("extdata",package="flowWorkspaceData")    library(flowCore)  file.name <-system.file("extdata","0877408774.B08",package="flowCore")  fcs.dir <-system.file("extdata","compdata","data",package="flowCore")

2. 讀取數據及可視化

首先我們看下flowCore包中的讀取數據的函數read.FCS和read.flowSet：

其中主要的參數：

Transformation 主要是指的原始數據的轉化，默認是線性轉化，一般我不需要轉化直接設為FALSE

Which.lines 主要是指的讀取的範圍，也就是行數eg：which.lines=2:5

實例：

X=read.FCS(file.name, transformation=FALSE)  Y=read.FCS(paste(dataDir,'/CytoTrol_CytoTrol_1.fcs',sep=''),transformation=FALSE)  exprs(X[1:10,])##獲取對應的訊號的值  parameters(X)$name###獲取X中的變數  keyword(X) ###獲取文件的關鍵資訊
library(ggcyto)  autoplot(X, "FL1-H","FL2-H")
fs <- read.flowSet(path = fcs.dir)    sampleNames(fs)##查看樣本編號資訊，多批實驗時可用。  exprs(fs[[1]][1:12,])##獲取數據  autoplot(fs, "FL1-H","FL2-H")

那如果我們的熒光素1對應的通道可以檢測其熒光訊號，那麼其它的通道也會產生熒光素1的訊號，那這些訊號就成為溢出，組成的熒光訊號值的矩陣也就是溢出矩陣，在flowCore中作者提到了對溢出矩陣的處理方法：

首先獲取溢出矩陣：

fcsfiles <- list.files(pattern ="CytoTrol",system.file("extdata",package="flowWorkspaceData"),full = TRUE)  fs <- read.flowSet(fcsfiles)  x <- fs[[1]]  comp_list <- spillover(x)  comp_list

接下來我們看下如何進行標準化處理： 

comp <- comp_list[[1]]  x_comp <- compensate(x, comp)##單批數據處理  comp <- fsApply(fs, function(x)spillover(x)[[1]], simplify=FALSE) ###多批處理  fs_comp <- compensate(fs, comp)
library(gridExtra)  transList <- estimateLogicle(x,c("V450-A","V545-A"))  p1 <- autoplot(transform(x, transList),"V450-A", "V545-A") +ggtitle("Before")  p2 <- autoplot(transform(x_comp,transList), "V450-A", "V545-A")+ggtitle("Before")  grid.arrange(as.ggplot(p1), as.ggplot(p2),ncol = 2)

從上圖我們可以看出在經過溢出矩陣處理後的效果，可以更好的將數據的差異分開。

3. 品質控制

我們主要是看下flowClean中clean函數進行對flowdata類型的數據進行品質控制並獲得質控之後的數據，當然這個函數有一個缺點，那就是細胞數必須大於30000，不然沒法進行質控：

主要的參數：

vectMarkers 指的要進行標準化的列。e.g. 'FSC-A'

nCellCutoff 指一個群體最小的細胞數。

Cutoff 指的閾值，可以是0-1的閾值，當然也可以是大於1的熒光值。

filePrefixWithDir 輸出的文件路徑或者文件名

Ext 文件擴展名

實例：

library(flowClean)  library(flowViz)  data(synPerturbed)  synPerturbed.c <- clean(synPerturbed,vectMarkers=c(5:16),filePrefixWithDir="sample_out",ext="fcs", diagnostic=TRUE)    library(grid)  library(gridExtra)  lgcl <- estimateLogicle(synPerturbed.c,unname(parameters(synPerturbed.c)$name[5:16])) ##獲取要進行轉化的列對應的名稱。  synPerturbed.cl <-transform(synPerturbed.c, lgcl) ##轉化後的數據集  p1 <- xyplot(`<V705-A>` ~ `Time`,data=synPerturbed.cl, abs=TRUE, smooth=FALSE, alpha=0.5, xlim=c(0, 100))  p2 <- xyplot(`GoodVsBad` ~ `Time`,data=synPerturbed.cl, abs=TRUE, smooth=FALSE, alpha=0.5, xlim=c(0, 100),ylim=c(0, 20000))    rg <-rectangleGate(filterId="gvb", list("GoodVsBad"=c(0, 9999)))  idx <- filter(synPerturbed.cl, rg) ###獲取過濾後的結果  synPerturbed.clean <-Subset(synPerturbed.cl, idx)  p3 <- xyplot(`<V705-A>` ~ `Time`,data=synPerturbed.clean, abs=TRUE, smooth=FALSE, alpha=0.5, xlim=c(0, 100))  grid.arrange(p1, p2, p3, ncol=3)

從上圖我們可以看出在過濾後，數據隨著時間的增加會清晰的分開。

我們看下flowCore包中的數據的標準化需要用到函數transfom，我們直接看一個對數轉換的實例：

fs <-read.flowSet(path=system.file("extdata", "compdata","data", package="flowCore"), name.keyword="SAMPLEID", phenoData=list(name="SAMPLE ID",Filename="$FIL"))  autoplot(transform(fs[[1]],`FL1-H`=log(`FL1-H`), `FL2-H`=log(`FL2-H`) ),"FL1-H","FL2-H")

當然我們的標準化形式不止這一種，還有下面的幾種：

那這些的使用就比較複雜了，我們直接上實例：

aTrans <-truncateTransform("truncate at 1", a=1) ##構造參數  myTrans <- transformList('FL1-H',aTrans)  transform(fs, myTrans)

對於一個數據集基本的流程如上，那麼如果多個數據集我們需要用到flowStats，此包主要用於數據的標準化和自動圈門。那麼我們看下具體的應用：

library(flowStats)  fcs.dir <-system.file("extdata", "compdata","data",package="flowCore")  frames <- lapply(dir(fcs.dir,full.names=TRUE), read.FCS)  names(frames) <-c("UNSTAINED", "FL1-H", "FL2-H","FL4-H", "FL3-H")  frames <- as(frames,"flowSet")##構建flowset，當然可以直接用flowCore讀取生成flowset  comp <- spillover(frames,unstained="UNSTAINED", patt = "-H",fsc = "FSC-H",ssc = "SSC-H",stain_match = "ordered")

當然這裡是不考慮通道的對應性時stain_match為ordered，當需要一一對應時則需要另一個參數regexpr.

comp <- spillover(frames,unstained="UNSTAINED", patt = "-H", fsc ="FSC-H", ssc = "SSC-H", stain_match = "regexpr")

那麼除了上面的方法，還提供了spillover_match函數，可以直接對目錄下的文件進行匹配：

comp_match <-system.file("extdata", "compdata", "comp_match",package="flowCore")  control_path <-system.file("extdata", "compdata", "data",package="flowCore")    matched_fs <-spillover_match(path=control_path, fsc = "FSC-H", ssc ="SSC-H", matchfile = comp_match)  comp <- spillover(matched_fs, fsc ="FSC-H", ssc = "SSC-H", prematched = TRUE)

3. 圈門以及可視化

首先我們看下通過演算法的自動圈門過程。我們看下flowCore中的幾種方法：

實例：

rectGate <-rectangleGate(filterId="Fluorescence Region", "FL1-H"=c(0,12), "FL2-H"=c(0, 12))  result = filter(fs[[1]],rectGate) ###獲取此區域的數據  summary(result)  ##Fluorescence Region+: 9811 of 10000events (98.11%)

當然多個數據集也是可以的：

result = filter(fs,rectGate) ###獲取此區域的數據  summary(result)

最後還能針對不同的樣本獲取不同的門：

rg2 <- rectangleGate(filterId="nonDebris","FSC-H"=c(300,Inf))  rg3 <-rectangleGate(filterId="nonDebris","FSC-H"=c(400,Inf))  flist <- list(rectGate, rg2, rg3)  names(flist) <- sampleNames(foo)  fr3 <- filter(foo, flist)

接下來我們看下在flowStats是如何實現這個過程的，其中主要的函數是norm2Filter，我們直接看下實例：

library(flowStats)  morphGate <- norm2Filter("FSC-H","SSC-H", filterId="MorphologyGate", scale=2)  smaller <- Subset(fs, morphGate)

接下來我們看下如何進一步對門數據進行k-mean聚類分析：

split(smaller[[1]],kmeansFilter("FSC-H"=c("Low","Medium","High"),filterId="myKMeans"))

當然也可以直接對所有的樣本進行直接的分組：

split(smaller,kmeansFilter("FSC-H"=c("Low", "Medium","High"), filterId="myKMeans"))

那麼如何自動獲取多個門的數據並進行可視化，需要用到下面的函數，我們直接用一個實例來說明：

fcsfiles <- list.files(pattern ="CytoTrol", system.file("extdata", package ="flowWorkspaceData"), full = TRUE)  fs <- read.flowSet(fcsfiles)    gs <- GatingSet(fs)#構造門數據集    fs_comp <- getData(gs)  autoplot(fs_comp[[1]], "B710-A")+ ggtitle("raw")

我們看到參數有-A和-H等具體是什麼意思。其實早有定義A是area的縮寫，H是height的縮寫，W是width的縮寫。那麼有時候需要增加一個門的話，其實此包還增加了這一部分功能：

rg1 <-rectangleGate("FSC-A"=c(50000, Inf), filterId="NonDebris")  gs_pop_add(gs, rg1, parent ="root")  gs_get_pop_paths(gs) ##獲取當下的節點  recompute(gs)  autoplot(gs, "NonDebris")
ggcyto(gs, aes(x = `FSC-A`)) +geom_density() + geom_gate("NonDebris")

那如何獲取每一部分的百分比呢，我們可以直接看下實例：

dataDir <-system.file("extdata",package="flowWorkspaceData")  suppressMessages(gs <-load_gs(list.files(dataDir, pattern = "gs_manual",full = TRUE)))    # get stats all nodes  dt <- gs_pop_get_stats(gs) #default is"count"    nodes <- c("CD4","CD8")  gs_pop_get_stats(gs, nodes,"percent")
plot(gs)##可視化所有門的組成
autoplot(gs[[1]])#可視化所有的門

當然還可以更加細化門的分類，以及門的構成，在此我們就不再深究了。還有兩個專門做自動圈門的包，我們就不細化講解了當然如果我們想讓可視乎更加的漂亮，也可以利用一些降維的方法操作原始數據，最終實現圈門。

歡迎大家學習交流！

VirMach 便宜 VPS

QNews

QNews