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

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

snp_dataContains 라는 파일이 있습니다.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 및 창 디렉토리에서 다음을 수행하십시오.

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

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

관련 정보