求大神帮助

13171082187

最近在做组蛋白修饰的匹配情况，自己编了个程序但是运行较慢，所以相求大神帮忙
H3K18ac.bed
chr1        +   65092
chr1        +   65099
chr1    -   65099
chr1        +   65089
chr1        +        363660
chr1        +        857222
promoter.txt
OR4F5        NM_001005484        chr1         +        65091        73090
OR4F3        NM_001005224        chr1         +        363659        371658
SAMD11        NM_152486        chr1         +        857121        865120
NOC2L        NM_015658        chr1         -        875583        883582
KLHL17        NM_198317        chr1           +        891967        899966
PLEKHN1        NM_032129        chr1         +        897877        905876
令$cnt1=0想在H3K18ac.bed的第一列、第二列分别与promoter.txt第三列、第四列相同的情况下，如果(promoter.txt的第五列=<H3K18ac.bed的第三列<=promoter.txt的第六列时，$cnt会自增）
想要的结果如下
OR4F5        NM_001005484        chr1        +        65091        73090    2
OR4F3        NM_001005224        chr1        +        363659        371658   1
SAMD11        NM_152486        chr1        +        857121        865120   1
NOC2L        NM_015658        chr1        -        875583        883582   0
KLHL17        NM_198317        chr1        +        891967        899966   0
PLEKHN1        NM_032129        chr1        +        897877        905876   0
编写的程序如下：
#!/usr/bin/perl  -w
use strict;
use warnings;
open FILE, '<','H3K18ac.bed';      #启动子区域这个文件输入到FILE中                    
open AFILE, '<','promoter.txt';      #H3K18ac.bed文件输入到AFILE中     
open OUT,'>','readH3K18ac.bed';     #最后的结果输出到readH3K18ac.bed中
my @array1 = <FILE>;
my @array2 = <AFILE>;
foreach my $array1(@array1){
   chomp($array1);
   my @a1=split /\s+/,$array1;  #按空格键分割并存到数组@a1中
   my $cnt1=0;   #$cnt1=0
   foreach my $array2(@array2) {
    chomp($array2);
          my @a2 = split /\s+/,$array2;#按空格键分割并存到数组@a2中
            if ($a1[0] eq $a2[2] and $a1[1] eq $a2[3]){#在文件1第一列、第二列与文件2的第三列和第四列相同的情形下
             $cnt1++ if ($a1[2]>= $a2[4] and $a1[2] <= $a2[5]);#如果文件1的第三列大于等于文件2的第五列且小于等于文件2的第六列时，$cnt1会自增
                  }
         }print OUT "@a1 $cnt1\n"
}

ntwarren

把数据文件发上来,估计你的数据文件很大吧?
要优化也是要看了文件之后,才好做出有针对性的优化.

13171082187

回复 2# ntwarren


H3K18ac.bed的数据有点太大了，我上传了一个类似的文件 新建 好压 7Z 压缩文件.rar (2.56 MB, 下载次数: 7)

2015-01-21 16:46 上传

点击文件名下载附件

ntwarren

用这个吧,在我的机器上大概13秒完成计算,前提是基于观察,你的数据是按数字递增的.如果没有递增的话,就先要排序

#!/usr/bin/perl  -w
use strict;
use warnings;

open FILE, '<','H3K4me1.bed';      #启动子区域这个文件输入到FILE中                    
open AFILE, '<','promoter.txt';      #H3K18ac.bed文件输入到AFILE中     
open OUT,'>','readH3K18ac.bed';     #最后的结果输出到readH3K18ac.bed中
my $t = time();

my @array1 = <FILE>;
my @array2 = <AFILE>;
my $i=0;

for my $array2 (@array2){
    chomp $array2;
    my @a2 = split /\s+/,$array2;
    my $cnt1 = 0;

    for(;$i<@array1;$i++){
        chomp $array1[$i];
        my @a1 = split /\s+/, $array1[$i];
        if($a1[2]>=$a2[4]){
            if( $a1[2]<=$a2[5]){
                $cnt1++ if $a1[0] eq $a2[2] and $a1[1] eq $a2[3];
            } else{
                last;
            }
        }
    }
    print OUT "@a2 $cnt1\n";
}


my $e = time();
print ($e - $t);


回复 3# 13171082187


   

13171082187

回复 4# ntwarren


    非常感谢大神了