还是老问题,请各位帮忙!
随机序列.rar
(1.39 MB)
这是500万长度的随机序列,由A、T、C、G四个字母组成。
八模体及其出现次数.zip
(233.39 KB)
这是对应附件得到的结果模式,4的8次幂个八模体及其各自的出现次数。8模体由8个碱基构成,(碱基即A、T、C、G)每个碱基有4种选择(或A、或T、或C、或G),这样就有4的8次幂种8模体(即65536种)。
附件1中的随机序列包含500万个字符,如果以1为间隔的话,应有 4999993个 8模体,因为每[i,i+7]位置区域上的8个碱基长的字符串都是一个8模体,i从0开始,直到4999992结束,i是i++逐渐增加的,附件2的结果就是在以1为间隔下得到的。 4999993个8模体由65536种8模体组成,附件2中的结果就是在4999993个8模体中每种8模体出现了多少次。
如果以2为间隔的话,应有 2499996个 8模体,因为每[i,i+7]位置区域上的8个碱基长的字符串都是一个8模体,i从0开始,直到4999992结束,i是i+=2逐渐增加的。
如果以k为间隔的话,应有 500万除k个 8模体,因为每[i,i+7]位置区域上的8个碱基长的字符串都是一个8模体,i从0开始,i是i+=k逐渐增加的,直到i+7>4999999结束。
现在我能计算以1为间隔的,但以k为间隔的程序实在弄不出来。我是学物理的学生,计算机这块有些棘手,麻烦各位了。
谢谢。
这是统计以1为间隔的程序。
#include <stdio.h>
#define N_A 9
#define N_B (1u<<(2*N_A))
int main()
{
static unsigned long octamer_count[N_B] = { 0 }; //当跑八模体以上时,因为数据量已经很大,要价个static,八模体以下不用它
////// 1 //////
{
FILE* file = fopen( "E:\\酵母\\光学与信息学\\碱基随机.txt", "rb" );
if( !file )
{
puts( "cannot open the input file.");
return 1;
}
{
unsigned long octamer;
unsigned long octamer_n = 0;
char c;
for( ; fread(&c,1,1,file)==1; )
{
switch( c )
{
case 'A':
octamer = (octamer*4+0)%N_B; ++octamer_n; break;
case 'T':
octamer = (octamer*4+1)%N_B; ++octamer_n; break;
case 'C':
octamer = (octamer*4+2)%N_B; ++octamer_n; break;
case 'G':
octamer = (octamer*4+3)%N_B; ++octamer_n; break;
case '\r':
case '\n':
continue;
default:
puts( "Fuck" );
case 'N':
octamer_n = 0;
continue;
}
if( octamer_n < N_A )
continue;
octamer_n = N_A;
++octamer_count[octamer];
}
}
fclose( file );
}
////// 2 //////
{
FILE* fileout = fopen( "E:\\酵母\\光学与信息学\\碱基随机9模体.txt", "wt" );
if( !fileout )
{
puts( "cannot open the output file.");
return 2;
}
{
unsigned int i, j;
for( i=0; i<sizeof(octamer_count)/sizeof(octamer_count[0]); ++i )
{
if( octamer_count[i] == 0 )
continue;
for( j=N_B/4; j!=0; j/=4 )
fprintf( fileout, "%c", "ATCG"[i/j%4] );
fprintf( fileout, "\t%u\n", octamer_count[i] );
}
}
fclose( fileout );
}
return 0;
}