搜索
查看: 2634|回复: 0

perl实现二分法查询最近的坐标

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-4-10 12:43:36 | 显示全部楼层 |阅读模式
如果给你一个坐标,比如chr1的27378位点,你能写程序查询它离哪个基因的转录起始位点最近吗?
hg38的基因的转录起始坐标是:http://www.biotrainee.com/jmzeng/tmp/hg38.tss

我写的程序如下:

[Perl] 纯文本查看 复制代码
$my_pos=$ARGV[0];
sub nearest_search{
    my($array,$value,$low_index,$high_index)=@_;
    $high_index=@{$array}-1 unless $high_index;
	$low_index=0 unless $low_index;
	## first check whether the arrary is sorted or not !
	if ($array->[$low_index] > $array->[$high_index]){
        print "not sort array !!!\n" ;
        last;
    }
	## then check 
    if ($value < $array->[$low_index] ){
        if($value < ($array->[$low_index] + $array->[$low_index-1] )/2){
			return $low_index-1
		}else{
			return $low_index
		}
		
    }elsif($value > $array->[$high_index]){
		 if($value < ($array->[$high_index] + $array->[$high_index-1] )/2){
			return $high_index 
		}else{
			return $high_index+1
		}
	}

    $mid_index=int (($low_index+$high_index)/2);
	
    if ($array->[$mid_index] == $value) { 
        return $mid_index;
    }
    elsif ($array->[$mid_index] < $value) {
        &nearest_search($array,$value,$mid_index+1,$high_index);
    }
    else{
        &nearest_search($array,$value,$low_index,$mid_index-1);
    }
}
open FH,"hg38.tss";
while(<FH>){
	chomp;
	@F=split;
	$chr{$F[1]}->{$F[2]}=[($F[0],$F[4])]
}
close FH;
foreach(keys %chr){
	$sort_pos{$_}=[sort {$a<=>$b}keys %{$chr{$_}}];
}
#print join"\n",@{$sort_pos{'chr1'}}[0..10],'\n';
$nearest_index=&nearest_search($sort_pos{'chr1'},$my_pos);
$nearest_pos=${$sort_pos{'chr1'}}[$nearest_index];
$up_pos=${$sort_pos{'chr1'}}[$nearest_index-1];
$down_pos=${$sort_pos{'chr1'}}[$nearest_index+1];
print "$nearest_index\t$up_pos\t$nearest_pos\t$down_pos\n";


简单的测试了一下,都是很准确的:
[AppleScript] 纯文本查看 复制代码
[jianmingzeng@jade hg38]$ perl tmp.pl  27371
2	15436	27370	28366
[jianmingzeng@jade hg38]$ perl tmp.pl  27370
2	15436	27370	28366
[jianmingzeng@jade hg38]$ perl tmp.pl  27378
2	15436	27370	28366
[jianmingzeng@jade hg38]$ perl tmp.pl  273785
8	185958	204597	449678
[jianmingzeng@jade hg38]$ perl tmp.pl  2737852
117	2631042	2787665	3019482
[jianmingzeng@jade hg38]$ perl tmp.pl  27378526
625	27372852	27381333	27390644
[jianmingzeng@jade hg38]$ perl tmp.pl  273785268
3421	248904243		
[jianmingzeng@jade hg38]$ perl tmp.pl  2737852688
3421	248904243	


希望你们可以看明白,我不喜欢讲代码



上一篇:perl实现hash和数组互相嵌套-多维hash
下一篇:请有root服务器的朋友务必学习这个
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-22 22:28 , Processed in 0.117536 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.