没啥说的,直接上代码

#!/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");

 

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据