☰

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_data的值和文件中和欄位的值決定的，例如，值為in的行對應於值為for的行，並且其行以 a 開頭的前綴。chromosomewindowchromosomesnp_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

重點是位置僅在每個染色體內是唯一的，因此我需要逐個染色體地比較這兩個文件（即分別針對每個染色體）。我怎樣才能做到這一點？

答案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 來完成此操作。它不運行伺服器，而是使用儲存在檔案中的 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

一般來說，對於大檔案來說這應該會更快。

相關內容