08

# 别人写的代码运行真快！！！

dat=matrix(rnorm(10000000),nrow = 50000)
dim(dat) #50000   200
system.time(
apply(dat,1,function(x){
t.test(x[1:100],x[101:200])\$p.value
})
)
#用户  系统  流逝
#29.29  0.04 30.64
library(pi0)
system.time(matrix.t.test(dat,1,100,100))
#用户 系统 流逝
#0.48 0.03 0.53

tmp = .C("tstatistic", dat = x, n1 = n1, n2 = n2, ntests = ntests,
MARGIN = MARGIN, pool = pool, tstat = rep(0, ntests),
df = rep(0, ntests), PACKAGE = "pi0")

24

# perl多重循环生成全排列精简

AAAAAAAAAAAA

AAAAAAAAAAAT

AAAAAAAAAAAC

AAAAAAAAAAAG

AAAAAAAAAATA

`````````````````````````

perl -le '{\$a{"00"}="A";\$a{"01"}="T";\$a{"10"}="C";\$a{"11"}="G";for(0..100){print join"",map{\$a{\$_}} unpack("a2"x12,sprintf("%024b\n",\$_))}}'

AAAAAAAAAAAA

AAAAAAAAAAAT

AAAAAAAAAAAC

AAAAAAAAAAAG

AAAAAAAAAATA

AAAAAAAAAATT

AAAAAAAAAATC

AAAAAAAAAATG

AAAAAAAAAACA

AAAAAAAAAACT

AAAAAAAAAACC

AAAAAAAAAACG

AAAAAAAAAAGA

AAAAAAAAAAGT

AAAAAAAAAAGC

AAAAAAAAAAGG

AAAAAAAAATAA

AAAAAAAAATAT

perl -le '{\$a=sprintf("%024b\n",1);print \$a}'

000000000000000000000001

perl -le '{@a=unpack("a2"x12,sprintf("%024b\n",100));print join":", @a}'

00:00:00:00:00:00:00:00:01:10:01:00

unpack这个函数很简单，就是把二进制的数值拆分成字符串数组，两个数字组成一个字符串。

[perl]

#!/usr/bin/perl

my \$arr = [];

\$arr->[\$_] = 0 for (0..11);
my @b = ("A","C","G","T");

while(1){
print join "", map {\$b[\$_]} @\$arr;
print "\n";
exit unless &count(\$arr,11);
}

sub count {
my \$arr = shift;
my \$x = shift;

return 0 if \$x < 0;

if (\$arr->[\$x] < 3){
\$arr->[\$x]++;
return 1;
}

else{
\$arr->[\$x] = 0;
return &count(\$arr,\$x - 1);
}
}

[/perl]

[perl]

#!/usr/bin/perl -w

my @A = qw(A T G C);
my \$N = 12;

foreach my \$n (0..4**\$N){
foreach my \$index (0..\$N){
print \$A[\$n%4];
\$n = \$n >>2;
}
print "\n";
}

[/perl]

28

# 自己动手写bowtie第三讲：序列查询。

[perl]

\$a='CGCTATGTACTGGATGCGCTGGCAAACGAGCCTGCCGTAAG';

while(<>){

chomp;

@F=split;

\$hash_count_atcg{\$F[0]}++;

\$hash{\$.}=\$_;

}

\$all_a=\$hash_count_atcg{'A'};

\$all_c=\$hash_count_atcg{'C'};

\$all_g=\$hash_count_atcg{'G'};

\$all_t=\$hash_count_atcg{'T'};

#print "\$all_a\t\$all_c\t\$all_g\t\$all_t\n";

\$len_a=length \$a;

\$end_a=\$len_a-1;

print "your query is \$a\n";

print "and the length of your query is \$len_a \n";

foreach (reverse (0..\$end_a)){

\$after=substr(\$a,\$_,1);

\$before=substr(\$a,\$_-1,1);

#对第一个字符进行找阈值的时候，我们需要人为的定义起始点！

if(\$_ == \$end_a){

if (\$after eq 'A') {

\$start=1;

\$end=\$all_a;

}

elsif (\$after eq 'C') {

\$start=\$all_a+1;

\$end=\$all_a+\$all_c;

}

elsif (\$after eq 'G') {

\$start=\$all_a+\$all_c+1;

\$end=\$all_a+\$all_c+\$all_g;

}

elsif (\$after eq 'T'){

\$start=\$all_a+\$all_c+\$all_g+1;

\$end=\$all_a+\$all_c+\$all_g+\$all_t;

}

else {print "error !!! we just need A T C G !!!\n";exit;}

}

#如果阈值已经无法继续分割，但是字符串还未查询完

if (\$_  > 0 && \$start == \$end) {

\$find_char=substr(\$a,\$_);

\$find_len=length \$find_char;

#这里需要修改，但是不影响完全匹配了

print "we can just find the last \$find_len char ,and it is \$find_char \n";

exit;

}

#如果进行到了最后一个字符

if (\$_ == 0) {

if (\$start == \$end) {

print "It is just one perfect match ! \n";

my @F_start=split/\s+/,\$hash{\$start};

print "The index is \$F_start[1]\n";

exit;

}

else {

print "we find more than one perfect match!!!\n";

#print "\$start\t\$end\n";

foreach  (\$start..\$end) {

my @F_start=split/\s+/,\$hash{\$_};

print "One of the index is \$F_start[1]\n";

}

exit;

}

}

(\$start,\$end)=&find_level(\$after,\$before,\$start,\$end);

}

sub find_level{

my(\$after,\$before,\$start,\$end)=@_;

my @F_start=split/\s+/,\$hash{\$start};

my @F_end=split/\s+/,\$hash{\$end};

if (\$before eq 'A') {

return (\$F_start[2],\$F_end[2]);

}

elsif (\$before eq 'C') {

return (\$all_a+\$F_start[3],\$all_a+\$F_end[3]);

}

elsif (\$before eq 'G') {

return (\$all_a+\$all_c+\$F_start[4],\$all_a+\$all_c+\$F_end[4]);

}

elsif (\$before eq 'T') {

return (\$all_a+\$all_c+\$all_g+\$F_start[5],\$all_a+\$all_c+\$all_g+\$F_end[5]);

}

else {print "sorry , I can't find the right match!!!\n";}

}

#perl -alne '{next if />/;\$all.=\$_;}END{print substr(\$all,308,10)}'   lambda_virus.fa

[/perl]

28

# 自己动手写bowtie第二讲：优化索引

time perl bwt_new_index.pl e-coli.fa >e-coli.index

real 0m43.071s

user 0m41.277s

sys 0m1.779s

[perl]

while (<>){

next if />/;

chomp;

\$a.=\$_;

}

\$len=length \$a;

open FH_F,">tmp_forward.txt";

open FH_R,">tmp_reverse.txt";

for(my \$i=0;\$i<=\$len-1;\$i+=20){

print FH_F  substr(\$a,\$i,20);

print FH_F "\n";

}

\$rev_a=reverse \$a;

for(my \$i=0;\$i<=\$len-1;\$i+=20){

print FH_R  substr(\$rev_a,\$i,20);

print FH_R "\n";

}

close FH_F;

close FH_R;

\$a='';

open FH_F,"tmp_forward.txt";

open FH_R,"tmp_reverse.txt";

#把前一行的所有20bp碱基当做后一行的头部信息

\$residue_F=<FH_F>;

\$residue_R=<FH_R>;

\$i=0;

#这样每次就需要处理20个碱基

foreach  (0..19) {

\$up  =substr(\$F_merge,\$_,20);

\$down=substr(\$R_merge,\$_,1);

\$hash{"\$up\t\$down"}=\$i;

\$i++;

}

#处理完毕之后再保存当行的20bp碱基做下一行的头部信息

}

#print "then we sort it\n";

\$count_a=0;

\$count_c=0;

\$count_g=0;

\$count_t=0;

foreach  (sort keys  %hash){

\$first=substr(\$_,0,1);

\$len=length;

\$last=substr(\$_,\$len-1,1);

#print "\$first\t\$last\t\$hash{\$_}\n";

\$count_a++ if \$last eq 'A';

\$count_c++ if \$last eq 'C';

\$count_g++ if \$last eq 'G';

\$count_t++ if \$last eq 'T';

print "\$last\t\$hash{\$_}\t\$count_a\t\$count_c\t\$count_g\t\$count_t\n";

}