R语言如何查询正弦波的波峰值

Python010

R语言如何查询正弦波的波峰值,第1张

用内置函数 optim()

optim(par,fun,lower,upper,method) 大致用到这5个参数

par是初始值,你选离你峰值差不远的x

fun是生成你正弦波的函数

lower和upper定义域

method用 "BFGS"牛顿迭代法,或者"L-BFGS-B"升级版牛顿迭代法。

以下是得到的结果,我用f(x)=x^2-2*x+1试了以下

>optim(3,fun,lower=-5,upper=5,method="BFGS")

$par

[1] 1 # x值

$value

[1] 0 # y值

$counts

function gradient # 不知道是什么

44 44

$convergence

[1] 52 # 算了多少次收敛

目录

在之前的文章里介绍了如何通过直方图来可视化等位杂合碱基的比例来判断物种的染色体倍数性。在本文里会继续向下挖掘,介绍如何可视化染色体上的拷贝数变化(CNVs)。

前文一样的操作,使用包自带的数据。

我们需要去除过高和过低深度的数据。和前文的操作一样,提取vcf文件里的深度数据"AD"。

然后过滤出10%~90%的数据,当然此处可以根据实际情况进行微调。然后对第一种出现频率最高的碱基进行可视化。(一般情况下一个位点上会有两种碱基,具体参考前文。)

同样也可以对出现频率第二高的碱基进行同样的操作,这里节约篇幅就省略了。

为了避免复杂的基于AD比例的模型假设,程序里设计了非参数估计法来计算峰值。计算完了以后可以直接对染色体进行拆分以后可视化进行校验。

根据尺寸把染色体分割成合适的大小

然后用 freq_peak 函数计算峰值。并对数据进行处理,去掉负数和Na值。

计算到此为止,可以可视化实际数据来验证计算的正确性。

仔细想一下,峰值计算的结果其实就是CNV的结果。这里根据窗口大小把染色体分成了若干段。(那么是不是可以给每一段 CDS进行细分然后计算出每一个CDS的具体数字呢????)

当然也可以把所有样本组合到一起。