|
#2
fall_bernana2020-05-26 17:17
回复 楼主 白鸦
你这个相当于每<RANGE> 就重新读取一边INPUT1 .INPUT2
应该先存放<RANGE>  程序代码：my %range; while(<RANGE>){ chomp; my @s=split; $range{$s[0]}{$s[1]}{$s[2]}=1;#或者把区间长度也当一个键值存入,超出区间长度的肯定就是不要的.反正就是把循环数降低就怎么来 } while(<INPUT1>){ chomp; next if(/#/); next if(/^(?!.*PASS)/); my @t=split; if (exists $range{$t[0]}){#此处能过滤掉大部分数据,相当于给$s[0]建了一个索引 foreach (keys %{$range{$t[0]}}){ if ......这里做你定义的操作 } } } seek INPUT1,0,0; seek RANGE,0,0; 如果是个反复查询的功能.可以存到数据库里.建个索引来查询更简单. 类似于 select count(*) from ** where *=$s[0] and *>$s[1] and *<$s[2]; [此贴子已经被作者于2020-5-26 17:48编辑过] |
效率实在太低了，请各位前辈给出一点宝贵的意见优化下程序，或是更好的方案 感激不尽!
以下是统计的部分,因为其它部分太长就不贴出了来了；
RANGE：捕获区间(第一列为染色号，如chr1 ，第二、三列为左闭右开的区间如2、5即为 [ 2 , 5 ) ； )
INPUT1、INPUT2:进行筛选与比对的vcf文件；
while(<RANGE>){
chomp;
my @s=split;
while(<INPUT1>){
chomp;
next if(/#/);
next if(/^(?!.*PASS)/);
my @t=split;
next if($s[0]ne$t[0]or$s[1]gt$t[1]or$t[1]ge$s[2]);
my $id=join "-",$t[0],$t[1],$t[3],$t[4];
if(/Dels=0.00/){
$hash_s{$id}=0;
$A_snp++;}
elsif(/^(?!.*Dels)/){
$hash_i{$id}=0;
$A_indel++;}
}
seek INPUT1,0,0;
}
seek RANGE,0,0;
while(<RANGE>){
chomp;
my @s=split;
while(<INPUT2>){
chomp;
next if(/#/);
next if(/^(?!.*PASS)/);
my @t=split;
next if($s[0]ne$t[0]or$s[1]gt$t[1]or$t[1]ge$s[2]);
my $id=join "-",$t[0],$t[1],$t[3],$t[4];
if(/Dels=0.00/){
$B_snp++;}
if(/^(?!.*Dels)/){
$B_indel++;}
if(exists $hash_s{$id}){
$S_share++;
}
elsif(exists $hash_i{$id}){
$I_share++;
}
}
seek INPUT2,0,0;
}
[此贴子已经被作者于2020-4-29 09:50编辑过]