使用perl 查找str序列
没啥说的,直接上代码
#!/usr/bin/perl -w use strict; use Bio::SeqIO; #use Bio::Tools::TandemRepeatsFinder; #Bio perl 需要事先安装 #安装过程如下,以ubuntu系统以为 # 1、安装cpanm # wget http://xrl.us/cpanm -O /usr/bin/cpanm; # chmod +x /usr/bin/cpanm # 2、安装Bio::perl # cpanm --mirror http://mirrors.163.com/cpan --mirror-only Bio::Perl sub found_str{ my %STRS; my $seqio_obj = Bio::SeqIO->new(-file => $_[0], -format => "fasta" ); while (my $seq_obj = $seqio_obj->next_seq ) { my @agcts=split('',$seq_obj->seq); #将记录转换成字符数组 my $agct_length= length $seq_obj->seq; my $start_index=0; for(my $j=$start_index;$j<$agct_length-1;$j=$start_index){ #从序列数组最左边的元素开始搜索str串,直到序列数组结束 my @end_indexs=(0,0,0,0,0); #记录本次搜索过程中找到的中止下标,由于str的周期为2-6,共5个周期, #故需要一个长度为5的数组,来记录当周期分别2、3、4、5、6时的中止下标 my @repeated_times=(0,0,0,0,0); #记录当周期分别2、3、4、5、6时的重复周期数 for(my $i=2;$i<=6;$i=$i+1){ #分别搜索当周期数为2、3、4、5、6时的str 长度与周期数 my $breaks=0; for(my $j=$start_index;$j<$agct_length-1;$j=$j+$i){ my $k; for($k=0;$k<$i;$k++){ #比较一个周期内的全部元素,只有当一个周期的内全部元素相同时,才能算完整的重复了一个周期 if($j+$k+$i>=$agct_length-1 || ($agcts[$j+$k] ne $agcts[$j+$k+$i])){ $breaks=1; last; } } if($breaks!=1){ $repeated_times[$i-2]=$repeated_times[$i-2]+1; $end_indexs[$i-2]=$j+$k+$i; }else{ last; } } } my $max_end_index=0; my $max_repeat_peroid=0; #找出本次搜索到的最长STR串 for(my $i=2;$i<=6;$i=$i+1){ if($end_indexs[$i-2]>$max_end_index){ $max_end_index=$end_indexs[$i-2]; $max_repeat_peroid=$i; } } if($max_repeat_peroid!=0 and $start_index>8 and $max_end_index+9< $agct_length and $repeated_times[$max_repeat_peroid-2]+1 >5){ #本次搜索到了STR 串,将找到的STR串放入最终结果 #本次搜索到了STR 串,右移下标 STR长度 my $repeated_times = $repeated_times[$max_repeat_peroid-2]+1; my $start_bp = substr($seq_obj->seq,$start_index-8,8); my $repeated_element = substr($seq_obj->seq,$start_index,$max_repeat_peroid); my $end_bp = substr($seq_obj->seq,$max_end_index,8); my $str = "$start_bp\_$repeated_element\_$end_bp"; if( exists $STRS{$str}) { push @{$STRS{$str}}, $repeated_times[$max_repeat_peroid-2]+1; } else { my @nums=($repeated_times[$max_repeat_peroid-2]+1); $STRS{$str}=\@nums; } $start_index=$max_end_index+1; }else{ #本次搜索没有找到STR串,则右移下标 $start_index=$start_index+1; } } } open(OUT,">test3.fasta.ssr.result"); select OUT; print ("need_seq\tfind_cishu\tmeicichongfucishu\n"); foreach my $key (keys %STRS){ my $count = @{$STRS{$key}}; print ($key," ",$count ," "); foreach my $num (@{$STRS{$key}}){ print ($num,";"); } print ("\n"); } } found_str("test3.fasta");