自己动手写bowtie第二讲：优化索引 | 生信菜鸟团

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

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

 

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向量，也就是BWT变换后的尾字符

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

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

 

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

查看源代码

打印帮助

1	while (<>){

2

3	next if />/;

4

5

chomp;

6

7

$a.=$_;

8

9

}

10

11	$len=length $a;

12

13	open FH_F,'>tmp_forward.txt';

14

15	open FH_R,'>tmp_reverse.txt';

16

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

18

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

20

21	print FH_F '\n';

22

23

}

24

25	$rev_a=reverse $a;

26

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

28

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

30

31	print FH_R '\n';

32

33

}

34

35	close FH_F;

36

37	close FH_R;

38

39

$a='';

40

41	open FH_F,'tmp_forward.txt';

42

43	open FH_R,'tmp_reverse.txt';

44

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

46

47	$residue_F=;

48

49	$residue_R=;

50

51

$i=0;

52

53	while  ($F_reads=){

54

55	$R_reads=;

56

57	$F_merge=$residue_F.$F_reads;

58

59	$R_merge=$residue_R.$R_reads;

60

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

62

63	foreach  (0..19) {

64

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

66

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

68

69	$hash{'$up\t$down'}=$i;

70

71

$i++;

72

73

}

74

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

76

77	$residue_F=$F_reads;

78

79	$residue_R=$R_reads;

80

81

}

82

83	#print 'then we sort it\n';

84

85	$count_a=0;

86

87	$count_c=0;

88

89	$count_g=0;

90

91	$count_t=0;

92

93	foreach  (sort keys  %hash){

94

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

96

97	$len=length;

98

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

100

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

102

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

104

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

106

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

108

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

110

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

112

113

}

114

115	unlink('tmp_forward.txt');

116