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的最大距离。任一点(x0lovebet体育官网,,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, ...)

  平板电脑,拟合表面是从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, ...)

  三星GALAXY Tab,拟合表面是从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的长度
lovebet体育官网 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")

lovebet体育官网 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")

lovebet体育官网 3

本来,也想看一下amplicon 长度len 对RC的影响,但是影响不大
lovebet体育官网 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的长度
lovebet体育官网 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")

lovebet体育官网 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")

lovebet体育官网 7

本来,也想看一下amplicon 长度len 对RC的震慑,不过影响不大
lovebet体育官网 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")