intersección entre 2 files (valores en el file 1 que caen en el range de valores en el file 2)

Tengo un file llamado snp_data contiene datos del cromosoma SNP (polymorphism de un solo nucleótido) . Este es un file CSV delimitado en blanco y espacio de 3 columnas que tiene el siguiente formatting:

 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 

Tenga en count que para cada fila el campo snp_id tiene la forma Chr{chromosome}__{position} para los valores correspondientes de chromosome y position .

Tengo otro file llamado window contiene datos auxiliares. Este es un file CSV delimitado en espacios blancos de 5 columnas que tiene el siguiente formatting:

 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 

Observe la correspondencia entre la window y los files snp_data determinados por el valor del campo del chromosome en el file de window y los valores de los campos chromosome y snp_id en el file snp_data , por ejemplo, las filas con un valor de "Chr1" en la window corresponden a las filas snp_data con un valor de 1 para el chromosome y cuyas filas snp_id comienzan con un prefijo de Chr01__ .

Para cada fila en el file snp_data (cada snp dentro de cada cromosoma), si el valor de position esa fila cae dentro del range dado por los valores de start y end de cualquiera de las filas en la window para ese cromosoma en particular, seqname el seqname de seqname de la window file a la fila del file snp_data .

Para la input dada arriba, este sería mi resultado deseado:

 user@host:~$ cat desinetworking_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 

El punto principal es que las posiciones son únicas solo dentro de cada cromosoma, por lo que necesito comparar estos 2 files cromosómicos por cromosoma (es decir, para cada cromosoma por separado). ¿Cómo puedo hacer esto?

Aquí hay una secuencia de commands de Python que debe hacer lo que desee:

 #!/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 # ie `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 desinetworking 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)) 

Tenga en count que este script asume que todas sus filas son filas de datos. Si sus files de input se llaman snp_data y window entonces puede ejecutarlos así:

 python intersect_snp.py snp_data window 

Si sus files tienen filas de encabezado, entonces puede usar tail para omitir / eliminar los encabezados y ejecutarlo así:

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

Supongamos que este es su file 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 

Y que este es tu file de 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 

Entonces, si ejecutamos este command:

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

Obtenemos el siguiente resultado:

 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 

Para evitar grandes times de espera, puede hacer esto con el motor SQLite de SQL minimalist que con frecuencia está preinstalado en Linux. No ejecuta un server y funciona con bases de datos SQL que están almacenadas en files.

En tu snp_data & window directory haz:

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

Esto normaliza los espacios entre campos y los prepara para la import.

A continuación, importe estos datos en SQL y ejecute la consulta para get el resultado:

 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 aquí para detener la input]

Y finalmente:

 cat task.sql | sqlite3 my_database.db 

En general, esto debería ser más rápido para files grandes.