标题:[新人求教]统计两个VCF文件在捕获区域中的一致性与SNP、INDEL数量
只看楼主
白鸦
Rank: 1
等 级:新手上路
帖 子:4
专家分:0
注 册:2020-4-28
 问题点数:0 回复次数:3 
[新人求教]统计两个VCF文件在捕获区域中的一致性与SNP、INDEL数量
为了能够得到INPUT1、2文件在捕获区间RANGE中的一致性,并统计其中的质量合格的SNP与INDEL数,于是写出了以下的一段程序...
效率实在太低了,请各位前辈给出一点宝贵的意见优化下程序,或是更好的方案  感激不尽!

以下是统计的部分,因为其它部分太长就不贴出了来了;
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编辑过]

搜索更多相关主题的帖子: 统计 if 文件 捕获 next 
2020-04-29 09:33
fall_bernana
Rank: 11Rank: 11Rank: 11Rank: 11
等 级:贵宾
威 望:17
帖 子:240
专家分:2086
注 册:2019-8-16
得分:0 
回复 楼主 白鸦
你这个相当于每<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编辑过]

2020-05-26 17:17
白鸦
Rank: 1
等 级:新手上路
帖 子:4
专家分:0
注 册:2020-4-28
得分:0 
回复 2楼 fall_bernana
恩,确实是这样的,这个程序跑起来效率非常低下
我刚学习所以还不太了解,请问应该怎么办才好呢
2020-05-26 17:23
白鸦
Rank: 1
等 级:新手上路
帖 子:4
专家分:0
注 册:2020-4-28
得分:0 
回复 2楼 fall_bernana
好的,非常感谢!
我马上测试一下
2020-05-26 17:52



参与讨论请移步原网站贴子:https://bbs.bccn.net/thread-501323-1-1.html




关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.047422 second(s), 7 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved