☰

2 つのファイル間の共通部分 (ファイル 1 の値がファイル 2 の値の範囲内にある)

tag-icon • tag-icon text-processing python perl bioinformatics
2 つのファイル間の共通部分 (ファイル 1 の値がファイル 2 の値の範囲内にある)

という名前のファイルがありsnp_dataますSNP（一塩基多型）染色体データ。これは、次の形式の 3 列の空白で区切られた CSV ファイルです。

user@host:~$ cat snp_data

snp_id    chromosome  position
Chr01__912 1 912 1
Chr01__944 1 944 1
Chr01__1107 1 1107 1
Chr01__1118 1 1118 1
Chr01__1146 1 1146 1
Chr01__1160 1 1160 1
...
...
...
Chr17__214708367 17 214708367
Chr17__214708424 17 214708424
Chr17__214708451 17 214708451
Chr17__214708484 17 214708484
Chr17__214708508 17 214708508

各行のフィールドには、およびの対応する値のsnp_id形式があることに注意してください。Chr{chromosome}__{position}chromosomeposition

補助データを含む という名前の別のファイルがありますwindow。これは、次の形式の 5 列の空白で区切られた CSV ファイルです。

user@host:~$ cat window

seqname chromosome start end width
1 Chr1 1 15000 15000
2 Chr1 15001 30000 15000 
3 Chr1 30001 45000 15000
4 Chr1 45001 60000 15000 
5 Chr1 60001 75000 15000 
6 Chr1 75001 90000 15000 
...
...
...
199954 Chr17 214620001 214635000 15000
199955 Chr17 214635001 214650000 15000
199956 Chr17 214650001 214665000 15000
199957 Chr17 214665001 214680000 15000
199958 Chr17 214680001 214695000 15000
199959 Chr17 214695001 214708580 13580

ファイル内の フィールドの値と、ファイル内の フィールドと フィールドwindowの値によって決まる、 および ファイル間の対応に注意してください。たとえば、内の の値を持つ行は、内のの値を持つ の行に対応し、その行は のプレフィックスで始まります。snp_datachromosomewindowchromosomesnp_idsnp_data"Chr1"windowsnp_data1chromosomesnp_idChr01__

ファイル内の各行snp_data（各染色体内の各 SNP）について、その行の値が、その特定の染色体の内のいずれかの行のおよび値positionによって指定された範囲内にある場合、ファイルからのをファイルからの 行に追加します。startendwindowseqnamewindowsnp_data

上記の入力に対して、私が望む出力は次のようになります。

user@host:~$ cat desired_output

snp_id  chromosome  position   window
Chr01__912  1   912      1
Chr01__944  1   944      1
Chr01__1107 1   1107     1
...
...
...
Chr17__214708367 17 214708367   199959
Chr17__214708424 17 214708424   199959
Chr17__214708451 17 214708451   199959
Chr17__214708484 17 214708484   199959
Chr17__214708508 17 214708508   199959

重要な点は、位置は各染色体内でのみ一意であるため、これら 2 つのファイルを染色体ごとに (つまり、各染色体を個別に) 比較する必要があることです。これはどのように実行できますか?

答え1

必要なことを実行する Python スクリプトは次のとおりです。

#!/usr/bin/env python2
# -*- coding: ascii -*-
"""intersect_snp.py"""

import sys

# Read data from the SNP file into a list
snp_list = []
with open(sys.argv[1], 'r') as snp_file:
    for line in snp_file:
        snp_row = line.split() 
        snp_list.append(snp_row)

# Read data from the "window" file into a dictionary of lists
win_dict = {} 
with open(sys.argv[2], 'r') as win_file:
    for line in win_file:
        seqnames, chromosome, start, end, width = win_row = line.split()
        if chromosome not in win_dict:
            win_dict[chromosome] = []
        win_dict[chromosome].append(win_row)

# Compare data and compute results
results = []

# Iterate over snp data rows
for snp_row in snp_list:

    # Extract the field values for each snp row
    snp_id, chromosome, position = snp_row

    # Convert the chromosome to match the format in the "window" file
    # i.e. `1` -> `Chr1`
    chromosome_name = "Chr{}".format(chromosome)

    # Iterate over the "window" rows corresponding to that chromosome
    for win_row in win_dict.get(chromosome_name, []):

        # Extract the field values for each "window" row
        seqnames, chromosome, start, end, width = win_row

        # Perform the desired comparison
        if int(start) <= int(position) <= int(end):

            # If the comparison returns true, construct the result row
            result = [snp_id, chromosome, position, seqnames]
            results.append(result)
            break

# Print the output column headers
columns = ["snp_id", "chromosome", "position", "window"]
print(" ".join(columns))

# Print the results
for row in results:
    print(' '.join(row))

このスクリプトでは、すべての行がデータ行であると想定していることに注意してください。入力ファイルに名前が付けられている場合はsnp_data、window次のように実行できます。

python intersect_snp.py snp_data window

ファイルにヘッダー行がある場合は、tailヘッダーをスキップ/削除し、代わりに次のように実行することができます。

python intersect_snp.py <(tail -n+2 snp_data) <(tail -n+2 window)

snp_dataこれがあなたのファイルだと仮定します:

snp_id              chromosome  position
Chr01__912          1           912
Chr01__944          1           944
Chr01__1107         1           1107
...
...
...
Chr17__214708367    17          214708367
Chr17__214708424    17          214708424
Chr17__214708451    17          214708451
Chr17__214708484    17          214708484
Chr17__214708508    17          214708508

そして、これがあなたのwindowファイルです:

seqnames    chromosome  start       end         width
1           Chr1        1           15000       15000
2           Chr1        15001       30000       15000
3           Chr1        30001       45000       15000
4           Chr1        45001       60000       15000
5           Chr1        60001       75000       15000
...
...
...
199954      Chr17       214620001   214635000   15000
199955      Chr17       214635001   214650000   15000
199956      Chr17       214650001   214665000   15000
199957      Chr17       214665001   214680000   15000
199958      Chr17       214680001   214695000   15000
199959      Chr17       214695001   214708580   13580

次に、このコマンドを実行します。

python intersect_snp.py <(tail -n+2 snp_data) <(tail -n+2 window)

次のような出力が得られます。

snp_id chromosome position window
Chr01__912 Chr1 912 1
Chr01__944 Chr1 944 1
Chr01__1107 Chr1 1107 1
...
...
...
Chr17__214708367 Chr17 214708367 199959
Chr17__214708424 Chr17 214708424 199959
Chr17__214708451 Chr17 214708451 199959
Chr17__214708484 Chr17 214708484 199959
Chr17__214708508 Chr17 214708508 199959

答え2

長い待ち時間を回避するには、Linux にプリインストールされていることが多い最小限の SQL エンジン SQLite を使用します。SQLite はサーバーを実行せず、ファイルに保存されている SQL データベースで動作します。

snp_data および window ディレクトリで次の操作を実行します。

cat snp_data | tr -s " " > snp_data.csv
sed 's#Chr##g' window | tr -s " " > window.csv

これにより、フィールド間のスペースが正規化され、インポートの準備が整います。

次に、このデータを SQL にインポートし、クエリを実行して出力を取得します。

cat > task.sql
CREATE TABLE snp_data (snp_id text,chromosome int, position int);
CREATE TABLE window (seqname int,chromosome int, c_start int , c_end int, c_width int);

.mode csv
.separator "  "
.import snp_data.csv snp_data
.import window.csv window

.mode column
.header on
SELECT D.snp_id, D.chromosome, D.position, W.seqname FROM snp_data D, window W WHERE W.chromosome=D.chromosome AND D.position BETWEN W.c_start AND W.c_end;

[ここでCTRL+Dを押すと入力が停止します]

そして最後に：

cat task.sql | sqlite3 my_database.db

一般的に、大きなファイルの場合はこの方が高速になります。

関連情報