EBImage是一个用于图形处理的R包,简洁优雅功能强大,可以完成很多计算机图形处理算法。
本文使用EBImage完成对一组细胞荧光图的定量分析,数据使用EBImage内置的测试图片。
先看一下EBImage中对图形的定义,EBImage使用readImage
函数读入R,读进来的对象是EBImage包定义的Image对象。
nuc = readImage(system.file('images', 'nuclei.tif', package='EBImage'))
cel = readImage(system.file('images', 'cells.tif', package='EBImage'))
其实这个Image对象就是基于一个多维数组,所以它也可以使用很多数组的函数。它的前两维度固定是图形的宽与高,后面的维度是图层:或者是R、G、B通道的颜色图层,或者是叠加的其他图片(如果是叠加彩色图,那么第三维是通道,第四维是叠加)。
这两种图层的区别在于是否会渲染为一张图片,比如说一张彩色图片是有三个颜色通道的,所以一张彩色图片的维度就是(X,Y,3),渲染是一张图片。如果是三个灰度图叠加在一起,那么其维度也是(X,Y,3),但是渲染的时候,它是三张图片。
在Image中是通过frames.render这个参数来区别的,有几张图片就有多少render图层,所以一张彩色图片的frames.render是1,但是三个灰度图叠加的frames.render是3。
nuc
#Image
# colorMode : Grayscale
# storage.mode : double
# dim : 510 510 4
# frames.total : 4
# frames.render: 4
#
#imageData(object)[1:5,1:6,1]
# [,1] [,2] [,3] [,4] [,5] [,6]
#[1,] 0.06274510 0.07450980 0.07058824 0.08235294 0.10588235 0.09803922
#[2,] 0.06274510 0.05882353 0.07843137 0.09019608 0.09019608 0.10588235
#[3,] 0.06666667 0.06666667 0.08235294 0.07843137 0.09411765 0.09411765
#[4,] 0.06666667 0.06666667 0.07058824 0.08627451 0.08627451 0.09803922
#[5,] 0.05882353 0.06666667 0.07058824 0.08235294 0.09411765 0.10588235
dim(nuc)
#[1] 510 510 4
通过一个例子来看一下:
red <- blue <- green <- matrix(0, nrow = 3, ncol = 3)
red[, 1] <- 1
blue[, 2] <- 1
green[, 3] <- 1
# 一张彩图
p_col <- rgbImage(red = red, blue = blue, green = green)
# 三个单一图形
p_combn <- combine(as.Image(red), as.Image(blue), as.Image(green))
display(p_col, method = "raster", all = TRUE)
display(p_combn, method = "raster", interpolate = FALSE, all = TRUE)
彩色图只渲染一张图,但是三个组合在一起的图层是渲染了三个灰度图,如下图所示。
在R中,颜色使用0-1的数值范围,另外matrix在转换为Image对象时,行数会转换为宽度,列数为转换为高度。 图2中为了区分三个图的边界,人工添加了两条边界线。
先将测试图片合并成一张图形:细胞核和细胞质的扫描图合并成一张彩色图。
从刚才nuc的输出信息可以知道,它是由四个图层组成,这里合并成彩色图时也是四张。
cells <- rgbImage(green=1.5*cel, blue=nuc) # cel是胞质, 乘以1.5代表提高了一下对比度
display(cells, method = "raster", all = TRUE)
Thresholding阈值化
在做图形处理的时候,一般都是现将图片转换成灰度图进行处理,再进行阈值化。
上面的显微镜扫描的原始图其实已经是灰度图了,需要再进行阈值化。
可以有两个策略:一个是直接寻找一个阈值,原图减去即可。另一个是本地化策略,比较前景和背景的差异决定其阈值化的策略。
原图减去阈值
# 只以第一个图层为例
nuc %>% getFrame(1) %>% EBImage::otsu()
#[1] 0.3535156
# 减去阈值
# combine函数可以将读个图层合并成一个Image对象
combine(nuc %>% getFrame(1) %>% {. > 0.35},
nuc %>% getFrame(1)
) %>%
display("raster", all = TRUE)
本地化阈值
为了获取一个图片的各个位置的背景表达,其实是通过一个滤镜实现的,然后再定义一个cutoff值代表超过背景值多少既可以认为有表达。
由于这个过程非常常用,所以EBImage封装了一个函数来执行这个过程:thresh函数。
一般做完阈值后,还会处理一下噪点或者断连的地方,再将闭合空腔填充完整。如下,opening用于去除噪点,fillHull用于填充闭合空腔。如果有断连的地方,可以使用closing函数处理一下。
细胞核分割
有两个函数可以用来做分割,一个是bwlabel,另一个是watershed。
如果细胞都是各自游离的位置,那么bwlabel函数就可以很好的分开。一旦有细胞核黏连,那么watershed进行切割会更加理想。
下面以watershed为例进行:
# 阈值化
nuc_th <- nuc %>% thresh(w=10, h=10, offset=0.05) %>% opening(kern = makeBrush(5, shape = "diamond")) %>% fillHull()
# 区分
nmask <- watershed( distmap(nuc_th), 2 )
有两种展示方式:
display(colorLabels(nmask))
display(paintObjects(nmask, cells, col="red"))
泰森多边形分割单细胞
可以使用propagate来对胞质进行单细胞分割,从细胞核位置到细胞质位置进行分析,屏蔽掉非细胞的空腔位置。
ctmask <- opening(cel>0.1, makeBrush(5, shape='disc')) #阈值化
cmask <- propagate(cel, seeds=nmask, mask=ctmask) #单细胞分割
查看一下分割的效果:
# segment visualization
segmented <- paintObjects(cmask, cells, col='#ff00ff')
segmented <- paintObjects(nmask, segmented, col='#ffff00')
computeFeature函数用于计算单细胞的各项指标:
# 根据cmask计算表达量
computeFeatures.basic(
x = cmask %>% getFrame(1),
ref = cel %>% getFrame(1)
) %>% head
# b.mean b.sd b.mad b.q001 b.q005 b.q05 b.q095 b.q099
#1 0.15659996 0.03345964 0.04360588 0.10196078 0.10588235 0.15686275 0.20529412 0.22380392
#2 0.27868010 0.07614482 0.05232706 0.10588235 0.12784314 0.30980392 0.36470588 0.38541176
#3 0.19299299 0.07739029 0.07558353 0.10196078 0.10588235 0.16470588 0.33098039 0.34901961
#4 0.12109329 0.01324227 0.01162824 0.09137255 0.10196078 0.12156863 0.14117647 0.15294118
#5 0.05376344 0.01155338 0.01744235 0.03529412 0.03529412 0.05490196 0.07058824 0.07058824
#6 0.17997949 0.04603025 0.05814118 0.10196078 0.11372549 0.17647059 0.25490196 0.27450980
#...
# 计算细胞形态
computeFeatures.shape(cmask %>% getFrame(1)) %>% head
# s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
#1 194 67 9.030723 4.0627810 2.269099 16.275598
#2 473 88 12.583776 3.0285855 8.505285 18.975435
#3 333 91 12.143212 4.9556619 2.991797 21.763092
#4 231 61 8.389792 1.9535046 4.813758 12.861499
#5 31 16 2.681693 0.4283942 2.000000 3.162278
#6 1587 158 22.357467 3.2960345 15.564528 28.457854
#...
# 计算细胞位置信息
computeFeatures.moment(cmask %>% getFrame(1)) %>% head
# m.cx m.cy m.majoraxis m.eccentricity m.theta
#1 124.139175 4.288660 31.563604 0.9523654 0.15538879
#2 202.209302 10.723044 31.644121 0.6983248 0.88456709
#3 224.177177 7.300300 38.638433 0.9238340 -0.41367500
#4 503.536797 8.627706 21.587927 0.7170447 0.94331665
#5 490.000000 8.000000 7.326488 0.6933752 1.57079633
#6 21.199118 23.676749 52.513408 0.6364721 1.37670424
#...
做教学我们是认真的,如果你对我们的马拉松授课(直播一个月互动教学)有疑问,可以看完我们从2000多个提问互动交流里面精选的200个问答! 2021第二期_生信入门班_微信群答疑整理,以及 2021第二期_数据挖掘班_微信群答疑笔记