dealing with a list of VEP files

Question

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

Answer 1

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

Answer 2

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