3
$\begingroup$

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

$\endgroup$
9
  • 2
    $\begingroup$ 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. $\endgroup$
    – Community Bot
    Jan 16, 2023 at 0:06
  • 1
    $\begingroup$ 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. $\endgroup$
    – terdon
    Jan 16, 2023 at 10:44
  • $\begingroup$ I've edited the question hope it's now more clear and also hope you can help me @terdon $\endgroup$
    – Joman
    Jan 17, 2023 at 21:20
  • $\begingroup$ 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. $\endgroup$
    – terdon
    Jan 18, 2023 at 0:09
  • $\begingroup$ now you will find the file @terdon $\endgroup$
    – Joman
    Jan 18, 2023 at 6:52

2 Answers 2

3
$\begingroup$

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 . 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
$\endgroup$
21
  • $\begingroup$ thank you in python code, what is the line code for determining the list of files of vep.txt.gz? @Steve $\endgroup$
    – Joman
    Jan 16, 2023 at 7:06
  • $\begingroup$ vep.txt.gz? ----> I suspect is just one file named vep.txt see: stackoverflow.com/questions/31028815/… ; unix.stackexchange.com/questions/93139/… $\endgroup$
    – pippo1980
    Jan 16, 2023 at 10:51
  • $\begingroup$ @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 your question to include some details? $\endgroup$
    – Steve
    Jan 16, 2023 at 11:45
  • $\begingroup$ I've edited the question hope it's now more clear and also hope you can help me @Steve $\endgroup$
    – Joman
    Jan 17, 2023 at 21:21
  • 1
    $\begingroup$ @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. $\endgroup$
    – Steve
    Jan 26, 2023 at 2:56
1
$\begingroup$

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
$\endgroup$
9
  • 1
    $\begingroup$ 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. $\endgroup$
    – Steve
    Jan 18, 2023 at 12:26
  • $\begingroup$ @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! $\endgroup$
    – terdon
    Jan 18, 2023 at 12:38
  • $\begingroup$ 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. $\endgroup$
    – Steve
    Jan 18, 2023 at 13:30
  • $\begingroup$ @terdon do you know, how I can run these commands in Windows 11? for file in *vep.txt.gz; do...etc $\endgroup$
    – Joman
    Jan 21, 2023 at 19:42
  • $\begingroup$ No, @maj, sorry I don't use Windows. You should be able to do it using windows subshell for windows. $\endgroup$
    – terdon
    Jan 21, 2023 at 19:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.