1

我有一个看起来像这样的数据集

-9030   KIR3DX1
-75     SLC12A6
8005    C14orf79
-251    ARAP1
65994   EFNB1
-12111  SLC7A5
-11643  CAMK2G
-19749  PRPS2
-23324  MIR198
10012   LOC100506172
-77     CCDC88A
12171   MMP14

其中第 1 列表示第 2 列中的元素（基因）在任一方向上与 0 的距离（以碱基对为单位）。我想把这些数据放在一个 50 个碱基对的窗口中。

有什么建议么？谢谢

4

5 回答 5

2

程序（希望 Perl 没问题）：

#!/usr/bin/perl

# Create a histogram of some data

# Denis Howe 2012-07-03 - 2012-07-03 18:40

use strict;
use warnings;

# my $n_bins = 10;
my $width = 50;

# Read lines into @d
my @d = <DATA>; chomp @d;

# Split each line containing a digit into a pair
@d = map [split(/\s+/, $_)], grep /\d/, @d;

# Find range
my $min = 9E9; my $max = -9E9;
foreach (@d)
{
    $min = $_->[0] if ($_->[0] < $min);
    $max = $_->[0] if ($_->[0] > $max);
}

# Round down to multiple of $width
$min = int($min/$width) * $width;

# Ensure there's a bin for max value
# my $width = ($max*1.01 - $min) / $n_bins;
my $n_bins = int(($max - $min) / $width) + 1;

# Allocate data to bins
my @bin;
foreach (@d)
{
    push @{$bin[($_->[0]-$min)/$width]}, $_;
}

# Show content of each bin
foreach (0 .. $n_bins-1)
{
    next unless ($bin[$_]);             # Ignore empty bins
    printf "%6d - %6d", $min + $_*$width, $min + ($_+1)*$width;
    print map("  " . $_->[0] . ":" . $_->[1], @{$bin[$_]}), "\n";
}

__DATA__
-9030   KIR3DX1
-75     SLC12A6
8005    C14orf79
-251    ARAP1
65994   EFNB1
-12111  SLC7A5
-11643  CAMK2G
-19749  PRPS2
-23324  MIR198
10012   LOC100506172
-77     CCDC88A
12171   MMP14
EOF

输出：

-23300 - -23250  -23324:MIR198
-19750 - -19700  -19749:PRPS2
-12150 - -12100  -12111:SLC7A5
-11650 - -11600  -11643:CAMK2G
 -9050 -  -9000  -9030:KIR3DX1
  -300 -   -250  -251:ARAP1
  -100 -    -50  -75:SLC12A6  -77:CCDC88A
  8000 -   8050  8005:C14orf79
 10000 -  10050  10012:LOC100506172
 12150 -  12200  12171:MMP14
 65950 -  66000  65994:EFNB1

HtH

于 2012-07-03T17:47:22.330 回答
1

非常简单的单行：

perl -MPOSIX=floor -anE'push@{$f{floor($F[0]/50)}},$F[1]}{$,=" ";for(sort{$a<=>$b}keys%f){$i=$_*50;say"$i -",$i+49,": @{$f{$_}}"}'

请注意，此单行代码为测试数据生成正确的输出（例如，查看 -23324 MIR198，它肯定在 -23350 - -23301 中）：

-23350 - -23301 : MIR198
-19750 - -19701 : PRPS2
-12150 - -12101 : SLC7A5
-11650 - -11601 : CAMK2G
-9050 - -9001 : KIR3DX1
-300 - -251 : ARAP1
-100 - -51 : SLC12A6 CCDC88A
8000 - 8049 : C14orf79
10000 - 10049 : LOC100506172
12150 - 12199 : MMP14
65950 - 65999 : EFNB1
于 2012-07-03T20:33:43.917 回答
0 
s='''
-9030   KIR3DX1
-75     SLC12A6
8005    C14orf79
-251    ARAP1
65994   EFNB1
-12111  SLC7A5
-11643  CAMK2G
-19749  PRPS2
-23324  MIR198
10012   LOC100506172
-77     CCDC88A
12171   MMP14
'''

import re
from collections import defaultdict

bin = defaultdict( list )
for distance, gene in re.findall('^(\S+)\s+(\S+)',s,re.M):
    bin[int(distance)//50].append(gene)

print( bin )
于 2012-07-03T16:03:10.187 回答
0

这行得通吗？

 from collections import defaultdict

 binner=defaultdict(list)
 with open(datafilename) as f:
     for line in f:
         i=int(line.split()[0])
         binner[i//50].append(line)

我只是将整条线装箱，因为我不知道您实际上想要保留哪些信息（对于混乱的数据集，这有点不清楚）......

于 2012-07-03T16:10:18.093 回答
0

我将假设您希望基因名称按碱基对中的距离进行分类：

from collections import defaultdict, Counter

bins = defaultdict(Counter)
binsize = 50

with open(datafile) as inf:
    for line in inf:
        data = line.split('<', 1)[0]
        offset, name = data.split()
        bins[int(offset)//binsize][name] += 1

然后

keys = sorted(bins)
for key in keys:
    values = ', '.join('{1} {0}'.format(a,b) for a,b in bins[key].most_common())
    print('{:>7} - {:>7} : {}'.format(binsize*key, binsize*(key+1)-1, values))

在您的样本数据上导致

 -23350 -  -23301 : 1 MIR198
 -19750 -  -19701 : 1 PRPS2
 -12150 -  -12101 : 1 SLC7A5
 -11650 -  -11601 : 1 CAMK2G
  -9050 -   -9001 : 1 KIR3DX1
   -300 -    -251 : 1 ARAP1
   -100 -     -51 : 1 CCDC88A, 1 SLC12A6
   8000 -    8049 : 1 C14orf79
  10000 -   10049 : 1 LOC100506172
  12150 -   12199 : 1 MMP14
  65950 -   65999 : 1 EFNB1
于 2012-07-03T16:15:15.497 回答