snp_data
Contains 라는 파일이 있습니다.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}
chromosome
position
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가 붙습니다.chromosome
window
chromosome
snp_id
snp_data
"Chr1"
window
snp_data
1
chromosome
snp_id
Chr01__
snp_data
파일의 각 줄(각 염색체 내의 각 snp) 에 대해 해당 줄의 값 position
이 해당 특정 염색체에 대한 임의 줄의 값의 합으로 주어진 범위 내에 속하면 파일의 줄에 추가합니다.start
end
window
seqname
window
snp_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
일반적으로 대용량 파일의 경우 이 방법이 더 빠릅니다.