用程序怎么計(jì)算下面公式

zhagnqiang829

怎么計(jì)算D2…
    ，（i,j =a,c,g,t）
其中 1/ln(2) = 1.442695040889
一  下面是所有出現(xiàn)的概率的情況( 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

這不就是把這個(gè) ∑ 展開(kāi)么?
如果我的數(shù)學(xué)老師教得沒(méi)錯(cuò)的話, 不就是個(gè) 加法么? 兩層的循環(huán)加一下好啦. 結(jié)果再除以 ln(2) 不就 OK了?

zhagnqiang829

你幫我改下代碼吧，我運(yùn)行總出錯(cuò)。
#!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];
}
}


回復(fù) 2# q1208c


   

zhagnqiang829

唉，弄不出來(lái)。我還是那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:\\計(jì)算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

復(fù)制代碼

zhagnqiang829

謝謝你的熱心幫助。。。回復(fù) 5# zhlong8


   

kk861123

zhlong8 發(fā)表于 2012-11-24 16:23

第一次看到glob這種用法，學(xué)習(xí)了

rubyish

一段比較精簡(jiǎn)的perl6代碼實(shí)現(xiàn):

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);

復(fù)制代碼

zhlong8

回復(fù) 7# kk861123


    做排列組合那個(gè)模塊名字想不起來(lái)了主要