0

我有从 .ab1 文件中提取的测序运行信息（4 个核苷酸）。我希望能够将四个多峰高斯分布拟合到数据中（对应于 4 个不同的核苷酸） 数据是一个 csv 文件，有五列 - 索引列和其他四列对应于来自四个核苷酸的读取 -A，T， G和C。

 x=data.frame(read.csv(file.choose()))
 smooth1=ksmooth(x$index,x$A,kernel="normal",bandwidth=2)
 smooth2=ksmooth(x$index,x$C,kernel="normal",bandwidth=2)
 smooth3=ksmooth(x$index,x$G,kernel="normal",bandwidth=2)
 smooth4=ksmooth(x$index,x$T,kernel="normal",bandwidth=2)
 dsmooth1=diff(smooth1$y)
 dsmooth2=diff(smooth2$y)
 dsmooth3=diff(smooth3$y)
 dsmooth4=diff(smooth4$y)
 locmax1<-sign(c(0,dsmooth1))>0 & sign(c(dsmooth1,0))<0
 locmax2<-sign(c(0,dsmooth2))>0 & sign(c(dsmooth2,0))<0
 locmax3<-sign(c(0,dsmooth3))>0 & sign(c(dsmooth3,0))<0
 locmax4<-sign(c(0,dsmooth4))>0 & sign(c(dsmooth4,0))<0
 plot(x$index,x$A,xlim=c(900,950))
 lines(smooth1)
 lines(smooth2,col="green")
 lines(smooth3,col="blue")
 lines(smooth4,col="red")
 points(smooth1$x[locmax1],smooth1$y[locmax1],cex=3,c=2)
 points(smooth2$x[locmax2],smooth2$y[locmax2],cex=3,c=2)
 points(smooth3$x[locmax3],smooth3$y[locmax3],cex=3,c=2)
 points(smooth4$x[locmax4],smooth4$y[locmax4],cex=3,c=2)

为了进一步定位峰，我使用了以下

peaks=function(x) {
   modes=NULL
   for ( i in 2:(length(x)-1) ){
   if ( (x[i] > x[i-1]) & (x[i] > x[i+1]) ) {
     modes=c(modes,i)
   }
   }
   if ( length(modes) == 0 ) {
   modes = 'This is a monotonic distribution'
   }
   return(modes)
   }

  x$A[peaks(x$A)] #similarly, for T,G and C

某些点有多个峰值，我需要编写一个代码来找到这样的位置，这些位置具有多个高斯分布的峰值（对应于来自多个核苷酸的信号）。有没有办法在 R 中做到这一点？

4

1 回答 1

1

你基本上很合适混合模型拟合到您的数据：四个高斯的混合。我建议你阅读这些。有比平滑和检测峰值更复杂的方法来处理这些问题（这可能在很大程度上取决于您的内核宽度 - 所以如果您进行平滑，您应该进行一些敏感性分析并检查您的结果如何随不同的内核和内核宽度变化）。

R的mixtools包应该很有用。

于 2014-12-17T08:24:06.650 回答