亚洲av成人无遮挡网站在线观看,少妇性bbb搡bbb爽爽爽,亚洲av日韩精品久久久久久,兔费看少妇性l交大片免费,无码少妇一区二区三区

Chinaunix

標(biāo)題: 用程序怎么計算下面公式 [打印本頁]
作者: zhagnqiang829    時間: 2012-11-23 20:24
標(biāo)題: 用程序怎么計算下面公式
怎么計算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    時間: 2012-11-23 23:55
這不就是把這個 ∑ 展開么?
如果我的數(shù)學(xué)老師教得沒錯的話, 不就是個 加法么? 兩層的循環(huán)加一下好啦. 結(jié)果再除以 ln(2) 不就 OK了?
作者: zhagnqiang829    時間: 2012-11-24 08:41
你幫我改下代碼吧，我運(yùn)行總出錯。
#!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    時間: 2012-11-24 12:01
唉，弄不出來。我還是那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    時間: 2012-11-24 16:23
  1. #!/usr/bin/perl -w

  3. use strict;

  5. my %p;
  6. while (<DATA>) {
  7.     my($key, $val) = /^P(\w+)\s+(\S+)$/;
  8.     $p{$key} = $val;
  9. }

  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. }

  17. print $sum / log(2);

  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    時間: 2012-11-24 18:04
謝謝你的熱心幫助。。。回復(fù) 5# zhlong8


   
作者: kk861123    時間: 2012-11-25 09:36
zhlong8 發(fā)表于 2012-11-24 16:23

第一次看到glob這種用法，學(xué)習(xí)了
作者: rubyish    時間: 2012-11-25 14:39
一段比較精簡的perl6代碼實(shí)現(xiàn):
  1. #!/usr/bin/perl6
  2. # 0.105694756594057

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

  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    時間: 2012-11-25 19:17
回復(fù) 7# kk861123


    做排列組合那個模塊名字想不起來了主要



歡迎光臨 Chinaunix (http://www.72891.cn/) Powered by Discuz! X3.2