R语言通过loess去除某个变量对数据的影响。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之间进行比。

  当我们怀念研究不同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局部加权多项式回归

  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回归曲线。

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值的拍卖

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含量对扩增子深度的熏陶。

实例

  生物数据解析着,我们怀念查PCR扩增出来的扩增子的测序深度以内的差异,但不同的扩增子的扩增效率受到GC含量之熏陶,因此我们首先应排除掉GC含量对扩增子深度的影响。

数据

amplicon
测序数据,处理后拿走的每个amplicon的深度,每个amplicon的GC含量,每个amplicon的长
大红鹰葡京会 1
先用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")

大红鹰葡京会 2

取残差,去除GC含量对纵深的震慑

#sustract the influence of GC
resi <- log(RC_DT$RC+0.01)-predictions1
RC_DT$RC <- resi
setkey(RC_DT,GC)

此时RC_DT$RC就是normalize之后的RC
绘画图显示nomalize之后的RC,并拿拟合的loess曲线和normalize之后的数量保存

#plot scatter and line using Norm GC data
plot(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")

大红鹰葡京会 3

本来,也想看一下amplicon 长度len 对RC的熏陶,不过影响不要命
大红鹰葡京会 4

不折不扣代码如下(经过修改,可能同方了配合):

library(data.table)

load("/home/ywliao/project/Gengyan/RC_DT.Rdata")
RRC_DT <- RC_DT[Type=="WBC" & !is.na(RC),]

lst <- list()
for (Samp in unique(RC_DT$Sample)){
RC_DT <- RRC_DT[Sample==Samp]
####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 GC
resi <- log(RC_DT$RC+0.01)-predictions1
RC_DT$NRC <- resi
setkey(RC_DT,GC)
#plot scatter and line using Norm GC data
#plot(RC_DT$GC,RC_DT$NRC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))
gcCount.loess <- loess(NRC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
predictions2 <- predict(gcCount.loess,RC_DT$GC)
#lines(RC_DT$GC,predictions2,col="red")
lst[[Samp]] <- RC_DT
}
NRC_DT <- rbindlist(lst)
save(RC_DT,file="/home/ywliao/project/Gengyan/NRC_DT.Rdata")

####loess len vs RC###
setkey(RC_DT,Len)
len.loess <- loess(RC_DT$NRC~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$NRC,cex=0.1,xlab="Length",ylab=expression(paste("log(RC"["GC"],"+0.01)",sep="")))
lines(RC_DT$Len,predictions2,col = "red")

数据

amplicon
测序数据,处理后拿走的每个amplicon的深,每个amplicon的GC含量,每个amplicon的长
大红鹰葡京会 5
优先用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")

大红鹰葡京会 6

取残差,去除GC含量对纵深的影响

#sustract the influence of GC
resi <- log(RC_DT$RC+0.01)-predictions1
RC_DT$RC <- resi
setkey(RC_DT,GC)

此时RC_DT$RC就是normalize之后的RC
打图显示nomalize之后的RC,并将拟合的loess曲线和normalize之后的数量保存

#plot scatter and line using Norm GC data
plot(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")

大红鹰葡京会 7

理所当然,也想看一下amplicon 长度len 对RC的震慑,不过影响不充分
大红鹰葡京会 8

全副代码如下(经过改,可能同方了配合):

library(data.table)

load("/home/ywliao/project/Gengyan/RC_DT.Rdata")
RRC_DT <- RC_DT[Type=="WBC" & !is.na(RC),]

lst <- list()
for (Samp in unique(RC_DT$Sample)){
RC_DT <- RRC_DT[Sample==Samp]
####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 GC
resi <- log(RC_DT$RC+0.01)-predictions1
RC_DT$NRC <- resi
setkey(RC_DT,GC)
#plot scatter and line using Norm GC data
#plot(RC_DT$GC,RC_DT$NRC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))
gcCount.loess <- loess(NRC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)
predictions2 <- predict(gcCount.loess,RC_DT$GC)
#lines(RC_DT$GC,predictions2,col="red")
lst[[Samp]] <- RC_DT
}
NRC_DT <- rbindlist(lst)
save(RC_DT,file="/home/ywliao/project/Gengyan/NRC_DT.Rdata")

####loess len vs RC###
setkey(RC_DT,Len)
len.loess <- loess(RC_DT$NRC~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$NRC,cex=0.1,xlab="Length",ylab=expression(paste("log(RC"["GC"],"+0.01)",sep="")))
lines(RC_DT$Len,predictions2,col = "red")

相关文章

admin

网站地图xml地图