自己动手写bowtie第三讲：序列查询。 | 生信菜鸟团

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

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

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

这里我简单讲讲我的程序

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

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

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

每次接受四个参数（左右两端的碱基，上下的阈值），并返回两个参数（新的上下两个阈值）

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

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

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

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

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

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

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

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

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

查看源代码

打印帮助

1	$a='CGCTATGTACTGGATGCGCTGGCAAACGAGCCTGCCGTAAG';

2

3	while(<>){

4

5

chomp;

6

7

@F=split;

8

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

10

11	$hash{$.}=$_;

12

13

}

14

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

16

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

18

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

20

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

22

23	#print '$all_a\t$all_c\t$all_g\t$all_t\n';

24

25	$len_a=length $a;

26

27	$end_a=$len_a-1;

28

29	print 'your query is $a\n';

30

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

32

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

34

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

36

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

38

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

40

41	if($_ == $end_a){

42

43	if ($after eq 'A') {

44

45

$start=1;

46

47	$end=$all_a;

48

49

}

50

51	elsif ($after eq 'C') {

52

53	$start=$all_a+1;

54

55	$end=$all_a+$all_c;

56

57

}

58

59	elsif ($after eq 'G') {

60

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

62

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

64

65

}

66

67	elsif ($after eq 'T'){

68

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

70

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

72

73

}

74

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

76

77

}

78

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

80

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

82

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

84

85	$find_len=length $find_char;

86

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

88

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

90

91

exit;

92

93

}

94

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

96

97	if ($_ == 0) {

98

99	if ($start == $end) {

100

101	print 'It is just one perfect match ! \n';

102

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

104

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

106

107

exit;

108

109

}

110

111

else {

112

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

114

115	#print '$start\t$end\n';

116

117	foreach  ($start..$end) {

118

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

120

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

122

123

}

124

125

exit;

126

127

}

128

129

}

130

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

132

133

}

134

135	sub find_level{

136

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

138

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

140

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

142

143	if ($before eq 'A') {

144

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

146

147

}

148

149	elsif ($before eq 'C') {

150

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

152

153

}

154

155	elsif ($before eq 'G') {

156

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

158

159

}

160

161	elsif ($before eq 'T') {

162

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

164

165

}

166

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

168

169

}

170

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