Pulling a numbered chromosome range file given a gene location from a lookup table ideally from command line or R

Question

have a folder with roughly 1000 vcf files which have divided the human genome into chunks, the folder looks like this:

- main_programme_aggregated_chr1_1_194789237.vcf.gz
....
- main_programme_aggregated_chrX_14235_8759493845.vcf.gz
etc 

I will divide this file name into:

*_chr[1-22, X,Y]_A_B.* 

I have a separate tab divided file with the names of genes and their genomic coordinates (all the genes in the human genome on build 38) - file1:

chr1 19853 90835 WASH7P ENSG0000022454
...
chr10 38732 390853 JHBF ENSG0000578382
..
chrX 9532 908032 LOTY ENSG00005847

etc

I would like to use the gene name from the 4th column of file1 to find the relevant chunked file it falls within. So I would need to take column 1 of file 1, match it to chr[1-22, X or Y] and then find the chunk were the range of columns2+3 fall between (A-B) in the file name.

I do not need to go into the files themselves, just use their names.

I have been just manually looking up the gene coordinates and the relevant chunk but am sure this can be automated.

An example would be that I want to get the chunk file that the gene WASH7P falls under, so the relevant file would be main_programme_aggregated_chr1_1_194789237.vcf.gz.

Ideally I could do this from the command line and the relevant chunk file then be inputted as a variable into some downstream workflows I have.

Would anyone be able to help? Many thanks for your time

Answer 1

Here is a solution with bash:

#!/bin/bash

# 1 - create an array of coordinate files

# use nullglob in case there are no matching files,
# not that important for this toy case
shopt -s nullglob

# make a(n) array/list of vcf file in the current (./) directory
coordinate_files=(./*.vcf.gz)

# going through each line of "file1"
while read line; do
    
    # 2 - get the gene related info from "file1"
    # 2.1 - change line to an array (line split by whitespace)
    line_array=( $line )
    # 2.2 - get the 4th column, gene name and chr, the first column
    gene=${line_array[3]}
    chr=${line_array[0]}
    # 2.3 get coordinates
    start=${line_array[1]}
    end=${line_array[2]}
    
    # loop through each file in the array created at step 1
    for ((i=0; i<${#coordinate_files[@]}; i++)); do
        
        # parse the info hardcoded in the filename
        chr_in_filename=$(echo ${coordinate_files[$i]} | awk -F_ '{print $4}')
        start_in_filename=$(echo ${coordinate_files[$i]} | awk -F_ '{print $5}')
        end_in_filename=$(echo ${coordinate_files[$i]} | awk -F_ '{print $6}')
        end_in_filename=${end_in_filename%.vcf.gz}

        if [ "$chr" == "$chr_in_filename" ]
        then
            if [ "$start" -ge "$start_in_filename" ] && [ "$end" -le "$end_in_filename" ]
            then
                echo $gene
                echo ${coordinate_files[$i]}
            fi
        fi
    done
    
done <file1

The input is the two example coordinate files + "file1" in the question and the output is:

WASH7P
chr1 19853 90835 WASH7P ENSG0000022454
./main_programme_aggregated_chr1_1_194789237.vcf.gz

Answer 2

Below is a Python solution:

import os

files = {}

for root, dirs, files in os.walk('/path/to/vcf/files'):
    for file in files:
        # Example: main_programme_ver2_chr22_2874210_2947732_VEPannot.vcf.gz
        chrom = file.split('_')[3]
        start = int(file.split('_')[4])
        end = int(file.split('_')[5])
        files[os.path.abspath(file)] = {'chrom': chrom, 'start': start, 'end': end}
    break

genes = {}

with open(file1) as f:
    for line in f:
        fields = line.strip().split('\t')
        chrom = fields[0]
        start = int(fields[1])
        end = int(fields[2])
        gene = fields[3]
        genes[gene] = {'chrom': chrom, 'start': start, 'end': end}

mapping = {}

for gene in genes:
    for file in files:
        if (genes[gene]['chrom'] == files[file]['chrom'] and
            genes[gene]['start'] >= files[file]['start'] and
            genes[gene]['end'] <= files[file]['end']):
            mapping[gene] = file
            break

with open('mappings.txt', 'w') as f:
    for gene in mapping:
        f.write(gene + '\t' + file + '\n')