用程序怎么计算下面公式

zhagnqiang829

怎么计算D2…
    ，（i,j =a,c,g,t）
其中 1/ln(2) = 1.442695040889
一  下面是所有出现的概率的情况( i,j 所有组合的16中情况概率，以及Pi,Pj的4中情况)
Pac        0.091116173
Pcc        0.054669704
Pag        0.102505695
Ptg        0.022779043
Pat        0.061503417
Ptc        0.06833713
Paa        0.141230068
Pga        0.104783599
Ptt        0.027334852
Pct        0.027334852
Pcg        0.045558087
Pgg        0.036446469
Pta        0.029612756
Pca        0.118451025
Pgc        0.034168565
Pgt        0.031890661

Pa        0.396355353
Pc        0.248291572
Pg        0.207289294
Pt        0.148063781

q1208c

这不就是把这个 ∑ 展开么?
如果我的数学老师教得没错的话, 不就是个 加法么? 两层的循环加一下好啦. 结果再除以 ln(2) 不就 OK了?

zhagnqiang829

你帮我改下代码吧，我运行总出错。
#!usr/bin/perl -w
use strict;
open IN,"file.txt";
while(<IN>){

if(/P(\w)(\w) (\d+\.\d+)){
%s=('a',0,'g',1,'c',2,'t',3);
$P[$s{$1}][$s{$2}]=$3;
}
}


for $i(0..$#P){
for $j(0..$#{$P[0]}){
$D2=1/1.44($P[$i][$j].$P[$i][$j]-2$P[$i][$j].$P[$i]$P[$j]+$P[$i]$P[$i]$P[$j]$P[$j])/$P[$i]$P[$j];
}
}


回复 2# q1208c


   

zhagnqiang829

唉，弄不出来。我还是那C++自己处理吧，由原始序列直接处理，分享下代码：


#include "stdio.h"
#include"stdlib.h"
#include"math.h"
#include"string.h"
#include"ctype.h"
#include"conio.h"

void main()
{
    char  a[300000],b[16][3]={"AA","AC","AG","AT","CA","CC","CG","CT","GA","GC","GG","GT","TA","TC","TG","TT"},c[3]={'\0','\0','\0'};
   FILE *fp1,*fp2;
   
   long int i,j, k=0,m=0;
  
   double   d[16]={0.0},e[4]={0.0},n=2.0,g=0.0,D2=0,f[16]={0.0};

   fp1=fopen("F:\\计算D2的序列.txt","r");
   fp2=fopen("F:\\D2值.txt","w");
   
  

       

          

   do
   {
           fgets(a,300000,fp1);

      
          if(a[0]=='@')
                  break;
       if((a[0]!='>')&&(a[0]!='>'))
          
           {  for(i=0;i<16;i++)
                     d[i]=0.0;
         for(i=0;i<4;i++)
               e[i]=0.0;
               g=0.0;
                D2=0;
                 m=0;
                   for(i=0;i<300000;i++)
                   {
               
               if(a[i]=='\n')
                           {
                                   //fprintf(fp2,"%s",a);
                                   m=i;
                               break;

                           }
                                   if(a[i]=='A')
                    e[0]++;
                   if(a[i]=='C')
                    e[1]++;
                   if(a[i]=='G')
                    e[2]++;

                   if(a[i]=='T')
                    e[3]++;

                                   c[0]=a[i];
                                   c[1]=a[i+1];
                   //fprintf(fp2,"%s\n",c);

                                   for(j=0;j<16;j++)
                                   {
                                           if(strcmp(b[j],c)==0)
                                           {
                                                   d[j]++;
                                                   break;
                                           }
                                   }

                  
                                   c[0]='\0';
                               
                                   c[1]='\0';

                   }
                   fprintf(fp2,"%d\n",m);
                   for(i=0;i<4;i++)
                   {
                           e[i]=e[i]/m;

                     fprintf(fp2,"%f\t",e[i]);
                   }
                  
                 
           
           fprintf(fp2,"\n");
                         c[0]=a[m-1];
                                          c[1]=a[0];
                  // fprintf(fp2,"%s\n",c);

                                   for(j=0;j<16;j++)
                                   {
                                           if(strcmp(b[j],c)==0)
                                           {
                                                   d[j]++;
                    // fprintf(fp2,"%d\n",j);
                                                   break;
                                           }
                                   }

                  
                                   c[0]='\0';
                               
                                   c[1]='\0';

           for(i=0;i<16;i++)
                   {   d[i]=d[i]/m;
                           fprintf(fp2,"%f\t",d[i]);
                   }
              fprintf(fp2,"\n");
            
                                   D2=0;
                                   fprintf(fp2,"%f\n",D2);
                          
                                  k=0;   
          for(i=0;i<4;i++)
            for(j=0;j<4;j++)
                        {
                          
                                g=((d[k]-e[i]*e[j])*(d[k]-e[i]*e[j]))/(e[i]*e[j]);
               fprintf(fp2,"%f\t",g);
               
                                D2+=g;
               fprintf(fp2,"%f\n",D2);
                                k++;
                        }
                        D2=D2/log(n);

                 fprintf(fp2,"%f\n",D2);
           k=0;
         for(i=0;i<4;i++)
            for(j=0;j<4;j++)
                        {  
            f[k]=(d[k]/(e[i]*e[j]))-1;
                       
             fprintf(fp2,"%f\n",f[k]);
                k++;
                        }

           }
                  
                  
   }
   while(a[0]!='@');
   //for(i=0;i<10;i++)
   //fprintf(fp2,"%ld    ",c[i]);
   fclose(fp1);
   fclose(fp2);
}

zhlong8

1. #!/usr/bin/perl -w
   
2. 
   
3. use strict;
   
4. 
   
5. my %p;
   
6. while (<DATA>) {
   
7.     my($key, $val) = /^P(\w+)\s+(\S+)$/;
   
8.     $p{$key} = $val;
   
9. }
   
10. 
    
11. my $sum = 0;
    
12. for my $ij (glob '{a,t,g,c}{a,t,g,c}') {
    
13.     my($i, $j) = split '', $ij;
    
14.     $sum += ($p{$ij} - $p{$i}*$p{$j})**2 / ($p{$i}*$p{$j});
    
15. }
    
16. 
    
17. print $sum / log(2);
    
18. 
    
19. __DATA__
    
20. Pac        0.091116173
    
21. Pcc        0.054669704
    
22. Pag        0.102505695
    
23. Ptg        0.022779043
    
24. Pat        0.061503417
    
25. Ptc        0.06833713
    
26. Paa        0.141230068
    
27. Pga        0.104783599
    
28. Ptt        0.027334852
    
29. Pct        0.027334852
    
30. Pcg        0.045558087
    
31. Pgg        0.036446469
    
32. Pta        0.029612756
    
33. Pca        0.118451025
    
34. Pgc        0.034168565
    
35. Pgt        0.031890661
    
36. Pa        0.396355353
    
37. Pc        0.248291572
    
38. Pg        0.207289294
    
39. Pt        0.148063781

复制代码

zhagnqiang829

谢谢你的热心帮助。。。回复 5# zhlong8


   

kk861123

zhlong8 发表于 2012-11-24 16:23

第一次看到glob这种用法，学习了

rubyish

一段比较精简的perl6代码实现:

1. #!/usr/bin/perl6
   
2. # 0.105694756594057
   
3. 
   
4. my @data = <
   
5.         Pac        0.091116173
   
6.         Pcc        0.054669704
   
7.         Pag        0.102505695
   
8.         Ptg        0.022779043
   
9.         Pat        0.061503417
   
10.         Ptc        0.06833713
    
11.         Paa        0.141230068
    
12.         Pga        0.104783599
    
13.         Ptt        0.027334852
    
14.         Pct        0.027334852
    
15.         Pcg        0.045558087
    
16.         Pgg        0.036446469
    
17.         Pta        0.029612756
    
18.         Pca        0.118451025
    
19.         Pgc        0.034168565
    
20.         Pgt        0.031890661
    
21.         Pa         0.396355353
    
22.         Pc         0.248291572
    
23.         Pg         0.207289294
    
24.         Pt         0.148063781 >;
    
25. 
    
26. my %p = map { s/^P//; $_ }, @data;
    
27. say ([+] map {(%p{$^a~$^b}-%p{$a}*%p{$b})**2/(%p{$a}*%p{$b})},
    
28.          [X] [<a c g t>] xx 2)/log(2);

复制代码

zhlong8

回复 7# kk861123


    做排列组合那个模块名字想不起来了主要