IT狗

R言语经过loess去除某个变量对数据的干扰

  当咱们想研讨不一样sample的某个变量A之间的差别时,常常会由于别的一些变量B对该变量的固有干扰,而干扰不一样sample变量A的对比,这个时间须要对sample变量A进行尺度化以后能力进行对比。尺度化的要领是对sample 的 A变量和B变量进行loess归来,拟合变量A关于变量B的函数 f(b),f(b)则暗示在B的干扰下A的实践取值,A-f(B)(A对f(b)残差)就可以去掉B变量对A变量的干扰,这时候残差值就可以作为尺度化的A值在不一样sample之间进行对比。

Loess个别地方加权多项式归来

  LOWESS最后由Cleveland 提出,后又被Cleveland&Devlin及其他很多人生长。在R中loess 函数是以lowess函数为底子的更庞大功用更壮大的函数。重要头脑为:在数据纠合的每一点用低维多项式拟合数据点的一个子集,并预计该点四周自变量数据点所对应的因变量值,该多项式是用加权最小二乘法来拟合;离该点越远,权重越小,该点的归来函数值就是这个个别地方多项式来获得,而用于加权最小二乘归来的数据子集是由最近邻要领肯定。
  最大长处:不须要事前设定一个函数来对一切数据拟合一个模子。而且可以对统一数据进行屡次不一样的拟合,先对某个变量进行拟合,再对另一变量进行拟合,以索求数据中大概存在的某种干系,这是平凡的归来拟合没法做到的。

LOESS光滑要领

  1. 以x0为中央肯定一个区间,区间的宽度可以机动掌握。具体来说,区间的宽度取决于q=fn。当中q是介入个别地方归来考察值的个数,f是加入个别地方归来考察值的个数占考察值个数的比重,n是考察值的个数。在实践运用中,常常先选定f值,再按照f和n肯定q的取值,一样平常状况下f的取值在1/3到2/3之间。q与f的取值一样平常没有肯定的原则。增大q值或f值,会引发光滑值光滑水平添加,关于数据中前在的纤细变更形式则分辨率低,但噪声小,而对数据中大的变更形式的施展阐发则对比好;小的q值或f值,曲线暗沉,分辨率高,但噪声大。没有一个尺度的f值,对比明智的方法是不断的调试对比。
  2. 界说区间内一切点的权数,权数由权数函数来肯定,比方立方加权函数weight = (1 - (dist/maxdist)^3)^3),dist为间隔x的间隔,maxdist为区间内间隔x的最大间隔。任一点(x0,y0)的权数是权数函数曲线的高度。权数函数应包罗以下三个方面特征:(1)加权函数上的点(x0,y0)具有最大权数。(2)当x离开x0(时,权数渐渐淘汰。(3)加权函数以x0为中央对称。
  3. 对区间内的散点拟合一条曲线y=f(x)。拟合的直线反应直线干系,靠近x0的点在直线的拟合中起到重要的感化,区间外的点它们的权数为零。
  4. x0的光滑点就是x0在拟合出来的直线上的拟合点(y0,f( x0))。
  5. 对一切的点求出光滑点,将光滑点毗连就获得Loess归来曲线。

R言语代码

 loess(formula, data, weights, subset, na.action, model = FALSE,       span = 0.75, enp.target, degree = 2,       parametric = FALSE, drop.square = FALSE, normalize = TRUE,       family = c("gaussian", "symmetric"),       method = c("loess", "model.frame"),       control = loess.control(...), ...)

  formula是公式,比方y~x,可以输出1到4个变量;
  data是放着变量的数据框,假如data为空,则在环境中寻觅;
  na.action指定对NA数据的处置,默许是getOption("na.action");
  model是不是前往模子框;
  span是alpha参数,可以掌握光滑度,相当于下面所述的f,关于alpha小于1的时间,区间包罗alpha的点,加权函数为立方加权,大于1时,运用一切的点,最大间隔为alpha^(1/p),p 为说明变量;
  anp.target,界说span的备选要领;
  normalize,对多变量normalize到统一scale;
  family,假如是gaussian则运用最小二乘法,假如是symmetric则运用双权函数进行再降低的M预计;
  method,是顺应模子大概仅仅采集模子框架;
  control进一步更高等的掌握,运用loess.control的参数;
  别的参数请本身拜见manual而且查找材料

loess.control(surface = c("interpolate", "direct"),          statistics = c("approximate", "exact"),          trace.hat = c("exact", "approximate"),          cell = 0.2, iterations = 4, ...)

  surface,拟合轮廓是从kd数进行插值仍是进行正确计较;
  statistics,统计数据是正确计较仍是近似,正确计较很慢
  trace.hat,要跟踪的光滑的矩阵正确计较或近似?倡议运用跨越1000个数据点迫近,
  cell,假如经过kd树最大的点进行插值的近似。大于cell floor(nspancell)的点被细分。
  robust fitting运用的迭代次数。

predict(object, newdata = NULL, se = FALSE,    na.action = na.pass, ...)

  object,运用loess拟合出来的对象;
  newdata,可选数据框,在内里寻觅变量并进行展望;
  se,是不是计较尺度偏差;
  对NA值的处置

实例

  生物数据分析中,咱们想检察PCR扩增出来的扩增子的测序深度之间的差别,但不一样的扩增子的扩增效力遭到GC含量的干扰,故此咱们起首应当清除掉GC含量对扩增子深度的干扰。

数据

amplicon 测序数据,处置后获得的每一个amplicon的深度,每一个amplicon的GC含量,每一个amplicon的长度

先用loess进行曲线的拟合

gcCount.loess <- loess(log(RC+0.01)~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)

画出拟合出来的曲线

predictions1<- predict (gcCount.loess,RC_DT$GC)#plot scatter and line plot(RC_DT$GC,log(RC_DT$RC+0.01),cex=0.1,xlab="GC Content",ylab=expression(paste("log(NRC"["lib"],"+0.01)",sep="")))lines(RC_DT$GC,predictions1,col = "red")

取残差,去除GC含量对深度的干扰

#sustract the influence of GCresi <- log(RC_DT$RC+0.01)-predictions1RC_DT$RC <- resisetkey(RC_DT,GC)

这时候RC_DT$RC就是normalize以后的RC
绘图显现nomalize以后的RC,并将拟合的loess曲线和normalize以后的数据保管

#plot scatter and line using Norm GC dataplot(RC_DT$GC,RC_DT$RC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))gcCount.loess <- loess(RC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)save(gcCount.loess,file="/home/ywliao/project/Gengyan/gcCount.loess.Robject")predictions2 <- predict(gcCount.loess,RC_DT$GC)lines(RC_DT$GC,predictions2,col="red")save(RC_DT,file="/home/ywliao/project/Gengyan/RC_DT.Rdata")

固然,也想看一下amplicon 长度len 对RC的干扰,不外干扰不大

悉数代码以下:

library(data.table)load("/home/ywliao/project/Gengyan/RC_DT.Rdata")####loess GC vs RC####gcCount.loess <- loess(log(RC+0.01)~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)predictions1<- predict (gcCount.loess,RC_DT$GC)#plot scatter and line plot(RC_DT$GC,log(RC_DT$RC+0.01),cex=0.1,xlab="GC Content",ylab=expression(paste("log(NRC"["lib"],"+0.01)",sep="")))lines(RC_DT$GC,predictions1,col = "red")#sustract the influence of GCresi <- log(RC_DT$RC+0.01)-predictions1RC_DT$RC <- resisetkey(RC_DT,GC)#plot scatter and line using Norm GC dataplot(RC_DT$GC,RC_DT$RC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))gcCount.loess <- loess(RC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)save(gcCount.loess,file="/home/ywliao/project/Gengyan/gcCount.loess.Robject")predictions2 <- predict(gcCount.loess,RC_DT$GC)lines(RC_DT$GC,predictions2,col="red")save(RC_DT,file="/home/ywliao/project/Gengyan/RC_DT.Rdata")####loess len vs RC###setkey(RC_DT,Len)len.loess <- loess(RC_DT$RC~RC_DT$Len, control = loess.control(surface = "direct"),degree=2)predictions2<- predict (len.loess,RC_DT$Len)#plot scatter and line plot(RC_DT$Len,RC_DT$RC,cex=0.1,xlab="Length",ylab=expression(paste("log(RC"["GC"],"+0.01)",sep="")))lines(RC_DT$Len,predictions2,col = "red")

此文由 IT狗 编辑,本网站所发布展示的作品/文章版权归原作者所有,任何商业用途均须联系作者!

相关推荐

评论 暂无评论