由于之前就简单的看了看bowtie作者的ppt，没有完全吃透就开始敲代码了，写了十几个程序最后我自己都搞不清楚进展到哪一步了，所以我现在整理一下，从新开始！！！
首先，bowtie的作用就是在一个大字符串里面搜索一个小字符串！那么本身就有一个非常笨的复杂方法来搜索，比如，大字符串长度为100万，小字符串为10，那么就依次取出大字符串的10个字符来跟小字符串比较即可，这样的算法是非常不经济的，我简单用perl代码实现一下。
[perl]
#首先读取大字符串的fasta文件
open FH ,"<$ARGV[0]";
$i=0;
while (<FH>) {
next if /^>/;
chomp;
$a.=(uc);
}
#print "$a\n";
#然后接受我们的小的查询字符串
$query=uc $ARGV[1];
$len=length $a;
$len_query=length $query;
$a=$a.'$'.$a;
#然后依次循环取大字符串来精确比较！
foreach (0..$len-1){
if (substr($a,$_,$len_query) eq $query){
print "$_\n";
#last;
}
}
[/perl]
这样在时间复杂度非常恐怖，尤其是对人的30亿碱基。
正是因为这样的查询效率非常低，所以我们才需要用bwt算法来构建索引，然后根据tally来进行查询
其中构建索引有三种方式，我首先讲最效率最低的那种索引构造算法，就是依次取字符串进行旋转，然后排序即可。
[perl]
$a=uc $ARGV[0];
$len=length $a;
$a=$a.'$'.$a;
foreach (0..$len){
$hash{substr($a,$_,$len+1)}=$_;
}
#print "$_\t$hash{$_}\n" foreach sort keys %hash;
print substr($_,-1),"\t$hash{$_}\n" foreach sort keys %hash;
[/perl]
这个算法从时间复杂度来讲是非常经济的，对小字符串都是瞬间搞定！！！
perl rotation_one_by_one.pl atgcgtanngtc 这个字符串的BWT矩阵索引如下！
C 12
T 6
$ 0
T 11
G 3
T 2
C 4
N 9
N 8
A 7
G 5
G 10
A 1
但同样的，它也有一个无法避免的弊端，就是内存消耗太恐怖。对于30亿的人类碱基来说，这样旋转会生成30亿乘以30亿的大矩阵，一般的服务器根本hold不住的。
最后我讲一下，这个BWT矩阵索引如何还原成原字符串，这个没有算法的差别，因为就是很简单的原理。
[perl]
#first read the tally !!!
#首先读取上面输出的BWT矩阵索引文件。
open FH,"<$ARGV[0]";
$hash_count{'A'}=0;
$hash_count{'C'}=0;
$hash_count{'G'}=0;
$hash_count{'T'}=0;
while(<FH>){
chomp;
@F=split;
$hash_count{$F[0]}++;
$hash{$.}="$F[0]\t$F[1]\t$hash_count{$F[0]}";
#print "$hash{$.}\n";
}
$all_a=$hash_count{'A'};
$all_c=$hash_count{'C'};
$all_g=$hash_count{'G'};
$all_t=$hash_count{'T'};
$all_n=$hash_count{'N'};
#start from the first char !
$raw='';
&restore(1);
sub restore{
my($num)=@_;
my @F=split/\t/,$hash{$num};
$raw.=$F[0];
my $before=$F[0];
if ($before eq 'A') {
$new=$F[2]+1;
}
elsif ($before eq 'C') {
$new=1+$all_a+$F[2];
}
elsif ($before eq 'G') {
$new=1+$all_a+$all_c+$F[2];
}
elsif ($before eq 'N') {
$new =1+$all_a+$all_c+$all_g+$F[2];
}
elsif ($before eq 'T') {
$new=1+$all_a+$all_c+$all_g+$all_n+$F[2];
}
elsif ($before eq '$') {
chop $raw;
$raw = reverse $raw;
print "$raw\n";
exit;
}
else {die "error !!! we just need A T C N G !!!\n"}
#print "$F[0]\t$new\n";
&restore($new);
}
[/perl]