1
$\begingroup$

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

$\endgroup$
2
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ Many thanks for taking the time to answer. I understand the principle behind this but the actual file chunks look like this: main_programme_ver2_chr22_2874210_2947732_VEPannot.vcf.gz (as an example), how would I adjust your script for this? $\endgroup$
    – tacrolimus
    Jun 28 at 12:19
  • 1
    $\begingroup$ The script looks for underscores in the filenames of the coordinate files, so $4 becomes the chromosome,chr_in_filename: so for main_programme_aggregated_chr1_1_194789237.vcf.gz, the chromosome is chr1 . You can play with the numbers in the for loop. But since the number of "fields" in the filenames in your original post and what you pasted in your comment is the same, the code should work fine. How do you run this code and what is the output? $\endgroup$
    – haci
    Jun 28 at 12:57
2
$\begingroup$

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')
$\endgroup$
3
  • $\begingroup$ Many thanks for taking the time to answer. I understand the principle behind this but the actual file chunks look like this: main_programme_ver2_chr22_2874210_2947732_VEPannot.vcf.gz (as an example), how would I adjust your script for this? $\endgroup$
    – tacrolimus
    Jun 28 at 12:19
  • 1
    $\begingroup$ @tacrolimus, if all the files look similar to your example (main_programme_ver2_chr22_2874210_2947732_VEPannot.vcf.gz), then above script should run just fine. Basically, the script will split each file's name by the underscore character ('_') and then picks the 3rd item as chromosome, 4th item as start, and 5th item as end (0-based indexing). Let me know in the comment if somehow it fails. $\endgroup$ Jun 28 at 23:22
  • $\begingroup$ Worked well thanks! I have accepted the shell one as that my more "native" language but really appreciate you taking the time to answer and has spurned me to look into python more! $\endgroup$
    – tacrolimus
    Jun 29 at 8:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.