두 파일 사이의 교차점(파일 2의 값 범위에 속하는 파일 1의 값)

Question 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

Answer

원하는 작업을 수행하는 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

Question 2

긴 대기 시간을 피하기 위해 Linux에 자주 사전 설치되는 최소 SQL 엔진 SQLite를 사용하여 이 작업을 수행할 수 있습니다. 서버를 실행하지 않으며 파일에 저장된 SQL 데이터베이스와 작동합니다.

snp_data 및 창 디렉토리에서 다음을 수행하십시오.

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

일반적으로 대용량 파일의 경우 이 방법이 더 빠릅니다.

Answer