Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-28T12:19:59.217Z
Answers.Text: 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?
Answers.UserDisplayName: tacrolimus
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-28T23:22:55.947Z
Answers.Text: @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.
Answers.UserDisplayName: Seung-been Steven Lee
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-29T08:36:51.580Z
Answers.Text: 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!
Answers.UserDisplayName: tacrolimus
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-24T20:54:05.577Z
Answers.LastActivityDate: 2021-06-28T23:19:41.007Z
Answers.LastEditDate: 2021-06-28T23:19:41.007Z
Answers.UserDisplayName: Seung-been Steven Lee
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-28T12:19:51.437Z
Answers.Text: 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?
Answers.UserDisplayName: tacrolimus
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-28T12:57:44.443Z
Answers.Text: 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?
Answers.UserDisplayName: haci
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-24T23:00:35.760Z
Answers.LastActivityDate: 2021-06-25T10:29:06.773Z
Answers.LastEditDate: 2021-06-25T10:29:06.773Z
Answers.UserDisplayName: haci
Body: <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>
  
Tags: <r><linux><shell><genomicranges>
Title: Pulling a numbered chromosome range file given a gene location from a lookup table ideally from command line or R
ThreadText: # 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 <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
  
  ```
  
  > **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**: @username_0, 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!  
  
  ---
  
  