【问题】对N个相似文本进行计算分析？5000积分求解决~

huang6894

嗯，因为一些奇葩的需求，想请教一下各位大神~
某目录下有N个以下格式的文本：

1.txt：

1. 
   
2. Region 1 : chr1 21 33
   
3. 0 0 0 0 0 1 1 6 26 26 27 27 28
   
4. Region 2 : chr1 73 108
   
5. 2 2 2 2 2 3 3 3 3 6 7 7 7 7 7 7 7 7 7 7 19 19 19 19 19 19 19 19 19 19 20 20 20 20 20 21
   
6. Region 3 : chr1 109 124
   
7. 48 48 48 48 48 48 47 53 53 53 53 53 56 57 55 55
   
8. Region 4 : chr1 21 33
   
9. 32 44 44 44 44 46 47 47 47 47 47 54 61
   
10. Region 5 : chr1 83 99
    
11. 26 26 27 27 27 27 28 28 33 33 33 33 33 33 34 34 35
    

复制代码

2.txt：

1. 
   
2. Region 1 : chr1 21 33
   
3. 1 2 0 4 5 1 1 7 6 6 2 27 28
   
4. Region 2 : chr1 73 108
   
5. 12 2 2 22 2 23 3 13 13 6 7 7 71 7 6 7 7 7 7 4 19 19 19 19 19 19 19 19 19 19 20 20 20 20 20 21
   
6. Region 3 : chr1 109 124
   
7. 8 8 8 8 8 8 7 53 53 53 53 53 56 57 55 55
   
8. Region 4 : chr1 21 33
   
9. 2 44 4 4 4 6 7 7 7 7 7 4 1
   
10. Region 5 : chr1 83 99
    
11. 6 6 7 7 7 7 8 8 3 3 33 33 33 33 34 34 35
    

复制代码

他们的单行‘Region ****’都是一样的，但是双行（双行表示region行最后两个数字之间每个点的质量值）是不同的，我现在需要把这N个文件每个点的质量值取平均，每个点的质量值（P)都与这个平均值（Q）做比较，如果P值位于0~0.3Q之间，该点定义为flag1，如果P值位于0.3Q~0.6Q之间，定义改点为flag2，否则为flag3.（不会啊~）

最难的是：只有在一个区域（最少连续5个点)中有60%(3个)以上相同的flag时，认为该flag可信，输出这个区域和flag，如chr1 2-6  flag1
----------------------------------------------

编辑一下：举个例子：

-----------------------------
比如我们有4个文件都是题目的格式，然后：
在于Region 2 : chr1 73 108上，表示73到77这5bp的质量分别是：
1.txt：

1 1 1 1 1

2.txt：

2 2 2 2 2

3.txt：

3 3 3 3 3

4.txt：

4 4 4 4 4

-------------------------
那么每个点的平均质量值应该是：
2.5 2.5 2.5 2.5 2.5
--------------------------------------
那么对比条件：每个点的质量值（P)都与这个平均值（Q）做比较，如果P值位于0~0.3Q之间，该点定义为flag1，如果P值位于0.3Q~0.6Q之间，定义改点为flag2，否则为flag3.
----------------------------------------------------------------
应该是：
1.txt：

flag2 flag2 flag2 flag2 flag2


2.txt：

flag3 flag3 flag3 flag3 flag3

3.txt：

flag3 flag3 flag3 flag3 flag3

4.txt：

flag3 flag3 flag3 flag3 flag3
-------------------------------------
于是我们得到结果是：
1.txt    chr1  73-77  flag2
2.txt    chr1  73-77  flag3
3.txt    chr1  73-77  flag3
4.txt    chr1  73-77  flag3

---------------------------------------------------
上面举的例子是很特殊的，一般的例子可能是：

同样四个文本：
chr1 1-11
对比条件后是：
1.txt：

flag2 flag1 flag2 flag2 flag2  flag1  flag3  flag3  flag3   flag3  flag3  


2.txt：

flag3 flag3 flag3 flag2 flag2  flag2  flag2  flag3  flag3   flag3  flag3  #1-5:flag3===6-11:flag3(4-8:flag2 ?)===>优先判断为flag1，然后是flag2，最后才是flag3，也就是说我们需要首先认同4-8是flag2，然后再认为1-3是flag3，9-11是flag3

3.txt：

flag3 flag1 flag3 flag1 flag1  flag3  flag1  flag1  flag1   flag1  flag1  

4.txt：

flag3 flag3 flag3 flag3 flag3  flag3  flag3  flag3  flag3   flag3  flag3  
-------------------------------------
于是我们得到结果是：
1.txt    chr1  1-5  flag2
1.txt    chr1  6-11  flag3
2.txt    chr1  1-3  flag3
2.txt    chr1  4-8  flag2
2.txt    chr1  9-11  flag3
3.txt    chr1  1-11  flag1
3.txt    chr1  1-11  flag3
---------------------------------------------------
说实话，我一个人研究，估计一年也搞不出来。。。。只能跪求各位大神了。。。。。

pitonas

这个问题, 我不放弃 ~ {:2_172:}
基本完成
~ {:2_168:}

• 使用猫 (cat)
• 我不知道 flag 这个计算规则，所以我用我的规则
• 在第 40 行有一个错误
  

1. #!/usr/bin/perl -w
   
2. use 5.010;
   
3. 
   
4. my @txt = qw/1.txt 2.txt/;          # glob or ...
   
5. my @tmp = map { "$_.tmp" } @txt;    # tmp files
   
6. my ( $fh, @fh ) = map { open my $f, '<', $_; $f } @txt;
   
7. my @write_fh    = map { open my $f, '>', $_; $f } @tmp;
   
8. my ( $CHR, $BEGIN );
   
9. 
   
10. while (<$fh>) {
    
11.     ( $CHR, $BEGIN ) = (split)[ 3, 4 ];
    
12.     my @pos   = [ split /\s/, <$fh> ];
    
13.     my $index = $#{ $pos[0] };
    
14. 
    
15.     for my $f (@fh) {
    
16.         readline $f;
    
17.         push @pos, [ split /\s/, <$f> ];
    
18.     }
    
19. 
    
20.     my @ave = map {
    
21.         my ( $i, $sum ) = $_;
    
22.         $sum += $pos[$_][$i] for 0 .. $#pos;
    
23.         $sum / @pos;
    
24.     } 0 .. $#{ $pos[0] };
    
25. 
    
26.     for my $i ( 0 .. $#pos ) {
    
27.         my @flag = map {
    
28.             my $F = $pos[$i][$_] / ( $ave[$_] || 1 );
    
29.             [ $F = $F <= 0.3 ? 0 : $F > 0.6 ? 2 : 1, 1 ];
    
30.         } 0 .. $index;
    
31. 
    
32.         separate( $i, @flag );
    
33.     }
    
34. }
    
35. 
    
36. sub separate {
    
37.     my $index = shift;
    
38.     my @flag = ( [], [], [] );
    
39. 
    
40.     for my $f ( 1 .. 3 ) {
    
41.         for ( my $i = 0 ; $i <= $#_ - 4 ; $i++ ) {
    
42.             next if $_[$i][0] ne $f;
    
43.             $i += 7 and next if !$_[ $i + 2 ][1];
    
44.             my $end  = $i + 4;
    
45.             my $zero = 0;
    
46.             for my $x ( reverse $end - 1 .. $end ) {
    
47.                 $_[$x][1] ? last : $zero++;
    
48.             }
    
49. 
    
50.             my $count = grep { $_[$_][0] == $f } $i .. $end - $zero;
    
51.             next if $count < 3;
    
52.             @{ $_[$_] } = ( $f, 0 ) for $i .. $end - $zero;
    
53.             $i += 4 + $zero ? 5 - $zero : 0;
    
54.         }
    
55.     }
    
56. 
    
57.     for my $i ( 0 .. $#_ ) {
    
58.         !$flag[ $_[$i][0] ][0] || ( $i - $flag[ $_[$i][0] ][-1][-1] != 1 )
    
59.           ? push @{ $flag[ $_[$i][0] ] }, [$i]
    
60.           : ( $flag[ $_[$i][0] ][-1][1] = $i );
    
61.     }
    
62.     my $i;
    
63.     for my $f (@flag) {
    
64.         $i++;
    
65.         for my $range (@$f) {
    
66.             @$range <= 1 || $range->[1] - $range->[0] < 2 and next;
    
67.             my $range = join '-', map { $_ + $BEGIN } @$range;
    
68.             say { $write_fh[$index] } "$txt[$index]\t$CHR\t$range\tf$i";
    
69.         }
    
70.     }
    
71. }
    
72. 
    
73. close $_ for @write_fh;
    
74. `cat @tmp > result.txt`;    # cat
    
75. 
    
76. __DATA__
    
77.     ~ 0.3
    
78. 0.3 ~ 0.6
    
79. 0.6 ~
    

复制代码

q1208c

没看明白你说的具体意思, 能不能手工写一下你要的结果是什么样的?

另外, 你的双行数据这么多, 全是 质量 么? 有的还是一串0 ...

pitonas

my $师傅 jason680:

1. 问题本不难，你又不举例，
   
2. 描述又乱写，让人费疑猜？

复制代码

难的不是问题 ~ {:2_172:}
最难的是看懂问题 ~ {:2_168:}

xiumu2280

生信的表示看懂了

xiumu2280

他们的单行‘Region ****’都是一样的， #两行为一个单位，第一行相当于一个ID,chr1 21 33，表示在chr1上21到33的位置，下面第二行是对应在21到33位置上碱基的质量值也就是每个点的P值。
但是双行（双行表示region行最后两个数字之间每个点的质量值）是不同的，
我现在需要把这N个文件每个点的质量值取平均， #按照位置在多个文件内来个交叉求和再算平均值，这个平均值就是Q值。
每个点的质量值（P)都与这个平均值（Q）做比较， #即用P/Q
如果P值位于0~0.3Q之间，该点定义为flag1，如果P值位于0.3Q~0.6Q之间，定义改点为flag2，否则为flag3.（不会啊~）

最难的是：只有在连续5个点中有3个以上相同的flag时，认为该flag可信， #就是连续的5个位置中有3个以上在相同的flag区域，则输出该区域
输出这5个点的区域和flag，如chr1 2-6  flag1

我来帮你解释一下回复 3# pitonas


   

huang6894

@q1208c@pitonas两位大神，我很抱歉啊~
比如我们有4个文件都是题目的格式，然后：
在于Region 2 : chr1 73 108上，表示73到77这5bp的质量分别是：
1.txt：

1 1 1 1 1

2.txt：

2 2 2 2 2

3.txt：

3 3 3 3 3

4.txt：

4 4 4 4 4

-------------------------
那么每个点的平均质量值应该是：
2.5 2.5 2.5 2.5 2.5
--------------------------------------
那么对比条件：每个点的质量值（P)都与这个平均值（Q）做比较，如果P值位于0~0.3Q之间，该点定义为flag1，如果P值位于0.3Q~0.6Q之间，定义改点为flag2，否则为flag3.
----------------------------------------------------------------
应该是：
1.txt：

flag2 flag2 flag2 flag2 flag2


2.txt：

flag3 flag3 flag3 flag3 flag3

3.txt：

flag3 flag3 flag3 flag3 flag3

4.txt：

flag3 flag3 flag3 flag3 flag3
-------------------------------------
于是我们得到结果是：
1.txt    chr1  73-77  flag2
2.txt    chr1  73-77  flag3
3.txt    chr1  73-77  flag3
4.txt    chr1  73-77  flag3

huang6894

回复 5# xiumu2280


    xiumu2280大神我爱死你了~{:2_178:}

-----------------------
各位大神请原谅一个体育老师教语文的孩子~

xiumu2280

能看懂问题，我感到很开心 {:2_168:} 回复 7# huang6894


   

xiumu2280

想了一下，哈希切片简直就是为这个问题准备的·····
看看晚上帮你写下回复 7# huang6894