【问题标题】:Compute transitive closure计算传递闭包
【发布时间】:2013-01-10 20:52:33
【问题描述】:

我的成对 DNA 序列数据以下列方式显示相似性..

AATGCTA|1   AATCGTA|2
AATCGTA|2   AATGGTA|3
AATGGTA|3   AATGGTT|8
TTTGGTA|4   ATTGGTA|5
ATTGGTA|5   CCTGGTA|9
CCCGGTA|6   GCCGGTA|7
GGCGGTA|10  AATCGTA|2
GGCGGTA|10  TGCGGTA|11
CAGGCA|12   GAGGCA|13

上面是一个示例输入文件,原始文件是几百万行。 我希望输出根据行之间的公共元素对重叠的 id 进行聚类,并将它们输出到每个聚类的一行,如下所示

AATGCTA|1   AATCGTA|2   AATGGTA|3   AATGGTT|8   GGCGGTA|10  TGCGGTA|11
TTTGGTA|4   ATTGGTA|5   CCTGGTA|9
CCCGGTA|6   GCCGGTA|7
CAGGCA|12   GAGGCA|13 

我目前正在尝试使用mclsilix 对它们进行集群,但我没有成功运行 silix。但是 mcl 目前正在进行中,我想知道在 awk 或 perl 中是否有其他聪明的方法可以做到这一点。我很感激一些解决方案,谢谢。 (这是我的第一篇文章,如有错误请见谅)

只是为了让它更简单.. 很容易说我的意见是,

1   2
2   3
3   8
4   5
5   9
6   7
10  2
10  11
12  13

我希望输出是,

1   2   3   8   10  11
4   5   9
6   7
12  13

【问题讨论】:

  • 每个序列的共同点是什么?
  • 每个序列都显示出与其他序列的一些相似性,这就是它们的分组方式,最后我想知道哪个组/序列簇形成一个组,我有兴趣根据成对对它们进行分组火柴。需要主帮助才能将它们分组为一行,每行代表一个簇
  • 请试着解释一下你所说的相似性/分组/共同的意思,不清楚,我无法挑出结构,第一行开始AAT..,但第5个元素开始GGC..为什么?
  • 我已经修改了我的问题,看到最后,是否有意义
  • 很高兴不仅仅是我没看懂,我还是看不到你是如何定义行和列的:|

标签: awk bioinformatics dna-sequence transitive-closure


【解决方案1】:

我认为这不是真的,但无论如何:

use strict;
use warnings;
my @rows;
my %indx;
while(<DATA>) {
  chomp;
  my @v = split (/\s+/);
  my $r = {};
  for my $k (@v) {
    $r = $indx{$k}[0] if defined $indx{$k};
  }
  $r->{$v[0]}++;
  $r->{$v[1]}++;
  # print join(",", @v), "\n";
  push(@{$indx{$v[0]}}, $r);
  push(@{$indx{$v[1]}}, $r);
  push(@rows,  $r);
}
my %seen;
for my $r (@rows) {
  print (join("\t", keys %$r), "\n") if not $seen{$r}++;
}

__DATA__
AATGCTA|1   AATCGTA|2
AATCGTA|2   AATGGTA|3
AATGGTA|3   AATGGTT|8
TTTGGTA|4   ATTGGTA|5
ATTGGTA|5   CCTGGTA|9
CCCGGTA|6   GCCGGTA|7
GGCGGTA|10  AATCGTA|2
GGCGGTA|10  TGCGGTA|11
CAGGCA|12   GAGGCA|13

输出:

GGCGGTA|10  AATGCTA|1   AATGGTT|8   AATCGTA|2   AATGGTA|3   TGCGGTA|11
CCTGGTA|9   TTTGGTA|4   ATTGGTA|5
CCCGGTA|6   GCCGGTA|7
CAGGCA|12   GAGGCA|13

【讨论】:

  • 谢谢,似乎工作正常。我将更改您的脚本以读取输入文件并查看它是否有效,
  • 如果您将CAGGCA|12 GAGGCA|1 附加到您的输入,这将不起作用,即,每一行应该至少有一个唯一元素。但我会尽快更新。
  • 在我的输入中总是有一个独特的元素,所以上面可能就足够了
【解决方案2】:

如你所愿,awk 解决方案来了:

awk 'BEGIN{f=1}{c=0;
        for(i=1;i<=f;i++){
                if(!a[i]){
                        a[i]=$1" "$2; c=1; break;
                }else if(a[i]~$1){
                        a[i]=a[i]" "$2; c=1; break;
                }else if(a[i]~$2){ a[i]=a[i]" "$1; c=1; break; }
        }
        if(!c){ a[++f]=$1" "$2; c=0; }
} END{for(x=1;x<=f;x++)print a[x]}' DnaFile

上面的代码已经用您更简单的输入文件和原始文件(使用 CCGGTTAA 等)进行了测试,两者都有效。输出省略。

【讨论】:

    猜你喜欢
    • 2018-06-04
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2011-11-17
    • 1970-01-01
    • 2012-11-26
    • 2012-10-24
    • 1970-01-01
    相关资源
    最近更新 更多