设为首页收藏本站
开启辅助访问
切换到窄版

 找回密码
 注册

QQ登录

只需一步,快速开始

搜索
查看: 4479|回复: 0

求一个计算测序reads对我特定的几个DNA区域覆盖度的程序

[复制链接]
发表于 2011-3-24 10:40:52 | 显示全部楼层 |阅读模式
10金钱
请教高手:我处理Solexa测序数据,我有一些感兴趣的dna片段区域,我要看这些区域分别被多少reads覆盖到。

reads文件
chr10   10000008        10000057        +
chr10   1000004 1000053 +
chr10   10000047        10000096        +
chr10   10000049        10000098        -
chr10   100000502       100000551       +
chr10   100000716       100000765       +
chr10   100000717       100000766       -

roi文件
chr10   100524501       100525000
chr10   100631501       100632000
chr10   100679001       100679500
chr10   100968501       100969000
chr10   100969001       100969500
chr10   101055501       101056000
我自己写的perl,没怎么学过数学,perl也刚开始学,有很多地方不太会写,还请高手指教,不胜感激。
#!/usr/bin/perl -w
use strict;
open READS, "reads_file" or die;
open ROIS, "roi_file" or die;
open OUT, ">coverage_file" or die;
my($read_flag,$roi_flag)=(1,1);
my($reads,$rois);
my($read_chr,$read_start,$read_end,$read_identfier);
my($roi_chr,$roi_start,$roi_end);
my $count = 0;
my $reads_num = 0;
my $overlap;
while(1){
if($roi_flag){
  $rois = <ROIS>; last unless $rois;
  chomp($rois);
  ($roi_chr,$roi_start,$roi_end)=split /\s+/, $rois;
  chomp($roi_chr,$roi_start,$roi_end);
  $roi_flag=0;
}
if($read_flag){
  $reads = <READS>; last unless $reads;
  chomp($reads);
  ($read_chr,$read_start,$read_end,$read_identfier)=split /\s+/, $reads;
  chomp($read_chr,$read_start,$read_end,$read_identfier);
  $read_flag=0;
}
if($roi_chr eq $read_chr){
  if($read_end<$roi_start){    #c01
   $read_flag=1;
  }elsif($read_start<$roi_start&&$read_end<$roi_end){ #c02
   $overlap = $read_end-$roi_start;
   if($overlap>30){
   $count++;
   }
   $read_flag=1;
  }elsif($read_start>=$roi_start&&$read_end<=$roi_end){ #0c3
   $count++;
   $read_flag=1;
  }elsif($read_start>$roi_start&&$read_end>$roi_end){ #c04
   $overlap = $read_end-$roi_start;
   if($overlap>30){
    $count++;
   }
    $read_flag=1;
  }elsif($read_start >= $roi_end){
   $reads_num=$count;
   print "$roi_chr\t$roi_start\t$roi_end\$reads_num\n";
   $count = 0;
   $roi_flag=1;
  }
}elsif($roi_chr gt $read_chr){
  $read_chr=1
}else{
  $roi_chr=1;
}
}
close ROIS;
close READS;
close OUT;

#                                    roi_start|*****************************|roi_end
#c01: read_start|*******|read_end
#c02:                        read_start|*******|read_end
#c03:                                             read_start|*******|read_end
#c04:                                                                          read_start|*******|read_end
#c05:                                                                                                 read_start|********|read_end

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

本版积分规则

QQ|申请友链|小黑屋|手机版|Archiver|生物信息学论坛 ( 蜀ICP备09031721号  

GMT+8, 2017-3-2 02:05 , Processed in 0.154308 second(s), 20 queries .

Powered by Discuz! X3

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表