<p>Below is a Python solution:</p>
<pre><code>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'] &gt;= files[file]['start'] and
            genes[gene]['end'] &lt;= files[file]['end']):
            mapping[gene] = file
            break

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

CC BY-SA 4.0
2021-06-28T12:19:59.217Z
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?
tacrolimus
CC BY-SA 4.0
2021-06-28T23:22:55.947Z
@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.
Seung-been Steven Lee
CC BY-SA 4.0
2021-06-29T08:36:51.580Z
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!
tacrolimus
CC BY-SA 4.0
2021-06-24T20:54:05.577Z
2021-06-28T23:19:41.007Z
2021-06-28T23:19:41.007Z
Seung-been Steven Lee
<p>Here is a solution with <code>bash</code>:</p>
<pre><code>#!/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 &quot;file1&quot;
while read line; do
    
    # 2 - get the gene related info from &quot;file1&quot;
    # 2.1 - change line to an array (line split by whitespace)
    line_array=( <span class="math-container">$line )
    # 2.2 - get the 4th column, gene name and chr, the first column
    gene=$</span>{line_array[3]}
    chr=<span class="math-container">${line_array[0]}
    # 2.3 get coordinates
    start=$</span>{line_array[1]}
    end=${line_array[2]}
    
    # loop through each file in the array created at step 1
    for ((i=0; i&lt;${#coordinate_files[@]}; i++)); do
        
        # parse the info hardcoded in the filename
        chr_in_filename=<span class="math-container">$(echo $</span>{coordinate_files[<span class="math-container">$i]} | awk -F_ '{print $4}')
        start_in_filename=$</span>(echo <span class="math-container">${coordinate_files[$i]} | awk -F_ '{print $5}')
        end_in_filename=$</span>(echo <span class="math-container">${coordinate_files[$i]} | awk -F_ '{print $6}')
        end_in_filename=$</span>{end_in_filename%.vcf.gz}

        if [ &quot;<span class="math-container">$chr" == "$</span>chr_in_filename&quot; ]
        then
            if [ &quot;<span class="math-container">$start" -ge "$</span>start_in_filename&quot; ] &amp;&amp; [ &quot;<span class="math-container">$end" -le "$</span>end_in_filename&quot; ]
            then
                echo <span class="math-container">$gene
                echo $</span>{coordinate_files[$i]}
            fi
        fi
    done
    
done &lt;file1
</code></pre>
<p>The input is the two example coordinate files + &quot;file1&quot; in the question and the output is:</p>
<pre><code>WASH7P
chr1 19853 90835 WASH7P ENSG0000022454
./main_programme_aggregated_chr1_1_194789237.vcf.gz
</code></pre>

CC BY-SA 4.0
2021-06-28T12:19:51.437Z
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?
tacrolimus
CC BY-SA 4.0
2021-06-28T12:57:44.443Z
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?
haci
CC BY-SA 4.0
2021-06-24T23:00:35.760Z
2021-06-25T10:29:06.773Z
2021-06-25T10:29:06.773Z
haci
<p>have a folder with roughly 1000 vcf files which have divided the human genome into chunks, the folder looks like this:</p>
<pre><code>- main_programme_aggregated_chr1_1_194789237.vcf.gz
....
- main_programme_aggregated_chrX_14235_8759493845.vcf.gz
etc 
</code></pre>
<p>I will divide this file name into:</p>
<pre><code>*_chr[1-22, X,Y]_A_B.* 
</code></pre>
<p>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:</p>
<pre><code>chr1 19853 90835 WASH7P ENSG0000022454
...
chr10 38732 390853 JHBF ENSG0000578382
..
chrX 9532 908032 LOTY ENSG00005847
</code></pre>
<p>etc</p>
<p>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.</p>
<p>I do not need to go into the files themselves, just use their names.</p>
<p>I have been just manually looking up the gene coordinates and the relevant chunk but am sure this can be automated.</p>
<p>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.</p>
<p>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.</p>
<p>Would anyone be able to help? Many thanks for your time</p>

<r><linux><shell><genomicranges>
Pulling a numbered chromosome range file given a gene location from a lookup table ideally from command line or R
# Pulling a numbered chromosome range file given a gene location from a lookup table ideally from command line or R

(Asked by: username_0 on 2021-06-24)

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

**Tags:** r, linux, shell, genomicranges

---

## 2 answers

### Answer by username_2 on 2021-06-24. Score: 2 | Accepted answer

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 
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

> **Comments:**
> **username_0**: 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?  
> **username_2**: 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?  

---

### Answer by username_1 on 2021-06-24. Score: 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')

> **Comments:**
> **username_0**: 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?  
> **username_1**: @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.  
> **username_0**: 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!  

---

