18

# Bowtie算法第六讲-tally法对bwt索引进行搜索

perl rotation_one_by_one.pl atgtgtcgtagctcgtnncgt

\$ATGTGTCGTAGCTCGTNNCGT 21

AGCTCGTNNCGT\$ATGTGTCGT 9

ATGTGTCGTAGCTCGTNNCGT\$ 0

CGT\$ATGTGTCGTAGCTCGTNN 18

CGTAGCTCGTNNCGT\$ATGTGT 6

CGTNNCGT\$ATGTGTCGTAGCT 13

CTCGTNNCGT\$ATGTGTCGTAG 11

GCTCGTNNCGT\$ATGTGTCGTA 10

GT\$ATGTGTCGTAGCTCGTNNC 19

GTAGCTCGTNNCGT\$ATGTGTC 7

GTCGTAGCTCGTNNCGT\$ATGT 4

GTGTCGTAGCTCGTNNCGT\$AT 2

GTNNCGT\$ATGTGTCGTAGCTC 14

NCGT\$ATGTGTCGTAGCTCGTN 17

NNCGT\$ATGTGTCGTAGCTCGT 16

T\$ATGTGTCGTAGCTCGTNNCG 20

TAGCTCGTNNCGT\$ATGTGTCG 8

TCGTAGCTCGTNNCGT\$ATGTG 5

TCGTNNCGT\$ATGTGTCGTAGC 12

TGTCGTAGCTCGTNNCGT\$ATG 3

TGTGTCGTAGCTCGTNNCGT\$A 1

TNNCGT\$ATGTGTCGTAGCTCG 15

T 21

T 9

\$ 0

N 18

T 6

T 13

G 11

A 10

C 19

C 7

T 4

T 2

C 14

N 17

T 16

G 20

G 8

G 5

C 12

G 3

A 1

G 15 perl search_char.pl    GTGTCG   tm.tally

your last char is G

start is 7 ; and end is 13

now it is number 5 and the char is C

start is 3 ; and end is 6

now it is number 4 and the char is T

start is 17 ; and end is 19

now it is number 3 and the char is G

start is 10 ; and end is 11

now it is number 2 and the char is T

start is 19 ; and end is 20

now it is number 1 and the char is G

start is 11 ; and end is 12

It is just one perfect match !

The index is 2

perl search_char.pl  GT  tm.tally

your last char is T

start is 15 ; and end is 22

now it is number 1 and the char is G

start is 8 ; and end is 13

we find more than one perfect match!!!

8 13

One of the index is 11

One of the index is 10

One of the index is 19

One of the index is 7

One of the index is 4

One of the index is 2

One of the index is 14

CTCGTNNCGT\$ATGTGTCGTAG 11

GCTCGTNNCGT\$ATGTGTCGTA 10

GT\$ATGTGTCGTAGCTCGTNNC 19

GTAGCTCGTNNCGT\$ATGTGTC 7

GTCGTAGCTCGTNNCGT\$ATGT 4

GTGTCGTAGCTCGTNNCGT\$AT 2

GTNNCGT\$ATGTGTCGTAGCTC 14

perl search_char.pl ACATGTGT tm.tally

your last char is T

start is 15 ; and end is 22

now it is number 7 and the char is G

start is 8 ; and end is 13

now it is number 6 and the char is T

start is 19 ; and end is 21

now it is number 5 and the char is G

start is 11 ; and end is 12

now it is number 4 and the char is T

start is 20 ; and end is 21

now it is number 3 and the char is A

start is 2 ; and end is 3

now it is number 2 and the char is C

start is 3 ; and end is 3

we can just find the last 6 char ,and it is ATGTGT

[perl]

\$a=\$ARGV;

\$a=uc \$a;

open FH,"<\$ARGV";

while(<FH>){

chomp;

@F=split;

\$hash_count_atcg{\$F}++;

\$hash{\$.}=\$_;

# the first line is \$ and the last char and the last index !

}

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

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

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

\$all_n=\$hash_count_atcg{'N'};

\$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";

\$after=substr(\$a,\$end_a,1);

#we fill search your query from the last char !

if (\$after eq 'A') {

\$start=2;

\$end=\$all_a+1;

}

elsif (\$after eq 'C') {

\$start=\$all_a+1;

\$end=\$all_a+\$all_c+1;

}

elsif (\$after eq 'G') {

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

\$end=\$all_a+\$all_c+\$all_g+1;

}

elsif (\$after eq 'T'){

\$start=\$all_a+\$all_c+\$all_g+\$all_n+1;

\$end=\$all_a+\$all_c+\$all_g+\$all_t+\$all_n+1;

}

else {die "error !!! we just need A T C G !!!\n"}

print "your last char is \$after\n ";

print "start is \$start ; and end is \$end \n";

foreach (reverse (1..\$end_a)){

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

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

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

print "now it is number \$_ and the char is \$before \n ";

print "start is \$start ; and end is \$end \n";

if (\$_  > 1 && \$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";

#return "miss";

last;

}

if (\$_ == 1) {

if ((\$end-\$start)==1) {

print "It is just one perfect match ! \n";

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

print "The index is \$F_start\n";

#return \$F_start;

last;

}

else {

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

print "\$start\t\$end\n";

foreach  ((\$start-1)..\$end) {

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

print "One of the index is \$F_start\n";

}

#return "multiple";

last;

}

}

}

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+1,\$F_end+1);

}

elsif (\$before eq 'C') {

return (\$all_a+\$F_start+1,\$all_a+\$F_end+1);

}

elsif (\$before eq 'G') {

return (\$all_a+\$all_c+1+\$F_start,\$all_a+\$all_c+1+\$F_end);

}

elsif (\$before eq 'T') {

return (\$all_a+\$all_c+\$all_g+\$all_n+1+\$F_start,\$all_a+\$all_c+\$all_g+1+\$all_n+\$F_end);

}

else {die "error !!! we just need A T C G !!!\n"}

}

[/perl]

18

# Bowtie 算法第四讲

[perl]

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

open FH ,"<\$ARGV";

\$i=0;

while (<FH>) {

next if /^>/;

chomp;

\$a.=(uc);

}

#print "\$a\n";

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

\$query=uc \$ARGV;

\$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]

[perl]

\$a=uc \$ARGV;

\$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

[perl]

#first read the tally !!!

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

open FH,"<\$ARGV";

\$hash_count{'A'}=0;

\$hash_count{'C'}=0;

\$hash_count{'G'}=0;

\$hash_count{'T'}=0;

while(<FH>){

chomp;

@F=split;

\$hash_count{\$F}++;

\$hash{\$.}="\$F\t\$F\t\$hash_count{\$F}";

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

my \$before=\$F;

if (\$before eq 'A') {

\$new=\$F+1;

}

elsif (\$before eq 'C') {

\$new=1+\$all_a+\$F;

}

elsif (\$before eq 'G') {

\$new=1+\$all_a+\$all_c+\$F;

}

elsif (\$before eq 'N') {

\$new =1+\$all_a+\$all_c+\$all_g+\$F;

}

elsif (\$before eq 'T') {

\$new=1+\$all_a+\$all_c+\$all_g+\$all_n+\$F;

}

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\t\$new\n";

&restore(\$new);

}

[/perl]

28

# 自己动手写bowtie第三讲：序列查询。    [perl]

\$a='CGCTATGTACTGGATGCGCTGGCAAACGAGCCTGCCGTAAG';

while(<>){

chomp;

@F=split;

\$hash_count_atcg{\$F}++;

\$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\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\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,\$F_end);

}

elsif (\$before eq 'C') {

return (\$all_a+\$F_start,\$all_a+\$F_end);

}

elsif (\$before eq 'G') {

return (\$all_a+\$all_c+\$F_start,\$all_a+\$all_c+\$F_end);

}

elsif (\$before eq 'T') {

return (\$all_a+\$all_c+\$all_g+\$F_start,\$all_a+\$all_c+\$all_g+\$F_end);

}

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]

23

# Snp-calling流程（BWA+SAMTOOLS+BCFTOOLS）

samtools faidx hg19.fa

Bowtie2-build hg19.fa hg19

Bowtie的使用方法详解见http://www.bio-info-trainee.com/?p=398

bowtie2 -x ../../../ref-database/hg19 -U  tmp1.fa -S tmp1.sam

samtools view -bS tmp1.sam > tmp1.bam

samtools sort tmp1.bam tmp1.sorted

samtools index tmp1.sorted.bam

samtools mpileup -d 1000  -gSDf   ../../../ref-database/hg19.fa  tmp1.sorted.bam |bcftools view -cvNg -  >tmp1.vcf

==> bwa_aln.sh <==

do

echo \$id

bwa aln hg19.fa \$id >\$id.sai

done <\$1

==> bwa_sampe.sh <==

do

echo \$id

bwa sampe hg19.fa \$id*sai \$id*single >\$id.sam

done <\$1

==> samtools.sh <==

do

echo \$id

samtools view -bS \$id.sam > \$id.bam

samtools sort \$id.bam \$id.sorted

samtools index \$id.sorted.bam

done <\$1

==> bcftools.sh <==

do

echo \$id

samtools mpileup -d 1000  -gSDf  ref.fa \$id*sorted.bam |bcftools view -cvNg -  >\$id.vcf

done <\$1

==> mpileup.sh <==

do

echo \$id

samtools mpileup -d 100000 -f hg19.fa \$id*sorted.bam >\$id.mpileup

done <\$1

20

# 自己动手写bowtie第一讲：BWT算法详解并建立索引

http://www.cnblogs.com/xudong-bupt/p/3763814.html

A1A2A3...An

A2A3...AnA1

A3...AnA1A2

...

AnA1A2A3...

n个串，每个字符串的长度都是n。  Perl程序如下

[perl]

while (<>){

next if />/;

chomp;

\$a.=\$_;

}

\$a.='\$';

\$len=length \$a;

\$i=0;

print "first we transform it !!!\n";

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

\$up=substr(\$a,0,\$_);

\$down=substr(\$a,\$_);

#print "\$down\$up\n";

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

\$key=substr("\$down\$up",0,20);

\$key=\$key.”\t”.substr("\$down\$up",\$len-1);

\$hash{\$key}=\$i;

\$i++;

}

print "then we sort it\n";

foreach  (sort keys  %hash){

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

\$len=length;

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

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

print "\$_\t\$hash{\$_}\n";

}

[/perl] http://tieba.baidu.com/p/1504205984

http://www.cnblogs.com/xudong-bupt/p/3763814.html

20

# bowtie简单使用

http://bowtie-bio.sourceforge.net/bowtie2/index.shtml   • 索引，bowtie2-build your-fastq-file.fa your-index-name • 然后比对，分两种，一是单端测序数据，二是双端数据

Bowtie –x tmp –U reads.fa –S hahahhha.sam

Bowtie –x tmp -1 reads1.fa -2 reads2.fa –S hahahha.sam

14