你好,游客 登录 注册 搜索
背景:
阅读新闻

用R语言搞定了Dynamica的第一题!

[日期:2015-03-13] 来源:科学网  作者:罗志达 [字体: ]

Coursera今年3月初开设了一门DynamicalModeling Methods for Systems Biology,怀着对建模的极大热忱二话不说就报了这个课程。按理说第一周的课无非就是水水课程简介什么的,可没想第一周的第一道练习题就立马给我个下马威!

 

计算机生成了可选文字:Question 1 
These exercises are meant to give you practice with the MATLAB programming environment and are designed tor you to demonstrate your 
competence perTorming oaslc computations using this sortware. 
Part 1:Array computations tor simple data analysis 
For this question you will read in data that were obtained experimentally using a laser scanning contocal microscope operating in line-scan 
mode. Tne manipulations required in tne assignment include:1 accessing onry a particular portion or the data (a denned spatial position). 2) 
spatial averaging to reduce noise. 3) normalization. and 4) converting the data into different units First, tor those of you who are interested, a 
few words about what me data represent In tms experiment a rat ventricular myocyte was loaded witn the sensitive indicator nuo-3. and 
a confocal microscope was used to record changes in intracellular (Caz• l. The microscope was operated in line scan

 

计算机生成了可选文字:The data appear as TollOWSj 
To perform tnis analysis. you have to take the following steps:a Oowtload flash4jpg trom here 
Read in tne Tile:data imread('flash4.jpg', Ipg') 
Take a 100K at it using tne imagesc function. Onent tne data so that time runs trom lett to right Transpose tne matrix it necessary. 
c, Average over the regvon ot the IJVffash Store this in the variable nasty This should have dimensions 1 x 634, 
Average over a control region that does not contain the rash. Store this in the variable nonash 
e. The fluorescence units are arbitrary, since the number depends on laser intensity. dye concentration. microscope detector gain, etc. Thus. 
we are interested in relative changes in fluorescence (F ) over the baseline value (Fo ) in a resting cell_ Convert trom raw fluorescence to units 
or F,F, by normalizing nasn and nonasn to tne average nuorescence in a region witn no activity (Hint between lines 70 and 100 is a good 
region). 
t TO a first approximation. one can assunne that the dye used in this experiment. fluo-3. only emits fluorescence when is bound to Ca2' 
Thus. where KO is the dissociation constant ot the dye. and (Flu031TOTAL is the dye fluorescence is proportional to; 
tCa2•

 

计算机生成了可选文字:If one makes a reasonable assumption for baseline I (eg_ 100 nm). the following equation can be used to convert from a ratio R (units Of 
) to I in units of concentration 
this equation to convert lash and notlash trom units otF,F, to units ot l, in nM_ 
You can assume - 100 nM. and 700 nM. Keep in mind that your variables nasn aru3 nonasn are not just scalars (numbers) 
but are defined at many time points 
Plot both versus time, on the same scale, in different colors 
wnicn of tne following plots can represent tne output plot if you assume - 150 nM and K' 
- 1000

 

 

什么乱七八糟啊()┻━┻

 

()b小白无人权,我忍!!

 

Google和百度,我又来啦|O|~~

 

来来来,先让我们先来了解一下计算机是如何储存和表示图像的:

 

==================================================================================

图像的存储原理:

 

从结构上,图形文件分为位图和矢量图。

  • 在 位图中,图像由许多的屏幕小点(像素)组成,这些小点对应显存中的“位”,而就是这些“位”决定了像素的图形属性,如像素的颜色、灰度、明暗对比度等。当 一个像素所占的位数(计算机物理存储空间)越多时,它所能表现的颜色就更丰富。从整体看,图像的色彩就更艳丽,分辨率更高。位图中所分的二位图、八位图等 正是指像素所占的位数。

  • 当位图被放大或缩小时,由于像素的数量没有改变,图像的分辨率就会降低,图像的外观自然大打折扣。从这点来说,位图的缺点显而易见——分辨率的固定导致大分辨率的清晰图像占用大量空间;像素的分散性使动态图像的表达显得困难(例如看VCD时出现马赛克,就是像素丢失造成的)。

  • 因 此一种新的图形格式应运而生——矢量图。顾名思义,矢量图就是用矢量代替位图中的“位”。简单来说,就是用矢量给图的几何部分做标记。它的表达方式是先用 语句调用调色板描述背景,再用带矢量的数学公式来描述圆圈的大小、形状等,这就使得图形的放大、缩小和移动变得十分简单,仅仅把公式中的矢量变量改一改就 可以。因此,矢量图有无限放缩而不失真,占用空间少等优点。

 

计算机生成了可选文字:点 阵 图 
斂 人 后 的 点 阵 图 
矢 量 图 
放 大 后 的 矢 量 图

 

 

这里我们回过头来说说位图是如何存储颜色的。

我们知道,自然界中所有颜色都可以由红、绿、蓝(RGB)组合而成。所谓的816位色中的位数,是指每个颜色(红、绿、蓝)中每个像素点可以存储的灰阶值。当一副图中每个像素被赋予不同的RGB值时,图像就能呈现出五彩缤纷的色彩。

 

当然,关于数学图像处理的学问同样博大精深,有兴趣的读者可以参考相关书籍。

 

=================================================================================

 

事实上,之所以自己还有点底气上这门课,是因为之前自己学了点R语言。所以尽管课程要求我们用matlab,不过我还是想用R去实现。但是这次的问题却是自己之前完全没遇过的!不过,用R进行图像处理真是一件很geek的事,这里特别推荐一篇很好玩的博文http://chenangliu.info/cn/use-r-for-fun-image/

 

进行图片处理的第一步自然是将图片读进来。CRAN上可供选择的有

CRAN - Package jpeg, CRAN - Package png, CRAN - Package bmp, CRAN - Package readbitmap。不过这四个包仅仅是读取图像文件,基本没有其他的图片处理扩展能力。其他的还有biOpsbioconductor上的EBImage等等。

 

那么下载安装并加载EBImage后,我们就可以开始读取图像了。

 

回到CEric的课程练习上。题目很有挑战性(_;)? [不懂] ,大意是在rat ventricular myocyte中加入了Ca2+ sensitive indicator, fluo-3,然后放在共聚焦显微镜下拍照,照片如下:

 

计算机生成了可选文字:

 

这个细胞同时加入了Ca2+的缓冲剂NP-EGTA,当它暴露在高强度的紫外线下时,它就不能结合Ca2+。因此紫外光可以作为开关:当照射紫外光后,NP-EGTA可以释放出大量游离的Ca2+进入细胞,触发荧光。紧接着,细胞自身的电传递会引起第二次的Ca2+浓度上升。我们的目的就是确定第一次Ca2+上升对第二次Ca2+的影响。

 

题目中特别强调了共聚焦显微镜是以“line scan”的模式成像,所以图片的一个维度是空间,另一个是时间。这句极具干扰性的话让从没接触过此技术的本人挠破脑袋不明所指。另外题目中对这张照片也缺乏描述,很难让人读懂这张图,因此只能求助于课程论坛。幸好还是有大神给出了提示:这是一张512x634的图片,其中634这个方向是时间轴(我想应该可以理解为在成像的时候显微镜打出一道线性的激光,沿着634这个方向边扫描边成像吧)。

 

这里当时我还遇到一个问题。前面说过这个课程实际上要求用matlab进行编程。不过因为在入门的时候学的只有R语言,所以碍于惰性也一心只想着看看能否在R上实现。但是接下来这道题目提示我们用matlabimagesc函数,这可愁坏了我——这个函数是干嘛的?不用这个函数这道题还做不做了?!

 

天无绝人之路!就在山穷水尽之时,我突然惊喜地发现R里面居然有一个叫matlab的包,这个包里同样也有imagesc函数!!

 

接下来就容易多了:

先用jpeg package中的readJPEG()函数读取保存的图片(图片保存在Rworkspace下的话就不用麻烦输入路径了);转置一下矩阵,让列表示时间轴。

> img <- readJPEG("flash4.jpg")

> T<-t(img)

 

EBImage package中的display(XXXmethod="raster")可以在R的作图窗口中显示。(不加“raster”的话,默认的显示方式是浏览器显示。)

调用matlabimagesc函数,原来不过如此!

> imagesc(T)

 

计算机生成了可选文字:

 

 

从这张图我们就可以很清楚的看到发荧光的部分。论坛的大神还特别说题目中其实已经暗示了那个独立在远处的一小块光斑就是紫外线照射的部分((_)好吧,碍于我拙劣的英语理解能力),然后线性的那条自然就是第二次发光。

 

那这个时候我们先把第一次发光的那部分截出来。因为这时候我们已经把图片读入并储存为矩阵形式,所以我自己粗略估计在第210-280行:

> display(T[210:280,],method="raster")

 

计算机生成了可选文字:

 

 

对这几行求平均,存入到一个1x634的变量flash_mean

flash_mean <-apply(T[210:280,],2,mean)

对除210-280的其他行求平均,存入到另一个1x634的变量noflash_mean中;

> noflash_mean <-apply(T[-(210:280),],2,mean)

 

那么这样处理后,图片就会变成这个样子((")丑是丑了点)!

> merge <-c(rep(noflash_mean,209),rep(flash_mean,71),rep(noflash_mean,232))

> new_T <-matrix(merge,ncol=634,byrow = TRUE)

>display(new_T,method="raster")

 

 

计算机生成了可选文字:

 

 

 

题目还特别友好地提示在第70-100行可以作为背景值,看原图的话其实就是穿过那条亮线下方的小缺口那一区域!所以我们同样求出这几行的均值,并存入到1x634的变量baseline中;

> baseline <-apply(T[70:100,],2,mean)

 

最后我们只要将flashnoflash分别除以baseline,就完成了normalization的过程。

> R_flash <-flash_mean/baseline
> R_noflash <- noflash_mean/baseline

 

接下来按照题目的提示,我们代入到公式

计算机生成了可选文字:[Сап 1-

其中KD=1000Ca2+_baseline=150

 

那么对两条曲线分别作图的话,就可以得到正确答案啦~\()/~

 

 

计算机生成了可选文字:100 
200 
300 
400 
500 
600 
700 
BOO

 

 

 





收藏 推荐 打印 | 录入:Cstor | 阅读:
本文评论   查看全部评论 (4)
表情: 表情 姓名: 字数
点评:
       
评论声明
  • 尊重网上道德,遵守中华人民共和国的各项有关法律法规
  • 承担一切因您的行为而直接或间接导致的民事或刑事法律责任
  • 本站管理人员有权保留或删除其管辖留言中的任意内容
  • 本站有权在网站内转载或引用您的评论
  • 参与本评论即表明您已经阅读并接受上述条款