Вычислить переходное замыкание

У меня есть данные о парных последовательностях ДНК, которые показывают сходство следующим образом.

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

Выше приведен пример входного файла, исходный файл состоит из нескольких миллионов строк. Я хочу, чтобы вывод объединял перекрывающиеся идентификаторы на основе общих элементов между строками и выводил их в одну строку для каждого кластера, как показано ниже.

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

В настоящее время я пытаюсь сгруппировать их, используя mcl, а также silix, мне не удалось запустить 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

bioinformatics awk dna-sequence transitive-closure
person bala    schedule 10.01.2013    source источник
comment
Что делает каждую последовательность общей?   -  person Chris Seymour    schedule 11.01.2013
comment
Каждая последовательность показывает некоторое сходство с другой последовательностью, то есть то, как они сгруппированы, в конце концов, я хочу знать, какая группа / кластер последовательностей формирует одну группу, мне интересно сгруппировать их на основе попарных совпадений. Основная помощь требуется, чтобы сгруппировать их в одну строку, каждая строка представляет собой кластер.   -  person bala    schedule 11.01.2013
comment
Пожалуйста, попробуйте объяснить, что вы имеете в виду под сходством / сгруппированными / общими, это неясно, и я не могу выделить структуру, первая строка начинается AAT.., а 5-й элемент начинается GGC.. почему?   -  person Chris Seymour    schedule 11.01.2013
comment
Я изменил свой вопрос, посмотрим ближе к концу, есть ли в нем смысл   -  person bala    schedule 11.01.2013
comment
как определены 6 столбцов   -  person Kent    schedule 11.01.2013
comment
Рад, что не только я не понимаю, я до сих пор не понимаю, как вы определили строки и столбцы :|   -  person Chris Seymour    schedule 11.01.2013
comment
вход имеет только 2 столбца во входном файле, они определяют пары. Я хочу создать объединение пар на основе пересекающегося идентификатора, например. 1 2 - это пара 2 3 - следующая пара, общий номер, который соединяется между строками, - это номер 2, поэтому я могу присоединиться к ним на основе перекрывающегося номера 2, поэтому мой результат в этом, например. будет 1 2 3 (не читать как 3 столбца), я могу назвать их одним кластером. Это лучше ?   -  person bala    schedule 11.01.2013
comment
В вашем упрощенном примере я вижу, что эти 1 2 3 8 связаны, но, учитывая эту серию, я не вижу, насколько 10 соответствует этому правилу. Удачи!   -  person shellter    schedule 11.01.2013
comment
1 2 в первой строке и 10 2 в седьмой строке, 2 является общим между ними, что сгруппирует 10 в кластере и, таким образом, 11 из восьмой строки 10 11 (потому что 10 - это перекрывающееся число)   -  person bala    schedule 11.01.2013
comment
Это не кластеризация. Вы хотите вычислить транзитивное замыкание вашего отношения идентичности. Прямой подход предполагает использование некоторых записей сопоставления хэш-карты со списком идентичных записей, а затем печать всех уникальных списков в конце.   -  person Has QUIT--Anony-Mousse    schedule 11.01.2013

Ответы (2)


arrow_upward
1
arrow_downward

Думаю, это не совсем то, но все же:

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
person perreal    schedule 10.01.2013
comment
Спасибо, вроде работает. Я изменю ваш скрипт, чтобы прочитать входной файл и посмотреть, работает ли он, - person bala; 11.01.2013
comment
Это не сработает, если вы добавите CAGGCA|12 GAGGCA|1 в свой ввод, т.е. в каждой строке должен быть хотя бы один уникальный элемент. Но я скоро обновлю это. - person perreal; 11.01.2013
comment
В моем вводе всегда есть один уникальный элемент, поэтому приведенного выше может быть достаточно - person bala; 11.01.2013

arrow_upward
1
arrow_downward

как вы и хотели, вот и решение 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 и т. д.) работал. Вывод опущен.

person Kent    schedule 10.01.2013