Bowtie?算法第四讲 | 生信菜鸟团

首先，bowtie的作用就是在一个大字符串里面搜索一个小字符串！那么本身就有一个非常笨的复杂方法来搜索，比如，大字符串长度为100万，小字符串为10，那么就依次取出大字符串的10个字符来跟小字符串比较即可，这样的算法是非常不经济的，我简单用perl代码实现一下。

查看源代码

打印帮助

1	#首先读取大字符串的fasta文件

2

3	open FH ,'<$ARGV[0]';

4

5

$i=0;

6

7	while () {

8

9	next if /^>/;

10

11

chomp;

12

13

$a.=(uc);

14

15

}

16

17	#print '$a\n';

18

19	#然后接受我们的小的查询字符串

20

21	$query=uc $ARGV[1];

22

23	$len=length $a;

24

25	$len_query=length $query;

26

27	$a=$a.'$'.$a;

28

29	#然后依次循环取大字符串来精确比较！

30

31	foreach (0..$len-1){

32

33	if (substr($a,$_,$len_query) eq $query){

34

35	print '$_\n';

36

37

#last;

38

39

}

40

41

}

 

这样在时间复杂度非常恐怖，尤其是对人的30亿碱基。

 

正是因为这样的查询效率非常低，所以我们才需要用bwt算法来构建索引，然后根据tally来进行查询

其中构建索引有三种方式，我首先讲最效率最低的那种索引构造算法，就是依次取字符串进行旋转，然后排序即可。

查看源代码

打印帮助

1	$a=uc $ARGV[0];

2

3	$len=length $a;

4

5	$a=$a.'$'.$a;

6

7	foreach (0..$len){

8

9	$hash{substr($a,$_,$len+1)}=$_;

10

11

}

12

13	#print '$_\t$hash{$_}\n' foreach sort keys %hash;

14

15	print  substr($_,-1),'\t$hash{$_}\n' foreach sort keys %hash;

这个算法从时间复杂度来讲是非常经济的，对小字符串都是瞬间搞定！！！

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矩阵索引如何还原成原字符串，这个没有算法的差别，因为就是很简单的原理。

查看源代码

打印帮助

1	#first read the tally !!!

2

3	#首先读取上面输出的BWT矩阵索引文件。

4

5	open FH,'<$ARGV[0]';

6

7	$hash_count{'A'}=0;

8

9	$hash_count{'C'}=0;

10

11	$hash_count{'G'}=0;

12

13	$hash_count{'T'}=0;

14

15

while(){

16

17

chomp;

18

19

@F=split;

20

21	$hash_count{$F[0]}++;

22

23	$hash{$.}='$F[0]\t$F[1]\t$hash_count{$F[0]}';

24

25	#print '$hash{$.}\n';

26

27

28

29

}

30

31	$all_a=$hash_count{'A'};

32

33	$all_c=$hash_count{'C'};

34

35	$all_g=$hash_count{'G'};

36

37	$all_t=$hash_count{'T'};

38

39	$all_n=$hash_count{'N'};

40

41	#start from the first char !

42

43

$raw='';

44

45	&restore(1);

46

47	sub restore{

48

49	my($num)=@_;

50

51	my @F=split/\t/,$hash{$num};

52

53	$raw.=$F[0];

54

55	my $before=$F[0];

56

57	if ($before eq 'A') {

58

59	$new=$F[2]+1;

60

61

}

62

63	elsif ($before eq 'C') {

64

65	$new=1+$all_a+$F[2];

66

67

}

68

69	elsif ($before eq 'G') {

70

71	$new=1+$all_a+$all_c+$F[2];

72

73

}

74

75	elsif ($before eq 'N') {

76

77	$new =1+$all_a+$all_c+$all_g+$F[2];

78

79

}

80

81	elsif ($before eq 'T') {

82

83	$new=1+$all_a+$all_c+$all_g+$all_n+$F[2];

84

85

}

86

87	elsif ($before eq '$') {

88

89	chop $raw;

90

91	$raw = reverse $raw;

92

93	print '$raw\n';

94

95

exit;

96

97

}

98

99	else {die 'error !!! we just need A T C N G !!!\n'}

100

101	#print '$F[0]\t$new\n';

102

103	&restore($new);

104

105

}