Chinaunix

标题: 字符串计算区间覆盖度问题 [打印本页]
作者: papagogogo    时间: 2014-08-18 12:10
标题: 字符串计算区间覆盖度问题
现在有一段DNA序列起始及终了位置信息，格式如下：
       （起始位置）（终了位置）
chr7        58507201        58510800
chr7        58662001        58674300
chr7        59197501        59208000
chr7        59698201        59729100
chr7        61172101        61188000
chr7        61500001        61630500
chr7        61836601        61842300
chr7        61869601        61884000
chr7        62198401        62207100
chr7        62308501        62339400
chr7        62492701        62540400
chr7        62596201        62645400
chr7        62691901        62708400

当然还有ch1，2，3......文件大小大概50M

将具体数值映射到宽度为100000的区间上，并计算每个区间的覆盖度。(......)(......)(......)(......)(......)(......)(......)(......)(......) .
我现在是这样做的，设置一个字符串00000000000000....000000，长度为最大的终了位置长度，然后用第一次substr，遍历整个数据，将字符串覆盖到的起始终了位置都设置成1，这样会得到诸如下列字符串：000001111100010111110000001111......111100000, 最后再按宽度框选起来，第二次用substr (0000)(0111)(1100)(0101)(1111)(0000)(0011)(11..)...(..11)(1100)(0000)依次输出每个区间的覆盖度即可。


但是现在更深的问题来了。针对上述的某个区间，假设落在区间上的情况如下：


                  
                   111111111   111111111111
              1111111111110011
          111
(000001111111111111100111111111111)
          22322333333333    332222222222   计数


括号是区间，上面的是具体数值映射到区间上，但是我希望设置一个计数器，计算每一位的数值，并且设置一个阈值，假设我希望大于2的保留，反之舍弃，那么将得到如下的结果：

                   111111111   111111111111
              1111111111110011
          111
(000001111111111111100111111111111)
          22322333333333    332222222222   计数

(000000010011111111100110000000000) 最终结果。

然后针对最终结果再进行覆盖度计算。可是问题来了，如果是个位数计数的话应该没问题好说，可是如果计数数值达到10位以上了，就会让整个字符串变长那自然是不合适的。而且我要处理的文件有五百多个，很可能重复的次数很高，请问如何实现呢。非常感谢。
作者: afukada    时间: 2014-08-18 13:36
可以不要用string

改用array來做嗎?
作者: q1208c    时间: 2014-08-18 13:42
楼主确定知道自己在说什么么?

你前面数据文件中的内容, 跟后面的二进制字符串有什么关系? 你最后又想得到一个什么结果?
作者: papagogogo    时间: 2014-08-18 13:49
回复 3# q1208c


    前面的数组，第二列是序列的起始位置，第三列是终了位置。一段序列，比如长度是10000，那我用一个长度为10000的0字符串去表示。我需要计算在这100000长度的序列中，每一个窗口（假设长度为100）中的覆盖度，这个覆盖度就是将其实位置和终了位置及中间的所有值映射在字符串上，我的做法是映射成1，然后最后框出每个区间，计算覆盖度。但是现在需求更复杂了，我就不会做了，因为需要加入一个阈值。
作者: pitonas    时间: 2014-08-18 13:49
听说有洞的才是高尔夫球场，没洞的就不算高尔夫球场。~ {:2_172:}
现在问题来了。     
问题已经难很了, 喔,

但理解问题比回答这个问题更困难 ... {:2_170:}
作者: q1208c    时间: 2014-08-18 13:55
回复 4# papagogogo

你的描述我本来明白了一点点, 但看了你的配的数字, 我不明白了.

我觉得你能不能慢点, 一步步来. 把问题分开, 一个个说.

因为我不是生物专业的, 无法理解你说的到底是什么.


从程序上看, 你的问题可以用一个 "与" 或者 "或" 操作来解决. 问题没有你想得那么复杂.

   
作者: afukada    时间: 2014-08-18 14:02
老實講

我相信你要達成的目的並不困難

但是你用的方法卻讓整件事變得相當困難

這個問題看了就知道是要算sequence的coverage

可能你有自己的想法來做

但是你現在的想法可能對你要的結果是花了很多不必要的功夫

還不一定可以得到

雖然建議大家在提問的時候要盡量把問題變成是所有人都可以懂得

但是以你現在這樣的提問法

不是bioinformatic的看不懂

bioinformatic的不知道你要幹嘛

這樣並不能解決問題

所以可以麻煩你直接說要計算的題目即便裡面包含了很多bioinformatic或是biology的術語

這樣才有機會幫到你
作者: papagogogo    时间: 2014-08-18 14:31
回复 3# q1208c


00000000000000000000000000000000000000000000000000000000000000000----->这是一段空白序列

---------<data1>--------
        起始    终了
chr1   4       8
chr1   11     12
chr1   17     19
chr1   22     31
chr1   39     46
.........

---------<data2>-----
chr1   6       7
chr1  15     23
chr1  29     34
chr1  38     49
chr1  51     53



         00000000000000000000000000000000000000000000000000000000000000000
data1 00011111001100001110011111111110000000111111110000000000000000000
data2 00000110000000111111111000001111110001111111111110111000000000000
      A 00011111001100111111111111111111110001111111111110111000000000000
      B 00011221001100112221122111112221110001222222221110111000000000000
      C 00000220000000002220022000002220000000222222220000000000000000000

将空白序列中标注出起始和终了的信息，A是我之前做的，只要有信息就记为1不考虑重复的问题。    B是现在需要的，需要计数，并且设置一个阈值假设上述例子中阈值为1，小于1的就不要了所以就得到序列C，我如何计数，并且如果这个数值大于10那么序列岂不是会变长？是不是就不能用substr去做了，谢谢。


while  (<INFILE>) {
        @a = split(" ");
      if($a[0] =~ m/chr1$/) {$chr1start[$count1] = $a[1];$chr1end[$count1]  =  $a[2];$count1 +=1;}
}

my $wide=100000;
$dna1="0"x$chr1end[-1];

foreach (0..$count1-1){
        $L1=$chr1end[$_]-$chr1start[$_]+1;
        substr($dna1,$chr1start[$_],$L1)="1"x$L1;
}
foreach $win(0..int($chr1end[-1]/$wide)){
  $D[$win]= substr($dna1,$wide*$win,$wide);
  $C[$win]=($D[$win]=~ tr/1g//);
  print OUTFILE1 $C[$win]/$wide*100;
  print OUTFILE1 "\n";

我原来是这样写的
作者: papagogogo    时间: 2014-08-18 14:33
回复 7# afukada

非常感谢你的回答，我也不是学生物出身，但是现在要做这类事情，所以描述起来可能比较被动。。很抱歉，看我楼上的解释不知道能否清楚？

   
作者: afukada    时间: 2014-08-18 15:17
如果你堅持要用string

還有a-z, A-Z這樣的可以來代表大於9的數字
作者: papagogogo    时间: 2014-08-18 15:54
已经改用hash做了，string貌似还是不科学

回复 11# afukada


   
作者: q1208c    时间: 2014-08-18 16:08
回复 8# papagogogo
希望可以解决你的问题.
  1. #!/usr/bin/perl
  2. #
  3. use warnings;
  4. use strict;
  5. use Number::Range;
  6. use Set::Scalar;

  8. my $file1 = "1.file";
  9. my $file2 = "2.file";
  10. my %hash1;
  11. my %hash2;

  13. open(FILE1, "<".$file1) or die("$file1, $!\n");
  14. while (my $line = <FILE1>) {
  15.     chomp $line;
  16.     my ($key, $start, $end) = split(/\s+/, $line);
  17.     if (defined $hash1{$key}) {
  18.         $hash1{$key} .= ",$start..$end";
  19.     } else {
  20.         $hash1{$key} = "$start..$end";
  21.     }
  22. }
  23. close(FILE1);

  25. open(FILE2, "<".$file2) or die("$file2, $!\n");
  26. while (my $line = <FILE2>) {
  27.     chomp $line;
  28.     my ($key, $start, $end) = split(/\s+/, $line);
  29.     if (defined $hash2{$key}) {
  30.         $hash2{$key} .= ",$start..$end";
  31.     } else {
  32.         $hash2{$key} = "$start..$end";
  33.     }
  34. }
  35. close(FILE2);

  37. foreach my $key (keys %hash1) {
  38.     my $range1 = Number::Range->new($hash1{$key});
  39.     my $range2 = Number::Range->new($hash2{$key});
  40.     my @number1 = $range1->range;
  41.     my @number2 = $range2->range;
  42.     my $set1 = Set::Scalar->new(@number1);
  43.     my $set2 = Set::Scalar->new(@number2);
  44.     my $setc = $set1->unique($set2);
  45.     my $setd = $set1->difference($setc);
  46.     print $set1, "\n";
  47.     print $set2, "\n";
  48.     print $setc, "\n";
  49.     print $setd, "\n";
  50. }
复制代码
作者: afukada    时间: 2014-08-18 16:17
回复 12# papagogogo

其實我會比較希望知道你真正整個流程到底想得到什麼東西

"比方"說你做上面這些事情是想知道有哪些區域有3個以上的coverage

那你現在這個方法就會浪費很多不必要的時間跟記憶體空間

會有其他的方法可以得到你想要的答案

所以我才會回覆你

把你想要得到的結果說出來

也許可以不用用這種作法
作者: chenhao392    时间: 2014-08-18 22:44
BEDtools, coverageBED



欢迎光临 Chinaunix (http://bbs.chinaunix.net/) Powered by Discuz! X3.2