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

CC BY-SA 4.0
2023-01-16T07:06:43.247Z
thank you in python code, what is the line code for determining the list of files of vep.txt.gz? @Steve
Joman
CC BY-SA 4.0
2023-01-16T10:51:17.267Z
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
pippo1980
CC BY-SA 4.0
2023-01-16T11:45:20.700Z
@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?
Steve
CC BY-SA 4.0
2023-01-17T21:21:25.127Z
I've edited the question hope it's now more clear and also hope you can help me @Steve
Joman
CC BY-SA 4.0
2023-01-18T04:28:59.963Z
@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.
Steve
CC BY-SA 4.0
2023-01-18T06:55:40.070Z
thanks a lot @Steve, sorry but I can not know how your code finds the files. is it for one files or for all?
Joman
CC BY-SA 4.0
2023-01-18T07:07:52.800Z
What is the command to locate the files? @Steve
Joman
CC BY-SA 4.0
2023-01-18T08:43:39.303Z
@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.
Steve
CC BY-SA 4.0
2023-01-18T10:56:32.380Z
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
Joman
CC BY-SA 4.0
2023-01-18T11:44:29.533Z
@Maj the last block is just the usage. You can print it using: `python filter_vep_txt.py -h`
Steve
CC BY-SA 4.0
2023-01-18T13:32:47.880Z
@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?
Steve
CC BY-SA 4.0
2023-01-22T07:08:56.537Z
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)
Joman
CC BY-SA 4.0
2023-01-22T12:32:59.773Z
@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.
Steve
CC BY-SA 4.0
2023-01-22T19:22:47.387Z
@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'))
Joman
CC BY-SA 4.0
2023-01-22T19:25:48.700Z
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
Joman
CC BY-SA 4.0
2023-01-22T19:26:46.187Z
fileobj = self.myfileobj = builtins.open(filename, mode or 'rb')
OSError: [Errno 22] Invalid argument: '*.vep.txt.gz'
@Steve sorry for many questions
Joman
CC BY-SA 4.0
2023-01-23T12:57:43.697Z
@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.
Steve
CC BY-SA 4.0
2023-01-25T20:19:16.377Z
@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
Joman
CC BY-SA 4.0
2023-01-26T02:56:30.220Z
@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.
Steve
CC BY-SA 4.0
2023-01-26T11:35:32.143Z
@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'?
Joman
CC BY-SA 4.0
2023-01-26T12:47:49.840Z
@Maj Sorry, typo. Yes I did indeed mean 'gz_file'.
Steve
CC BY-SA 4.0
2023-01-16T01:06:45.570Z
2023-01-26T12:46:18.530Z
2023-01-26T12:46:18.530Z
Steve
<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>

CC BY-SA 4.0
2023-01-18T12:26:07.883Z
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.
Steve
CC BY-SA 4.0
2023-01-18T12:38:57.863Z
@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!
terdon
CC BY-SA 4.0
2023-01-18T13:30:54.227Z
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.
Steve
CC BY-SA 4.0
2023-01-21T19:42:16.463Z
@terdon do you know, how I can run these commands in Windows 11?   for file in *vep.txt.gz; do...etc
Joman
CC BY-SA 4.0
2023-01-21T19:56:03.417Z
No, @maj, sorry I don't use Windows. You should be able to do it using windows subshell for windows.
terdon
CC BY-SA 4.0
2023-01-22T12:53:56.500Z
@Maj I think terdon means the [Windows Subsystem for Linux, WSL2](https://learn.microsoft.com/en-us/windows/wsl/install)
Steve
CC BY-SA 4.0
2023-01-22T19:23:55.967Z
I did indeed, thanks @Steve
terdon
CC BY-SA 4.0
2023-01-25T20:28:14.313Z
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?
Joman
CC BY-SA 4.0
2023-01-25T22:01:21.543Z
@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.
terdon
CC BY-SA 4.0
2023-01-18T11:21:15.047Z
2023-01-18T11:21:15.047Z
terdon
<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>

CC BY-SA 4.0
2023-01-16T00:06:15.610Z
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.
Community
CC BY-SA 4.0
2023-01-16T10:44:39.390Z
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.
terdon
CC BY-SA 4.0
2023-01-17T21:20:58.117Z
I've edited the question hope it's now more clear and also hope you can help me @terdon
Joman
CC BY-SA 4.0
2023-01-18T00:09:12.223Z
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.
terdon
CC BY-SA 4.0
2023-01-18T06:52:15.587Z
now you will find the file @terdon
Joman
CC BY-SA 4.0
2023-01-18T09:50:19.207Z
Should we assume the file is tab-separated?
terdon
CC BY-SA 4.0
2023-01-18T10:08:54.180Z
yeah, the format is tab-separated @terdon
Joman
CC BY-SA 4.0
2023-01-18T11:09:45.210Z
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.
terdon
CC BY-SA 4.0
2023-01-18T12:14:13.293Z
Let us [continue this discussion in chat](https://chat.stackexchange.com/rooms/142143/discussion-between-maj-and-terdon).
Joman
<r><python><filtering><vep>
dealing with a list of VEP files
# 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 @terdon  
> **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 @terdon  
> **username_2**: Should we assume the file is tab-separated?  
> **username_0**: yeah, the format is tab-separated @terdon  
> **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-terdon).  

---

## 2 answers

### Answer by username_3 on 2023-01-16. Score: 3

For these requirements, it might be possible to simply use filter_vep 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. 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 module in the standard library. These can be parsed line-by-line to avoid memory errors. Here's one way using the csv 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? @Steve  
> **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 @Steve  
> **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 @Steve, 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? @Steve  
> **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 @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  
> **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 @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)  
> **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**: @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'))  
> **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__ @Steve  
> **username_0**: fileobj = self.myfileobj = builtins.open(filename, mode or 'rb') OSError: [Errno 22] Invalid argument: '*.vep.txt.gz' @Steve 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**: @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  
> **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**: @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'?  
> **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**: @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!  
> **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**: @terdon 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 terdon means the [Windows Subsystem for Linux, WSL2](https://learn.microsoft.com/en-us/windows/wsl/install)  
> **username_2**: I did indeed, thanks @Steve  
> **username_0**: 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?  
> **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.  

---

