Perl (или R, или SQL): подсчитайте, как часто строка встречается в столбцах - PullRequest
9 голосов
/ 04 августа 2011

У меня есть текстовый файл, который выглядит так:

gene1   gene2   gene3
a       d       c
b       e       d
c       f       g
d       g       
        h
        i

(Каждый столбец является человеческим геном, и каждый содержит переменное количество белков (строки, обозначенные здесь буквами), которые могут связываться с этими генами).

То, что я хочу сделать, это подсчитать, в скольких столбцах представлена ​​каждая строка, вывести это число и все заголовки столбцов, например:

a   1   gene1
b   1   gene1
c   2   gene1 gene3
d   3   gene1 gene2 gene3
e   1   gene2
f   1   gene2
g   2   gene2 gene3
h   1   gene2
i   1   gene2

Я пытался выяснить, как это сделать в Perl и R, но пока безуспешно. Спасибо за любую помощь.

Ответы [ 5 ]

8 голосов
/ 04 августа 2011

Это решение кажется хакерским, но дает желаемый результат. Он основан на использовании пакетов plyr и reshape, хотя я уверен, что вы можете найти альтернативы base R. Хитрость в том, что функция melt позволяет нам сгладить данные в длинном формате, что позволяет легко (иш) манипулировать с этого момента.

library(reshape)
library(plyr)

#Recreate your data
dat <- data.frame(gene1 = c(letters[1:4], NA, NA),
                  gene2 = letters[4:9],
                  gene3 = c("c", "d", "g", NA, NA, NA)
                  )

#Melt the data. You'll need to update this if you have more columns
dat.m <- melt(dat, measure.vars = 1:3)

#Tabulate counts
counts <- as.data.frame(table(dat.m$value))

#I'm not sure what to call this column since it's a smooshing of column names
otherColumn <- ddply(dat.m, "value", function(x) paste(x$variable, collapse = " "))

#Merge the two together. You could fix the column names above, or just deal with it here
merge(counts, otherColumn, by.x = "Var1", by.y = "value")

Дает:

> merge(counts, otherColumn, by.x = "Var1", by.y = "value")
  Var1 Freq                V1
1    a    1             gene1
2    b    1             gene1
3    c    2       gene1 gene3
4    d    3 gene1 gene2 gene3
....
6 голосов
/ 04 августа 2011

В perl, при условии, что белки в каждом столбце не имеют дубликатов, которые необходимо удалить. (Если это так, вместо этого следует использовать хэш хэшей.)

use strict;
use warnings;

my $header = <>;
my %column_genes;
while ($header =~ /(\S+)/g) {
    $column_genes{$-[1]} = "$1";
}

my %proteins;
while (my $line = <>) {
    while ($line =~ /(\S+)/g) {
        if (exists $column_genes{$-[1]}) {
            push @{ $proteins{$1} }, $column_genes{$-[1]};
        }
        else {
            warn "line $. column $-[1] unexpected protein $1 ignored\n";
        }
    }
}

for my $protein (sort keys %proteins) {
    print join("\t",
        $protein,
        scalar @{ $proteins{$protein} },
        join(' ', sort @{ $proteins{$protein} } )
    ), "\n";
}

Читает из стандартного ввода, пишет в стандартный вывод.

5 голосов
/ 04 августа 2011

Один вкладыш (точнее 3 вкладыша)

ddply(na.omit(melt(dat, m = 1:3)), .(value), summarize, 
     len = length(variable), 
     var = paste(variable, collapse = " "))
1 голос
/ 04 августа 2011

В MySQL, например, так:

select protein, count(*), group_concat(gene order by gene separator ' ') from gene_protein group by protein;

при условии, что данные как:

create table gene_protein (gene varchar(255) not null, protein varchar(255) not null);
insert into gene_protein values ('gene1','a'),('gene1','b'),('gene1','c'),('gene1','d');
insert into gene_protein values ('gene2','d'),('gene2','e'),('gene2','f'),('gene2','g'),('gene2','h'),('gene2','i');
insert into gene_protein values ('gene3','c'),('gene3','d'),('gene3','g');
1 голос
/ 04 августа 2011

Если столбцов не много, вы можете сделать что-то подобное в sql.Вы в основном выровняете данные в 2-колоночную производную таблицу белок / ген, а затем суммируете их по мере необходимости.

;with cte as (
  select gene1 as protein, 'gene1' as gene
  union select gene2 as protein, 'gene2' as gene
  union select gene3 as protein, 'gene3' as gene
)

select protein, count(*) as cnt, group_concat(gene) as gene
from cte
group by protein
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...