Answers.Body: <p>For these requirements, it might be possible to simply use <a href="https://asia.ensembl.org/info/docs/tools/vep/script/vep_filter.html" rel="nofollow noreferrer">filter_vep</a> included in the VEP package for this. Note that multiple <code>--filter</code> flags are ANDed together (i.e. all filters must pass for a line to be printed). For example:</p>
  <pre><code>for i in *.vep.txt.gz; do
      filter_vep \
          -i &quot;${i}&quot; \
          -o &quot;${i%.vep.txt.gz}.filtered.vep.txt&quot; \
          --filter &quot;IMPACT is LOW&quot; \
          --filter &quot;CLIN_SIG matches likely_pathogenic&quot;
  done
  </code></pre>
  <p>If you then wanted to merge the output and include the filename, you could just use <a href="/questions/tagged/awk" class="post-tag" title="show questions tagged &#39;awk&#39;" aria-label="show questions tagged &#39;awk&#39;" rel="tag" aria-labelledby="awk-container">awk</a>. With the following in a file called <code>merge_filtered_vep_txt.awk</code>:</p>
  <pre><code>BEGIN {
      FS=OFS=&quot;\t&quot;
      f=1
  }
  
  /^##/ {
      next
  }
  
  sub(/^#/, &quot;&quot;) &amp;&amp; f {
      print &quot;ID&quot;, $0
      f=0
      next
  }
  
  {
      print FILENAME, $0
  }
  </code></pre>
  <p>Run using:</p>
  <pre><code>awk -f merge_filtered_vep_txt.awk /path/to/*.filtered.vep.txt
  </code></pre>
  <p>Python can also open gzip-compressed files for reading (and writing), using the <a href="https://docs.python.org/3/library/gzip.html" rel="nofollow noreferrer">gzip</a> module in the standard library. These can be parsed line-by-line to avoid memory errors. Here's one way using the <a href="https://docs.python.org/3/library/csv.html" rel="nofollow noreferrer">csv</a> module (also in the standard library):</p>
  <pre><code>import argparse
  import contextlib
  import csv
  import gzip
  import pathlib
  import sys
  </code></pre>
  <pre><code>def get_argument_parser():
  
      parser = argparse.ArgumentParser(add_help=False)
  
      group = parser.add_argument_group('input options')
      group.add_argument(
          &quot;-i&quot;,
          &quot;--input&quot;,
          dest='glob_pattern',
          type=str,
          default=&quot;**/*.vep.txt.gz&quot;,
          metavar=&quot;STR&quot;,
          help=&quot;The relative glob pattern (default: '%(default)s')&quot;,
      )
      group.add_argument(
          &quot;-d&quot;,
          &quot;--dir&quot;,
          dest='glob_path',
          type=pathlib.Path,
          metavar=&quot;DIR&quot;,
          default='.',
          help=&quot;Glob files in DIR (default: '%(default)s')&quot;,
      )
  
      group = parser.add_argument_group('output options')
      group.add_argument(
          &quot;-o&quot;,
          &quot;--out&quot;,
          type=pathlib.Path,
          metavar=&quot;FILE&quot;,
          default='-',
          help=&quot;Write the filtered output to FILE (default: stdout)&quot;,
      )
      group.add_argument(
          &quot;-f&quot;,
          &quot;--force&quot;,
          action='store_true',
          help=&quot;Overwrite the output file if it exists&quot;,
      )
  
      group = parser.add_argument_group('additional options')
      group.add_argument(
          &quot;-h&quot;,
          &quot;--help&quot;,
          action=&quot;help&quot;,
          help=&quot;Show this help message and exit&quot;,
      )
  
      return parser
  </code></pre>
  <pre><code>def main():
  
      parser = get_argument_parser()
      args = parser.parse_args()
  
      with contextlib.ExitStack() as stack:
          if args.out.name != '-':
              out = stack.enter_context(open(args.out, 'wt' if args.force else 'xt'))
          else:
              out = sys.stdout
  
          # Get a list of files from the directory specified by args.glob_path
  
          file_list = sorted(args.glob_path.glob(args.glob_pattern))
  
          # Build up a unique list of header fields from each input file
  
          header_cols = ['ID']
          for gz_file in file_list:
              with gzip.open(gz_file, 'rt') as infile:
                  for line in infile:
                      if not line.startswith('##'):
                          break
  
                  assert line.startswith('#')
  
                  for field in line[1:].rstrip().split('\t'):
                      if field not in header_cols:
                          header_cols.append(field)
  
          # Create a csv.DictWriter() and write the header fields
  
          writer = csv.DictWriter(out, fieldnames=header_cols, delimiter='\t')
          writer.writeheader()
  
          # Now re-open each file in the arguments list for reading
  
          for gz_file in file_list:
              with gzip.open(gz_file, 'rt') as infile:
  
                  # Skip the info header lines
                  line = next(infile)
                  while line.startswith('##'):
                      line = next(infile)
  
                  # Get the column headers
                  fieldnames = line.rstrip().removeprefix('#').split('\t')
  
                  # Parse each input file and store the desired rows in memory
                  reader = csv.DictReader(infile, fieldnames=fieldnames, delimiter='\t')
                  rows = []
                  for row in reader:
  
                      impact = row.get('IMPACT')
                      clinsig = row.get('CLINSIG')
  
                      if impact == 'LOW' and clinsig == 'likely pathogenic':
                          rows.append({'ID': gz_file, **row})
  
                  # Write the list of rows to the output file
                  writer.writerows(rows)
  </code></pre>
  <pre><code>if __name__ == '__main__':
      main()
  </code></pre>
  <p>With the above in a file called <code>filter_vep_txt.py</code>:</p>
  <pre><code>usage: filter_vep_txt.py [-i STR] [-d DIR] [-o FILE] [-f] [-h]
  
  input options:
    -i STR, --input STR  The relative glob pattern (default: '**/*.vep.txt.gz')
    -d DIR, --dir DIR    Glob files in DIR (default: '.')
  
  output options:
    -o FILE, --out FILE  Write the filtered output to FILE (default: stdout)
    -f, --force          Overwrite the output file if it exists
  
  additional options:
    -h, --help           Show this help message and exit
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-16T07:06:43.247Z
Answers.Text: thank you in python code, what is the line code for determining the list of files of vep.txt.gz? @Steve
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-16T10:51:17.267Z
Answers.Text: vep.txt.gz? ----> I suspect is just one file named vep.txt see: https://stackoverflow.com/questions/31028815/how-to-unzip-gz-file-using-python ; https://unix.stackexchange.com/questions/93139/can-i-zip-an-entire-folder-using-gzip
Answers.UserDisplayName: pippo1980
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-16T11:45:20.700Z
Answers.Text: @Maj Maybe I haven't understood correctly, but I think that's really a question for you. I.e. what files would you like to process? How would you like to come up with the list? Would you prefer some sort of glob pattern or would you prefer to supply a file of filenames, one per line? Perhaps you could [edit](https://bioinformatics.stackexchange.com/posts/20370/edit) your question to include some details?
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-17T21:21:25.127Z
Answers.Text: I've edited the question hope it's now more clear and also hope you can help me @Steve
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T04:28:59.963Z
Answers.Text: @Maj Please see the updated code above. It's not clear if each file contains the same header fields, so the above reads the list of input files twice to ensure we get a unique list of header columns. Run using `python filter_vep_txt.py -o out.txt *.vep.txt.gz` or similar.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T06:55:40.070Z
Answers.Text: thanks a lot @Steve, sorry but I can not know how your code finds the files. is it for one files or for all?
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T07:07:52.800Z
Answers.Text: What is the command to locate the files? @Steve
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T08:43:39.303Z
Answers.Text: @Maj You'll need to supply one or more files as positional arguments on the command line. If all of your files are in the current directory, you should be able to just use the above glob pattern. If these are under multiple directories, you could use: `python filter_vep_txt.py -o out.txt /path/to/files/*.vep.txt.gz /path/to/more/files/*.vep.txt.gz` for example.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T10:56:32.380Z
Answers.Text: ok thank you @Steve now it is clear, I saved the first 4 blocks of python code in filter_vep_txt.py but the last block code you wrote how can I run it? sorry I'm new in this field
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T11:44:29.533Z
Answers.Text: @Maj the last block is just the usage. You can print it using: `python filter_vep_txt.py -h`
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T13:32:47.880Z
Answers.Text: @Maj In my code above, you may need to change `'Likely pathogenic'` to `'likely pathogenic'`. Is your CLINSIG column actually a comma-separated list? I'm just curious, but what options did you use to run VEP? Maybe you could include this in your question?
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-22T07:08:56.537Z
Answers.Text: hi @Steve, when I run the code by python filter_vep_txt.py -o out.txt *.vep.txt.gz, It shows me this error>> OSError: [Errno 22] Invalid argument: '*.vep.txt.gz'.     NOTE (all the files are in the same directory)
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-22T12:32:59.773Z
Answers.Text: @Maj You just have a typo in your command, use: `python filter_vep_txt.py -o out.txt *.vep.txt.gz`. Note the star/asterisk in the glob pattern.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-22T19:22:47.387Z
Answers.Text: @Steve yeah I tried with one file and with all files but still I get the same error ((OSError: [Errno 22] Invalid argument: '*.vep.txt.gz'))
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-22T19:25:48.700Z
Answers.Text: Traceback (most recent call last):
    File "C:\Users \Desktop\files\filter_vep_txt.py", line 123, in <module>
      main()
    File "C:\Users\Desktop\files\filter_vep_txt.py", line 75, in main
      with gzip.open(vep_file, 'rt') as infile:
    File "C:\Program Files\WindowsApps\PythonSoftwareFoundation.Python.3.10_3.10.2544.0_x64__qbz5n2kfra8p0\lib\gzip.py", line 58, in open     binary_file = GzipFile(filename, gz_mode, compresslevel)
    File "C:\Program Files\WindowsApps\PythonSoftwareFoundation.Python.3.10_3.10.2544.0_x64__qbz5n2kfra8p0\lib\gzip.py", line 174, in __init__
  @Steve
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-22T19:26:46.187Z
Answers.Text: fileobj = self.myfileobj = builtins.open(filename, mode or 'rb')
  OSError: [Errno 22] Invalid argument: '*.vep.txt.gz'
  @Steve sorry for many questions
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-23T12:57:43.697Z
Answers.Text: @Maj I think the problem is just that your shell can't expand the glob pattern. I think the simpler solution is to just have the script find the files. Please see the changes above. Use like: `python filter_vep_txt.py -d /path/to/files -o out.txt`, and it will try to glob the files recursively.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-25T20:19:16.377Z
Answers.Text: @Steve the code is working now, but when I open the out.txt file, I found only the word (ID) and I did not find any other information, I ran this command  (python filter_vep_txt.py -d *.vep.txt.gz -o out.txt) in the same files path
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-26T02:56:30.220Z
Answers.Text: @Maj Use `-i` instead of `-d`: `python filter_vep_txt.py -i '*.vep.txt.gz' -o out.txt`. Be sure to use single quotes to prevent it from being expanded by your shell. Failing that, add a `print(line)` before the `assert` statement and post the output. I'm not expecting VEP output files to contain blank lines in the header, but maybe your output is slightly different.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-26T11:35:32.143Z
Answers.Text: @Steve thanks a lot for your reply... I get this error and when I open out.txt file I found only the header of files (columns names) ..... Traceback (most recent call last):
    File "C:\Users\t\Desktop\files\filter_vep_txt.py", line 122, in <module>
      main()
    File "C:\Users\t\Desktop\files\filter_vep_txt.py", line 116, in main
      rows.append({'ID': vep_file, **row})
  NameError: name 'vep_file' is not defined. Did you mean: 'gz_file'?
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-26T12:47:49.840Z
Answers.Text: @Maj Sorry, typo. Yes I did indeed mean 'gz_file'.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-16T01:06:45.570Z
Answers.LastActivityDate: 2023-01-26T12:46:18.530Z
Answers.LastEditDate: 2023-01-26T12:46:18.530Z
Answers.UserDisplayName: Steve
Answers.Body: <p>You don't need R or Python for this, it is a very straightforward text parsing problem and you can easily do it in an awk one-liner:</p>
  <pre><code>$ awk 'BEGIN{FS=OFS=&quot;\t&quot;} 
        /^[^#]/ &amp;&amp; $14==&quot;LOW&quot; &amp;&amp; $17=&quot;likely pathogenic&quot; ' file 
  11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW -   -1  likely pathogenic   
  </code></pre>
  <p><strong>Explanation</strong></p>
  <ul>
  <li><code>BEGIN{FS=OFS=&quot;\t&quot;}</code>: the <code>BEGIN{}</code> block is run before processing the file. Here we are setting the field separator (<code>FS</code>) and output field separator (<code>OFS</code>) to a tab so that awk will split each input line on tab characters.</li>
  <li><code>/^[^#]/ &amp;&amp; $14==&quot;LOW&quot; &amp;&amp; $17=&quot;likely pathogenic&quot;</code>: this will evaluate to true if the current line does not start with a <code>#</code>, and the 14th field is <code>LOW</code> and the 17th field is <code>likely pathogenic</code>. In awk, the default action when something evaluates to true is to print the current line, so this results in the expected output.</li>
  </ul>
  <p>Now, if you also want to add the file name as a field, you can do this to have it as the last field:</p>
  <pre><code>$ awk 'BEGIN{FS=OFS=&quot;\t&quot;} /^[^#]/ &amp;&amp; $14==&quot;LOW&quot; &amp;&amp; $17=&quot;likely pathogenic&quot;{$(NF+1)=FILENAME; print} ' file 
  11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW -   -1  likely pathogenic       yourfilename
  </code></pre>
  <p>Next, if your files are compressed and you want to do this to all of them, you could simply do:</p>
  <pre><code>for file in *vep.txt.gz; do
    zcat &quot;$file&quot; | awk -v fname=&quot;${file%.vep.txt.gz}&quot; '
                         BEGIN{FS=OFS=&quot;\t&quot;} 
                         /^[^#]/ &amp;&amp; $14==&quot;LOW&quot; &amp;&amp; $17=&quot;likely pathogenic&quot;{
                            $(NF+1)=FILENAME; 
                            print
                         }'
  done &gt; allfiles.filtered
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T12:26:07.883Z
Answers.Text: I would just use the included `filter_vep` script and merge the output files using AWK if needed. I'm not sure if VEP always produces it's output fields in exactly the same order, so using an array to refer to the IMPACT and CLINSIG fields might be preferable.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T12:38:57.863Z
Answers.Text: @Steve sure, but you already have given a good answer using the included tool, and I personally find it easier and more intuitive to use general scripting tools instead of taking the time to learn tool-specific helper scripts each time. Fair point about the array, but I admit I would be shocked if VEP is so badly written that it changes the order of its fields!
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T13:30:54.227Z
Answers.Text: Assuming it's been run the same way for each file you're probably fine. But if you were to re-run it with a different set of options, the order of the non-default fields (incl. IMPACT and CLIN_SIG) could change. For me, using a recent VEP, the non-default fields actually end up in an "Extra" column consisting of semicolon-separated key=value pairs.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-21T19:42:16.463Z
Answers.Text: @terdon do you know, how I can run these commands in Windows 11?   for file in *vep.txt.gz; do...etc
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-21T19:56:03.417Z
Answers.Text: No, @maj, sorry I don't use Windows. You should be able to do it using windows subshell for windows.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-22T12:53:56.500Z
Answers.Text: @Maj I think terdon means the [Windows Subsystem for Linux, WSL2](https://learn.microsoft.com/en-us/windows/wsl/install)
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-22T19:23:55.967Z
Answers.Text: I did indeed, thanks @Steve
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-25T20:28:14.313Z
Answers.Text: thank you @terdon  The whole code worked correctly except for the command to add a column for the filename as an identifier, do you know what is the problem?
Answers.UserDisplayName: Joman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-25T22:01:21.543Z
Answers.Text: @Maj sorry but I am not sure what you mean. If you are referring to the `FILENAME` part, it is possibel that you are using an awk that doesn't support it. GNU awk (`gawk`), which is the default on Linux systems, supports it but I'm not sure about others.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-01-18T11:21:15.047Z
Answers.LastActivityDate: 2023-01-18T11:21:15.047Z
Answers.UserDisplayName: terdon
Body: <p>I'm new to dealing with big and compressed data so, I have some questions regarding that</p>
  <p>I have 1000 files, each file has the extension (vep.txt.gz) I want to read these files, then filter them. The size of these files is very large. I tried to read these files in Python and R but I am facing a memory error issue.
  I need to read these files and extract some rows from each file for example I have two columns. In the first column, I want the rows with low values in the impact column and the rows with likely pathogenic values from the clinSIG column. Also, I have to add a new column to put the files' names as IDs for each sample.</p>
  <p>So, how can I read a list of big (vep.txt.gz) files in R or Python, then filter the content of the files based on two columns (impact = low) and (clin_SIG = Likely pathogenic) ?</p>
  <p>Should I decompress the files or just read them?</p>
  <p>The following is an example of the files but each file is compressed and all the files have the same header (you can find the original file, which is tab separated, [here][1]):</p>
  <pre><code>## ENSEMBL VARIANT EFFECT PREDICTOR v108.0
  ## Output produced at 2017-03-21 14:51:27
  ## Connected to homo_sapiens_core_108_38 on ensembldb.ensembl.org
  ## Using cache in /homes/user/.vep/homo_sapiens/108_GRCh38
  ## Using API version 108, DB version 108
  ## polyphen version 2.2.2
  ## sift version sift5.2.2
  ## COSMIC version 78
  ## ESP version 20141103
  ## gencode version GENCODE 25
  ## genebuild version 2014-07
  ## HGMD-PUBLIC version 20162
  ## regbuild version 16
  ## assembly version GRCh38.p7
  ## ClinVar version 201610
  ## dbSNP version 147
  ## Column descriptions:
  ## Uploaded_variation : Identifier of uploaded variant
  ## Location : Location of variant in standard coordinate format (chr:start or chr:start-end)
  ## Allele : The variant allele used to calculate the consequence
  ## Gene : Stable ID of affected gene
  ## Feature : Stable ID of feature
  ## Feature_type : Type of feature - Transcript, RegulatoryFeature or MotifFeature
  ## Consequence : Consequence type
  ## cDNA_position : Relative position of base pair in cDNA sequence
  ## CDS_position : Relative position of base pair in coding sequence
  ## Protein_position : Relative position of amino acid in protein
  ## Amino_acids : Reference and variant amino acids
  ## Codons : Reference and variant codon sequence
  ## Existing_variation : Identifier(s) of co-located known variants
  ## Extra column keys:
  ## IMPACT : Subjective impact classification of consequence type
  ## DISTANCE : Shortest distance from variant to transcript
  ## STRAND : Strand of the feature (1/-1)
  ## FLAGS : Transcript quality flags
  
  #Uploaded_variation  Location   Allele  Gene             Feature          Feature_type  Consequence                            cDNA_position  CDS_position  Protein_position  Amino_acids  Codons   Existing_variation  IMPACT    DISTANCE  STRAND  CLINSIG 
  11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000525319  Transcript    missense_variant                       742            716           239               S/I          aGc/aTc  -                   LOW  -         -1      likely pathogenic  
  11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000534381  Transcript    downstream_gene_variant                -              -             -                 -            -        -                   MODIFIER  1674      -1      -
  11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000529055  Transcript    downstream_gene_variant                -              -             -                 -            -        -                   HIGH  134       -1      -
  11_224585_G/A        11:224585  A       ENSG00000142082  ENST00000529937  Transcript    intron_variant,NMD_transcript_variant  -              -             -                 -            -        -                   MODIFIER  -         -1      -
  </code></pre>
  <p>From the previous file, the expected output will be one row which is the first row because the impact is low &amp; CLINSIG is likely pathogenic so i want to find all rows with this condition for 1000 files</p>
  <pre><code>ID  #Uploaded_variation  Location   Allele  Gene    Feature Feature_type  Consequence   cDNA_position   CDS_position    Protein_position  Amino_acids  Codons   Existing_variation  IMPACT    DISTANCE  STRAND  CLINSIG 
  M1457.vep   11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW  -  -1  likely pathogenic
  
  
  </code></pre>
  <p>Thanks in advance</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-16T00:06:15.610Z
Comments.Text: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
Comments.UserDisplayName: Community
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-16T10:44:39.390Z
Comments.Text: What is in `vep.txt.gz` exactly? How did you try to read this in Python? What, exactly, do you want to filter and how? What are your filtering criteria? If you show us an example input (you can decompress the file and post a few lines here) and the output you want from it, then we can help you but without any information about the file or its format, we cannot really do much.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-17T21:20:58.117Z
Comments.Text: I've edited the question hope it's now more clear and also hope you can help me @terdon
Comments.UserDisplayName: Joman
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-18T00:09:12.223Z
Comments.Text: Thanks but I'm afraid we can't do anything with an image. We need to have an example that we can use to test our solutions. Please add the file as text (paste it into your question, and use the `{}` button to format it as code) and also show us the output you expect from it so we can check if we get it right.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-18T06:52:15.587Z
Comments.Text: now you will find the file @terdon
Comments.UserDisplayName: Joman
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-18T09:50:19.207Z
Comments.Text: Should we assume the file is tab-separated?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-18T10:08:54.180Z
Comments.Text: yeah, the format is tab-separated @terdon
Comments.UserDisplayName: Joman
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-18T11:09:45.210Z
Comments.Text: Hang on, what output do you want? You say "the expected output will be [. . . ] the first row" but earlier you said "I have to add a new column to put the files' names as IDs for each sample". So what do you actually want as output? Please ***edit*** your question and add your expected output.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-01-18T12:14:13.293Z
Comments.Text: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/142143/discussion-between-maj-and-terdon).
Comments.UserDisplayName: Joman
Tags: <r><python><filtering><vep>
Title: dealing with a list of VEP files
ThreadText: # dealing with a list of VEP files
  
  (Asked by: username_0 on 2023-01-15)
  
  I'm new to dealing with big and compressed data so, I have some questions regarding that
  
  I have 1000 files, each file has the extension (vep.txt.gz) I want to read these files, then filter them. The size of these files is very large. I tried to read these files in Python and R but I am facing a memory error issue.
  I need to read these files and extract some rows from each file for example I have two columns. In the first column, I want the rows with low values in the impact column and the rows with likely pathogenic values from the clinSIG column. Also, I have to add a new column to put the files' names as IDs for each sample.
  
  So, how can I read a list of big (vep.txt.gz) files in R or Python, then filter the content of the files based on two columns (impact = low) and (clin\_SIG = Likely pathogenic) ?
  
  Should I decompress the files or just read them?
  
  The following is an example of the files but each file is compressed and all the files have the same header (you can find the original file, which is tab separated, [here][1]):
  
  ```
  ## ENSEMBL VARIANT EFFECT PREDICTOR v108.0
  ## Output produced at 2017-03-21 14:51:27
  ## Connected to homo_sapiens_core_108_38 on ensembldb.ensembl.org
  ## Using cache in /homes/user/.vep/homo_sapiens/108_GRCh38
  ## Using API version 108, DB version 108
  ## polyphen version 2.2.2
  ## sift version sift5.2.2
  ## COSMIC version 78
  ## ESP version 20141103
  ## gencode version GENCODE 25
  ## genebuild version 2014-07
  ## HGMD-PUBLIC version 20162
  ## regbuild version 16
  ## assembly version GRCh38.p7
  ## ClinVar version 201610
  ## dbSNP version 147
  ## Column descriptions:
  ## Uploaded_variation : Identifier of uploaded variant
  ## Location : Location of variant in standard coordinate format (chr:start or chr:start-end)
  ## Allele : The variant allele used to calculate the consequence
  ## Gene : Stable ID of affected gene
  ## Feature : Stable ID of feature
  ## Feature_type : Type of feature - Transcript, RegulatoryFeature or MotifFeature
  ## Consequence : Consequence type
  ## cDNA_position : Relative position of base pair in cDNA sequence
  ## CDS_position : Relative position of base pair in coding sequence
  ## Protein_position : Relative position of amino acid in protein
  ## Amino_acids : Reference and variant amino acids
  ## Codons : Reference and variant codon sequence
  ## Existing_variation : Identifier(s) of co-located known variants
  ## Extra column keys:
  ## IMPACT : Subjective impact classification of consequence type
  ## DISTANCE : Shortest distance from variant to transcript
  ## STRAND : Strand of the feature (1/-1)
  ## FLAGS : Transcript quality flags
  
  #Uploaded_variation  Location   Allele  Gene             Feature          Feature_type  Consequence                            cDNA_position  CDS_position  Protein_position  Amino_acids  Codons   Existing_variation  IMPACT    DISTANCE  STRAND  CLINSIG 
  11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000525319  Transcript    missense_variant                       742            716           239               S/I          aGc/aTc  -                   LOW  -         -1      likely pathogenic  
  11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000534381  Transcript    downstream_gene_variant                -              -             -                 -            -        -                   MODIFIER  1674      -1      -
  11_224088_C/A        11:224088  A       ENSG00000142082  ENST00000529055  Transcript    downstream_gene_variant                -              -             -                 -            -        -                   HIGH  134       -1      -
  11_224585_G/A        11:224585  A       ENSG00000142082  ENST00000529937  Transcript    intron_variant,NMD_transcript_variant  -              -             -                 -            -        -                   MODIFIER  -         -1      -
  
  ```
  
  From the previous file, the expected output will be one row which is the first row because the impact is low & CLINSIG is likely pathogenic so i want to find all rows with this condition for 1000 files
  
  ```
  ID  #Uploaded_variation  Location   Allele  Gene    Feature Feature_type  Consequence   cDNA_position   CDS_position    Protein_position  Amino_acids  Codons   Existing_variation  IMPACT    DISTANCE  STRAND  CLINSIG 
  M1457.vep   11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW  -  -1  likely pathogenic
  
  ```
  
  Thanks in advance
  
  **Tags:** r, python, filtering, vep
  
  > **Comments:**
  > **username_1**: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.  
  > **username_2**: What is in `vep.txt.gz` exactly? How did you try to read this in Python? What, exactly, do you want to filter and how? What are your filtering criteria? If you show us an example input (you can decompress the file and post a few lines here) and the output you want from it, then we can help you but without any information about the file or its format, we cannot really do much.  
  > **username_0**: I've edited the question hope it's now more clear and also hope you can help me @username_2  
  > **username_2**: Thanks but I'm afraid we can't do anything with an image. We need to have an example that we can use to test our solutions. Please add the file as text (paste it into your question, and use the `{}` button to format it as code) and also show us the output you expect from it so we can check if we get it right.  
  > **username_0**: now you will find the file @username_2  
  > **username_2**: Should we assume the file is tab-separated?  
  > **username_0**: yeah, the format is tab-separated @username_2  
  > **username_2**: Hang on, what output do you want? You say "the expected output will be [. . . ] the first row" but earlier you said "I have to add a new column to put the files' names as IDs for each sample". So what do you actually want as output? Please ***edit*** your question and add your expected output.  
  > **username_0**: Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/142143/discussion-between-maj-and-username_2).  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_3 on 2023-01-16. Score: 3
  
  For these requirements, it might be possible to simply use [filter\_vep](https://asia.ensembl.org/info/docs/tools/vep/script/vep_filter.html) included in the VEP package for this. Note that multiple `--filter` flags are ANDed together (i.e. all filters must pass for a line to be printed). For example:
  
  ```
  for i in *.vep.txt.gz; do
      filter_vep \
          -i "${i}" \
          -o "${i%.vep.txt.gz}.filtered.vep.txt" \
          --filter "IMPACT is LOW" \
          --filter "CLIN_SIG matches likely_pathogenic"
  done
  
  ```
  
  If you then wanted to merge the output and include the filename, you could just use [awk](/questions/tagged/awk "show questions tagged 'awk'"). With the following in a file called `merge_filtered_vep_txt.awk`:
  
  ```
  BEGIN {
      FS=OFS="\t"
      f=1
  }
  
  /^##/ {
      next
  }
  
  sub(/^#/, "") && f {
      print "ID", $0
      f=0
      next
  }
  
  {
      print FILENAME, $0
  }
  
  ```
  
  Run using:
  
  ```
  awk -f merge_filtered_vep_txt.awk /path/to/*.filtered.vep.txt
  
  ```
  
  Python can also open gzip-compressed files for reading (and writing), using the [gzip](https://docs.python.org/3/library/gzip.html) module in the standard library. These can be parsed line-by-line to avoid memory errors. Here's one way using the [csv](https://docs.python.org/3/library/csv.html) module (also in the standard library):
  
  ```
  import argparse
  import contextlib
  import csv
  import gzip
  import pathlib
  import sys
  
  ```
  
  ```
  def get_argument_parser():
  
      parser = argparse.ArgumentParser(add_help=False)
  
      group = parser.add_argument_group('input options')
      group.add_argument(
          "-i",
          "--input",
          dest='glob_pattern',
          type=str,
          default="**/*.vep.txt.gz",
          metavar="STR",
          help="The relative glob pattern (default: '%(default)s')",
      )
      group.add_argument(
          "-d",
          "--dir",
          dest='glob_path',
          type=pathlib.Path,
          metavar="DIR",
          default='.',
          help="Glob files in DIR (default: '%(default)s')",
      )
  
      group = parser.add_argument_group('output options')
      group.add_argument(
          "-o",
          "--out",
          type=pathlib.Path,
          metavar="FILE",
          default='-',
          help="Write the filtered output to FILE (default: stdout)",
      )
      group.add_argument(
          "-f",
          "--force",
          action='store_true',
          help="Overwrite the output file if it exists",
      )
  
      group = parser.add_argument_group('additional options')
      group.add_argument(
          "-h",
          "--help",
          action="help",
          help="Show this help message and exit",
      )
  
      return parser
  
  ```
  
  ```
  def main():
  
      parser = get_argument_parser()
      args = parser.parse_args()
  
      with contextlib.ExitStack() as stack:
          if args.out.name != '-':
              out = stack.enter_context(open(args.out, 'wt' if args.force else 'xt'))
          else:
              out = sys.stdout
  
          # Get a list of files from the directory specified by args.glob_path
  
          file_list = sorted(args.glob_path.glob(args.glob_pattern))
  
          # Build up a unique list of header fields from each input file
  
          header_cols = ['ID']
          for gz_file in file_list:
              with gzip.open(gz_file, 'rt') as infile:
                  for line in infile:
                      if not line.startswith('##'):
                          break
  
                  assert line.startswith('#')
  
                  for field in line[1:].rstrip().split('\t'):
                      if field not in header_cols:
                          header_cols.append(field)
  
          # Create a csv.DictWriter() and write the header fields
  
          writer = csv.DictWriter(out, fieldnames=header_cols, delimiter='\t')
          writer.writeheader()
  
          # Now re-open each file in the arguments list for reading
  
          for gz_file in file_list:
              with gzip.open(gz_file, 'rt') as infile:
  
                  # Skip the info header lines
                  line = next(infile)
                  while line.startswith('##'):
                      line = next(infile)
  
                  # Get the column headers
                  fieldnames = line.rstrip().removeprefix('#').split('\t')
  
                  # Parse each input file and store the desired rows in memory
                  reader = csv.DictReader(infile, fieldnames=fieldnames, delimiter='\t')
                  rows = []
                  for row in reader:
  
                      impact = row.get('IMPACT')
                      clinsig = row.get('CLINSIG')
  
                      if impact == 'LOW' and clinsig == 'likely pathogenic':
                          rows.append({'ID': gz_file, **row})
  
                  # Write the list of rows to the output file
                  writer.writerows(rows)
  
  ```
  
  ```
  if __name__ == '__main__':
      main()
  
  ```
  
  With the above in a file called `filter_vep_txt.py`:
  
  ```
  usage: filter_vep_txt.py [-i STR] [-d DIR] [-o FILE] [-f] [-h]
  
  input options:
    -i STR, --input STR  The relative glob pattern (default: '**/*.vep.txt.gz')
    -d DIR, --dir DIR    Glob files in DIR (default: '.')
  
  output options:
    -o FILE, --out FILE  Write the filtered output to FILE (default: stdout)
    -f, --force          Overwrite the output file if it exists
  
  additional options:
    -h, --help           Show this help message and exit
  
  ```
  
  > **Comments:**
  > **username_0**: thank you in python code, what is the line code for determining the list of files of vep.txt.gz? @username_3  
  > **username_4**: vep.txt.gz? ----> I suspect is just one file named vep.txt see: https://stackoverflow.com/questions/31028815/how-to-unzip-gz-file-using-python ; https://unix.stackexchange.com/questions/93139/can-i-zip-an-entire-folder-using-gzip  
  > **username_3**: @Maj Maybe I haven't understood correctly, but I think that's really a question for you. I.e. what files would you like to process? How would you like to come up with the list? Would you prefer some sort of glob pattern or would you prefer to supply a file of filenames, one per line? Perhaps you could [edit](https://bioinformatics.stackexchange.com/posts/20370/edit) your question to include some details?  
  > **username_0**: I've edited the question hope it's now more clear and also hope you can help me @username_3  
  > **username_3**: @Maj Please see the updated code above. It's not clear if each file contains the same header fields, so the above reads the list of input files twice to ensure we get a unique list of header columns. Run using `python filter_vep_txt.py -o out.txt *.vep.txt.gz` or similar.  
  > **username_0**: thanks a lot @username_3, sorry but I can not know how your code finds the files. is it for one files or for all?  
  > **username_0**: What is the command to locate the files? @username_3  
  > **username_3**: @Maj You'll need to supply one or more files as positional arguments on the command line. If all of your files are in the current directory, you should be able to just use the above glob pattern. If these are under multiple directories, you could use: `python filter_vep_txt.py -o out.txt /path/to/files/*.vep.txt.gz /path/to/more/files/*.vep.txt.gz` for example.  
  > **username_0**: ok thank you @username_3 now it is clear, I saved the first 4 blocks of python code in filter_vep_txt.py but the last block code you wrote how can I run it? sorry I'm new in this field  
  > **username_3**: @Maj the last block is just the usage. You can print it using: `python filter_vep_txt.py -h`  
  > **username_3**: @Maj In my code above, you may need to change `'Likely pathogenic'` to `'likely pathogenic'`. Is your CLINSIG column actually a comma-separated list? I'm just curious, but what options did you use to run VEP? Maybe you could include this in your question?  
  > **username_0**: hi @username_3, when I run the code by python filter_vep_txt.py -o out.txt *.vep.txt.gz, It shows me this error>> OSError: [Errno 22] Invalid argument: '*.vep.txt.gz'.     NOTE (all the files are in the same directory)  
  > **username_3**: @Maj You just have a typo in your command, use: `python filter_vep_txt.py -o out.txt *.vep.txt.gz`. Note the star/asterisk in the glob pattern.  
  > **username_0**: @username_3 yeah I tried with one file and with all files but still I get the same error ((OSError: [Errno 22] Invalid argument: '*.vep.txt.gz'))  
  > **username_0**: Traceback (most recent call last):   File "C:\Users \Desktop\files\filter_vep_txt.py", line 123, in <module>     main()   File "C:\Users\Desktop\files\filter_vep_txt.py", line 75, in main     with gzip.open(vep_file, 'rt') as infile:   File "C:\Program Files\WindowsApps\PythonSoftwareFoundation.Python.3.10_3.10.2544.0_x64__qbz5n2kfra8p0\lib\gzip.py", line 58, in open     binary_file = GzipFile(filename, gz_mode, compresslevel)   File "C:\Program Files\WindowsApps\PythonSoftwareFoundation.Python.3.10_3.10.2544.0_x64__qbz5n2kfra8p0\lib\gzip.py", line 174, in __init__ @username_3  
  > **username_0**: fileobj = self.myfileobj = builtins.open(filename, mode or 'rb') OSError: [Errno 22] Invalid argument: '*.vep.txt.gz' @username_3 sorry for many questions  
  > **username_3**: @Maj I think the problem is just that your shell can't expand the glob pattern. I think the simpler solution is to just have the script find the files. Please see the changes above. Use like: `python filter_vep_txt.py -d /path/to/files -o out.txt`, and it will try to glob the files recursively.  
  > **username_0**: @username_3 the code is working now, but when I open the out.txt file, I found only the word (ID) and I did not find any other information, I ran this command  (python filter_vep_txt.py -d *.vep.txt.gz -o out.txt) in the same files path  
  > **username_3**: @Maj Use `-i` instead of `-d`: `python filter_vep_txt.py -i '*.vep.txt.gz' -o out.txt`. Be sure to use single quotes to prevent it from being expanded by your shell. Failing that, add a `print(line)` before the `assert` statement and post the output. I'm not expecting VEP output files to contain blank lines in the header, but maybe your output is slightly different.  
  > **username_0**: @username_3 thanks a lot for your reply... I get this error and when I open out.txt file I found only the header of files (columns names) ..... Traceback (most recent call last):   File "C:\Users\t\Desktop\files\filter_vep_txt.py", line 122, in <module>     main()   File "C:\Users\t\Desktop\files\filter_vep_txt.py", line 116, in main     rows.append({'ID': vep_file, **row}) NameError: name 'vep_file' is not defined. Did you mean: 'gz_file'?  
  > **username_3**: @Maj Sorry, typo. Yes I did indeed mean 'gz_file'.  
  
  ---
  
  ### Answer by username_2 on 2023-01-18. Score: 1
  
  You don't need R or Python for this, it is a very straightforward text parsing problem and you can easily do it in an awk one-liner:
  
  ```
  $ awk 'BEGIN{FS=OFS="\t"} 
        /^[^#]/ && $14=="LOW" && $17="likely pathogenic" ' file 
  11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW -   -1  likely pathogenic   
  
  ```
  
  **Explanation**
  
  * `BEGIN{FS=OFS="\t"}`: the `BEGIN{}` block is run before processing the file. Here we are setting the field separator (`FS`) and output field separator (`OFS`) to a tab so that awk will split each input line on tab characters.
  * `/^[^#]/ && $14=="LOW" && $17="likely pathogenic"`: this will evaluate to true if the current line does not start with a `#`, and the 14th field is `LOW` and the 17th field is `likely pathogenic`. In awk, the default action when something evaluates to true is to print the current line, so this results in the expected output.
  
  Now, if you also want to add the file name as a field, you can do this to have it as the last field:
  
  ```
  $ awk 'BEGIN{FS=OFS="\t"} /^[^#]/ && $14=="LOW" && $17="likely pathogenic"{$(NF+1)=FILENAME; print} ' file 
  11_224088_C/A   11:224088   A   ENSG00000142082 ENST00000525319 Transcript  missense_variant    742 716 239 S/I aGc/aTc -   LOW -   -1  likely pathogenic       yourfilename
  
  ```
  
  Next, if your files are compressed and you want to do this to all of them, you could simply do:
  
  ```
  for file in *vep.txt.gz; do
    zcat "$file" | awk -v fname="${file%.vep.txt.gz}" '
                         BEGIN{FS=OFS="\t"} 
                         /^[^#]/ && $14=="LOW" && $17="likely pathogenic"{
                            $(NF+1)=FILENAME; 
                            print
                         }'
  done > allfiles.filtered
  
  ```
  
  > **Comments:**
  > **username_3**: I would just use the included `filter_vep` script and merge the output files using AWK if needed. I'm not sure if VEP always produces it's output fields in exactly the same order, so using an array to refer to the IMPACT and CLINSIG fields might be preferable.  
  > **username_2**: @username_3 sure, but you already have given a good answer using the included tool, and I personally find it easier and more intuitive to use general scripting tools instead of taking the time to learn tool-specific helper scripts each time. Fair point about the array, but I admit I would be shocked if VEP is so badly written that it changes the order of its fields!  
  > **username_3**: Assuming it's been run the same way for each file you're probably fine. But if you were to re-run it with a different set of options, the order of the non-default fields (incl. IMPACT and CLIN_SIG) could change. For me, using a recent VEP, the non-default fields actually end up in an "Extra" column consisting of semicolon-separated key=value pairs.  
  > **username_0**: @username_2 do you know, how I can run these commands in Windows 11?   for file in *vep.txt.gz; do...etc  
  > **username_2**: No, @maj, sorry I don't use Windows. You should be able to do it using windows subshell for windows.  
  > **username_3**: @Maj I think username_2 means the [Windows Subsystem for Linux, WSL2](https://learn.microsoft.com/en-us/windows/wsl/install)  
  > **username_2**: I did indeed, thanks @username_3  
  > **username_0**: thank you @username_2  The whole code worked correctly except for the command to add a column for the filename as an identifier, do you know what is the problem?  
  > **username_2**: @Maj sorry but I am not sure what you mean. If you are referring to the `FILENAME` part, it is possibel that you are using an awk that doesn't support it. GNU awk (`gawk`), which is the default on Linux systems, supports it but I'm not sure about others.  
  
  ---
  
  