0

別のデータ操作に関する質問があります。したがって、タブ区切りデータのこの .gtf ファイルがあり、特定の機能を抽出する必要があります。以前は、各遺伝子の「エクソン」タイプごとに遺伝子 ID、POS1、および POS2 を抽出するだけでよかったので、これは簡単でした。同じことを行う必要がありますが、最初に遺伝子内の位置に関連する各エクソンの POS1 と POS2 を見つける必要があります。現在、列 POS1 と POS2 は、ゲノム全体の TYPE の位置に基づいて番号が付けられています (これが、番号が非常に高い理由です)。ストランドが - の場合、これは逆になります。PITG_00002 を見ると、停止コドンが開始コドンの前にあることがわかります。これは、すべてが + (テンプレート) ストランドに対して相対的に番号付けされているためです。データシートのサンプルは次のとおりです。

GENE ID     TYPE        POS1    POS2    STRAND
PITG_00003  start_codon 38775   38777   +   0
PITG_00003  stop_codon  39069   39071   +   0
PITG_00003  exon        38775   39071   +   .
PITG_00003  CDS         38775   39068   +   0
PITG_00004  start_codon 39526   39528   +   0
PITG_00004  stop_codon  41492   41494   +   0
PITG_00004  exon        39526   40416   +   .
PITG_00004  CDS         39526   40416   +   0
PITG_00004  exon        40486   40771   +   .
PITG_00004  CDS         40486   40771   +   0
PITG_00004  exon        40827   41494   +   .
PITG_00004  CDS         40827   41491   +   2
PITG_00002  start_codon 10520   10522   -   0
PITG_00002  stop_codon  10097   10099   -   0
PITG_00002  exon        10474   10522   -   .
PITG_00002  CDS         10474   10522   -   0
PITG_00002  exon        10171   10433   -   .
PITG_00002  CDS         10171   10433   -   2
PITG_00002  exon        10097   10114   -   .
PITG_00002  CDS         10100   10114   -   0

したがって、各遺伝子について、「開始コドン」TYPE の位置に対して 1 から番号をやり直す必要があります。残念ながら、STRAND にリストされている遺伝子の番号は逆になっています (たとえば、PITG_00002)。したがって、これらのケースでは、番号付けは start_codon の POS2 を基準にして 1 から開始し、exon の POS1 で終了する必要があります。

そのため、エクソンごとに新しい POS1 と POS2 を取得する必要があり、これを POSA と POSB と呼びます。

各エクソンの POSA を取得するには、次のようにします。

POS1 of "exon" - POS1 of "start_codon" + 1 = POSA

各エクソンのPOSBを取得するには、次のようにします。

POS2 of "exon" - POS1 of "start_codon" + 1 = POSB

例として PITG_00004 を使用します。

POSA = 39526-39526 + 1 = 1
POSB = 40416 - 39526 + 1 = 891

そして、各遺伝子の各エクソンに対して同じことを行い、その遺伝子の start_codon 位置を使用して番号付けをリセットします。マイナス ストランドの場合を除き、その場合は次のことを行う必要があります。

各エクソンの POSA を取得するには、次のようにします。

POS2 of "start_codon" - POS2 of "exon" + 1 = POSA

各エクソンのPOSBを取得するには、次のようにします。

POS1 of "start_codon" - POS1 of "exon" + 1 = POSB

最終的に私はこれを取得したいと思います：

PITG_00002 exon 1 49
PITG_00002 exon 90 352
PITG_00002 exon 409 426
PITG_00003 exon 1 297
PITG_00004 exon 1 891
PITG_00004 exon 961 1246
PITG_00004 exon 1302 1969

+ストランドに対してこれを1つの方法で、-ストランドに対して別の方法で行う方法がよくわかりません。私は最近、より頻繁に python を使用していますが、perl も使用できます。

4

2 に答える 2

1

パールソリューション。ハッシュを使用して、各遺伝子に関する情報を保存します。@idxs配列は、数式の繰り返しを避けるために使用されます。

#!/usr/bin/perl
use warnings;
use strict;
use feature qw(switch);

my %hash;
<>;                   # Skip header.
while (<>) {
    my ($id, $type, $pos1, $pos2, $strand, undef) = split;
    given ($type) {
        when ('start_codon') {
            $hash{$id}{start}  = [$pos1, $pos2];
            $hash{$id}{strand} = $strand;
        }
        when ('stop_codon') {
            $hash{$id}{stop}  = [$pos1, $pos2];
        }
        when ('exon') {
            push @{ $hash{$id}{exons} }, [$pos1, $pos2];
        }
    }
}

for my $id (sort keys %hash) {
    my @idxs = '+' eq $hash{$id}{strand} ? (0, 1) : (1, 0);
    for my $exon (@{ $hash{$id}{exons} }) {
        my $posa = 1 + abs $hash{$id}{start}[$idxs[0]] - $exon->[$idxs[0]];
        my $posb = 3 + abs $hash{$id}{start}[$idxs[1]] - $exon->[$idxs[1]];
        print "$id exon $posa $posb\n";
    }
}
于 2013-01-22T16:39:26.900 に答える
0

OK、これがあなたの問題の解決策です（少なくとも私が理解していることについては）：これは、Pythonでのデータ分析の現在のゴールデンスタンダードであるpandasライブラリ（ http://pandas.pydata.org/ ）に基づいています。

まず最初にデータをロードします:

data = pd.read_csv('genetest.csv', sep='\t',
                   converters={'STRAND': lambda s: s[0]})

変換されたものは、ストランド列から余分な文字を取り除き、+ または - のみを残しています。

次に、groupby 関数を使用して、配列を鎖の方向と遺伝子名で分離します。

groups = data.groupby(['STRAND', 'GENE_ID'])

これにより、同じストランド方向と遺伝子名を持つデータセットが返され、それぞれを個別に処理できます。したがって、それらを辞書の項目 (キーと値のペアのリスト) として反復処理し、操作します。

corrected = []
for (direction, gene_name), group in groups:
    print direction,gene_name
    # take the index of the element you are going to subtract to the others
    start_exon = group.index[group.TYPE=='start_codon'][0]
    # now you perform your normalization and put it back into your group
    if direction == '+':
        group['POSA'] = 1 + group.POS1 - group.POS1[start_exon]
        group['POSB'] = 1 + group.POS2 - group.POS1[start_exon]
    else:
        group['POSA'] = 1 - group.POS2 + group.POS2[start_exon]
        group['POSB'] = 1 - group.POS1 + group.POS2[start_exon]
    print group
    # put into the result array
    corrected.append(group)
# join them together to obtain the whole dataset with the POSA and POSB
new_data = pd.concat(corrected)
print new_data

そして、これはあなたが得るものです：

    GENE_ID     TYPE    POS1    POS2    STRAND  POSA    POSB
0   PITG_00003  start_codon     38775   38777   +   1   3
1   PITG_00003  stop_codon  39069   39071   +   295     297
2   PITG_00003  exon    38775   39071   +   1   297
3   PITG_00003  CDS     38775   39068   +   1   294
4   PITG_00004  start_codon     39526   39528   +   1   3
5   PITG_00004  stop_codon  41492   41494   +   1967    1969
6   PITG_00004  exon    39526   40416   +   1   891
7   PITG_00004  CDS     39526   40416   +   1   891
8   PITG_00004  exon    40486   40771   +   961     1246
9   PITG_00004  CDS     40486   40771   +   961     1246
10  PITG_00004  exon    40827   41494   +   1302    1969
11  PITG_00004  CDS     40827   41491   +   1302    1966
12  PITG_00002  start_codon     10520   10522   -   1   3
13  PITG_00002  stop_codon  10097   10099   -   424     426
14  PITG_00002  exon    10474   10522   -   1   49
15  PITG_00002  CDS     10474   10522   -   1   49
16  PITG_00002  exon    10171   10433   -   90  352
17  PITG_00002  CDS     10171   10433   -   90  352
18  PITG_00002  exon    10097   10114   -   409     426
19  PITG_00002  CDS     10100   10114   -   409     423

ところで、あなたは質問に間違った距離補正を書き留めましたが、そうすべきです

POS1 of "start_codon" - POS2 of "exon" + 1 = POSB

反転した文字列が幾何学的な意味を持つようにする (そして投稿した値を取得する)

于 2013-01-22T16:40:54.477 に答える