08

别人写的代码运行真快!!!

最近需要做十几万个差异基因分析,每个分析都是对大约5万个探针,200个样本的数据量进行批量T检验计算P值
然后发现自己无论怎么用R来写,每个分析都要耗时半分钟左右,因为我必须循环所有的探针,即使不用for,而用R推荐的apply系列函数,也快不到哪里去,但是我搜索时候发现有一个package里面自带了矩阵T检验,直接对5万个探针进行T检验,而不需要循环处理它们
看下面代码!
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
差距真的是非常的明显呀!!!
然后,我解析了它的代码,发现里面调用了C写的代码,我想这就是问题所在咯,可是他们到底怎么写,才能把速度搞这么快???
  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")

源码在这个package的github里面可以找到,有兴趣的童鞋可以研究一下

 
 

 

24

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

我想输出ATCG四个字符,组成一个12个字符长度字符串的全排列。共4^12=16777216种排列,

AAAAAAAAAAAA

AAAAAAAAAAAT

AAAAAAAAAAAC

AAAAAAAAAAAG

AAAAAAAAAATA

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

按照正常的想法是通过多重循环来生成全排列,所以我写下了下面这样的代码,但是有个问题,它不支持多扩展性,如果100个全排列,那么得写一百次循环嵌套。

4

所以我求助了perl-china群里的高手,其中一个给了一个结合二进制的巧妙解决方案。

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

这个sprintf是根据024b的格式把我们的十进制数字转为二级制,但是补全为24位。

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这个函数很简单,就是把二进制的数值拆分成字符串数组,两个数字组成一个字符串。

最后通过map把$a{"00"}="A";$a{"01"}="T";$a{"10"}="C";$a{"11"}="G";对应成hash值并且输出即可!

然后这个大神又提成了一个递归的办法来解决

[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]

还有另外一个大神也提成了一个递归的算法解决,算法之道还是蛮有趣的

2

最后还补充一个也是计算机技巧的解决方案,也是群里的朋友提出来的。

[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]

这个是位运算,如果C语言基础学的好的同学很快就能看懂的。

 

28

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

查询需要根据前面建立的索引来做。

BWT算法第二讲532

这是一个比较复杂的过程,我也是看了bowtie的作者的ppt才慢慢弄懂的,感觉自己也不可能三言两语就说清楚,一般都是辅助图片,动画,再经过多方交流才能慢慢理解。

所以大家呢,就自己去看ppt,看懂那个查询算法。(ppt及代码在我的群里面有共享,欢迎大家加群交流)

这里我简单讲讲我的程序

首先读取索引文件,统计好A,C,G,T的总数

然后把查询序列从最后一个字符往前面回溯。

我创建了一个子函数,专门来处理回溯的问题

每次接受四个参数(左右两端的碱基,上下的阈值),并返回两个参数(新的上下两个阈值)

自己动手写bowtie之三序列查询261

大家要看懂阈值是如何更新迭代,这样动态的一个个回溯字符串,一个个迭代阈值。

直到四种临界情况的出现。

自己动手写bowtie之三序列查询477

第一是上下阈值已经相等了,但是我们还没有回溯完全,那就说明字符串只能查找后几个字符,前面还有字符是无法匹配的

第二种情况是上下阈值已经相等了,正巧我们也回溯到了最后一个字符串,那么我们就找到了精确匹配。

第三种情况是已经进行到了最后一个字符串,但是上下阈值还有差值,那么就找到了多个精确匹配点。

最后一种情况是各种非法字符。

然后我简单的测序了一下在病毒的5K基因组里面的精确匹配情况,好像效果还挺好的

自己动手写bowtie之三序列查询518

但是在酵母里面还有一个问题没有解决,就是取前二十个字符串排序的问题,不够精确,需要重新审视排序结果进行局部优化,可能是需要用堆排序发,具体我还得考虑一个星期,只能等下周上课再看看了,平时太忙了,基本没时间码代码。

这里贴上我的代码给大家看看,

[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第二讲:优化索引

BWT算法第二讲192

其中第一讲我提到了一个简单的索引产生方式,因为是课堂就半个小时想的,很多细节没有考虑到,对病毒那种几K大小的基因组来说是很简单的,速度也非常快,但是我测试了一下酵母,却发现好几个小时都没有结果,我只好kill掉重新改写算法,我发现之前的测序最大的问题在于没有立即substr函数的实现方式,把一个5M的字符串不停的截取首尾字符串好像是一个非常慢的方式。

BWT算法第二讲241

所以我优化了那个字符串的函数,虽然代码量变多了,实现过程也繁琐了一点,但是速度提升了几千倍。

 

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

测试了一下我的脚本,对酵母这样的5M的基因组,索引耗费时间是43秒

real 0m43.071s

user 0m41.277s

sys 0m1.779s

输出的index矩阵如下,我简单的截取头尾各10行给大家看,一点要看懂这个index。

BWT算法第二讲532

首先第一列就是我们的BWT向量,也就是BWT变换后的尾字符

第二列是之前的顺序被BWT变换后的首字符排序后的打乱的顺序。

第三,四,五,六列分别是A,C,G,T的计数,就是在当行之前累积出现的A,C,G,T的数量,是对第一列的统计。

 

这个索引文件将会用于下一步的查询,这里贴上我新的索引代码,查询见下一篇文章

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

while  ($F_reads=<FH_F>){

        $R_reads=<FH_R>;

        $F_merge=$residue_F.$F_reads;

        $R_merge=$residue_R.$R_reads;

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

foreach  (0..19) {

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

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

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

                        $i++;

}

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

$residue_F=$F_reads;

$residue_R=$R_reads;

}

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

}

unlink("tmp_forward.txt");

unlink("tmp_reverse.txt");

[/perl]