R语言:DNA序列比对后计算遗传距离(P-distance)

Python015

R语言:DNA序列比对后计算遗传距离(P-distance),第1张

在R语言中找到了计算遗传距离的函数dist.dna()但是不知道在R里面如何利用循环批量处理文件计算遗传距离。想到了利用python来调用R函数的方法,查找相关教程发现需要用到rpy2模块。

easy_install rpy2 报错(看不懂报错内容);

pip install rpy2 报错(提示需要更新pip到pip19.0.3);

利用 python -m pip install --upgrade pip 更新pip报错(看不懂报错内容);

利用 https://pip.pypa.io/en/stable/installing/ 教程安装pip成功更新。

使用 pip install rpy2 安装依旧报错(看不懂报错内容);

尝试教程 https://blog.csdn.net/suzyu12345/article/details/51476321 安装rpy2,提示 rpy2-2.9.5-cp37-cp37m-win_amd64.whl is not a supported wheel on this platform

在rpy2主页 https://rpy2.bitbucket.io/ 发现一句话 Releasend source packages are available on PyPi. Installing should be as easy * as

(*:except on Windows)

这意思是在windows系统使用pip安装不太容易吗?

找到了教程 rpy2:在python中调用R函数的一个实例 ;发现其中rpy2的安装使用的是conda,自己也尝试在windows的DOS窗口下使用 conda install rpy2 成功。但是结尾处提示了一句 此时不应有do 不明白是什么意思

在python中加载R包查到可以使用

加载R语言自带的包时没有遇到问题;但是加载需要额外安装的包时遇到了报错

按照教程 https://blog.csdn.net/arcers/article/details/79109535 使用 conda install -c r r-ggplot2 安装需要用到的包解决问题

一个简便描述序列分歧大小的测度是两条核苷酸序列中不同核苷酸位点的比例 P = nd/n;

nd为检测两条序列间不同核苷酸数;n为配对总数;P成为核苷酸间的p距离

关于谱能量,有这样一种解释,你可以试着去算一算信号可以分成能量信号与功率信号,非周期能量信号具有能量谱密度,是傅立叶变换的平方,功率信号具有功率谱密度,其与自相关函数是一对傅立叶变换对,等于傅立叶变换的平方/区间长度。不能混淆。能量信号是没有功率谱的。胡广书老师的书上找到这么一段话,“随机信号在时间上是无限的,在样本上也是无穷多,因此随机信号的能量是无限的,它应是功率信号。功率信号不满足付里叶变换的绝对可积的条件,因此其付里叶变换是不存在的。如确定性的正弦函数的付里叶变换是不存在,只有引入了冲激函数才求得其付里叶变换。因此,对随机信号的频谱分析,不再简单的是频谱,而是功率谱。”对于确定性信号而言,里面存在能量信号,是没有功率谱密度的,也存在功率信号,是有功率谱密度的。所以信号的频谱与是否是确定性信号没有必然联系。以下论点来源于研学论坛:频谱是信号的傅立叶变换。它描述了信号在各个频率上的分布大小。频谱的平方(当能量有限,平均功率为0时称为能量谱)描述了信号能量在各个频率上的分布大小。计算过程中,都是通过样本数据的快速傅立叶变换来计算。但不同的是,信号的频谱是复数,包含幅频响应和相频响应,重复计算时的结果基本相同。而随机信号的功率谱也可以对数据进行FFT,但必须计算模值的平方,因为功率谱是实数。而且换一组样本后,计算的结果略有不同,因为随机信号的样本取值不同。要得到真实的功率谱必须进行多次平均,次数越多越好。根据parseval定理,信号傅氏变换模平方被定义为能量谱,即单位频率范围内包含的信号能量。自然,能量跟功率有一个时间平均的关系,所以,能量谱密度在时间上平均就得到了功率谱。matlab实现经典功率谱估计fft做出来是频谱,psd做出来是功率谱;功率谱丢失了频谱的相位信息;频谱不同的信号其功率谱是可能相同的;功率谱是幅度取模后平方,结果是个实数matlab中自功率谱密度直接用psd函数就可以求,按照matlab的说法,psd能实现Welch法估计,即相当于用改进的平均周期图法来求取随机信号的功率谱密度估计。psd求出的结果应该更光滑吧。1、直接法:直接法又称周期图法,它是把随机序列x(n)的N个观测数据视为一能量有限的序列,直接计算x(n)的离散傅立叶变换,得X(k),然后再取其幅值的平方,并除以N,作为序列x(n)真实功率谱的估计。Matlab代码示例:clearFs=1000%采样频率n=0:1/Fs:1%产生含有噪声的序列xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))window=boxcar(length(xn))%矩形窗nfft=1024[Pxx,f]=periodogram(xn,window,nfft,Fs)%直接法plot(f,10*log10(Pxx))2、间接法:间接法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。Matlab代码示例:clearFs=1000%采样频率n=0:1/Fs:1%产生含有噪声的序列xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))nfft=1024cxn=xcorr(xn,'unbiased')%计算序列的自相关函数CXk=fft(cxn,nfft)Pxx=abs(CXk)index=0:round(nfft/2-1)k=index*Fs/nfftplot_Pxx=10*log10(Pxx(index+1))plot(k,plot_Pxx)3、改进的直接法:对于直接法的功率谱估计,当数据长度N太大时,谱曲线起伏加剧,若N太小,谱的分辨率又不好,因此需要改进。3.1、Bartlett法Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。Matlab代码示例:clear;Fs=1000n=0:1/Fs:1xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))nfft=1024window=boxcar(length(n))%矩形窗noverlap=0%数据无重叠p=0.9%置信概率[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p)index=0:round(nfft/2-1)k=index*Fs/nfftplot_Pxx=10*log10(Pxx(index+1))plot_Pxxc=10*log10(Pxxc(index+1))figure(1)plot(k,plot_Pxx)pausefigure(2)plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc])3.2、Welch法Welch法对Bartlett法进行了两方面的修正,一是选择适当的窗函数w(n),并再周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时,可使各段之间有重叠,这样会使方差减小。Matlab代码示例:clearFs=1000n=0:1/Fs:1xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n))nfft=1024window=boxcar(100)%矩形窗window1=hamming(100)%海明窗window2=blackman(100)%blackman窗noverlap=20%数据无重叠range='half'%频率间隔为[0 Fs/2],只计算一半的频率[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range)[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range)[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range)plot_Pxx=10*log10(Pxx)plot_Pxx1=10*log10(Pxx1)plot_Pxx2=10*log10(Pxx2)figure(1)plot(f,plot_Pxx)pausefigure(2)plot(f,plot_Pxx1)pausefigure(3)plot(f,plot_Pxx2)

这里给出以前写的一个示例。

下面以使用刀切法估计正态分布N(0,25)的方差为例来说明刀切法的实现和优势。这里,我们将样本的修正方差视为关于总体方差的一个估计量。

##1)产生随机变量x~N(0,25),产生100个样本

set.seed(206)

x<-rnorm(100,0,5)

##2)运用刀切法进行重抽样

jack<-function(x){

jackknife<-0

for(i in 1:length(x)){

jackknife[i]=length(x)*var(x)-(length(x)-1)/length(x)*sum(var(x[-i]))

}

jackknife

}

##3)计算刀切法得到的抽样方差

mean(jack(x))/length(x)

##4)全样本简单方差估计进行比较

var(x)

结果如下:

>mean(jack(x))/length(x)

[1] 23.26443

>var(x)

[1] 23.49705

由以上结果可知,刀切法得到的方差估计比全样本得到的方差估计更接近真实值。