Best practice to perform a functional analysis on a set of VCFs from a cohort of patients

Question

I am confused about how to perform a functional analysis after a variant calling on a cohort of patients. I have anotated all vcfs and filter them in order to get only those deleterious variants.

I was thinking on merging those vcf: I tried to merge them with bcftools merge and also using VariantAnnotation and transforming each sample vcf into a dataframe containing the variant lines to then row bind all those dataframes together.

As an example, Ill show some fields of teh first row of a dataframe:

CHROM   START REF ALT FILTER DP  GENEID  SIFT
  17     1000   T   A   PASS   62   TP53    0

Then, I would extract the genes affected by those variants and perform a gene set enrichment analysis, a pathway analysis similar to those performed after an RNA-seq analysis (with enrichR, for example). But I am confused about the viability of this steps. For example, I would like that, during the pathways analysis, the program I use not just reads the list of mutated genes in my cohort, but also takes into account that, if a gene is mutated in more than one sample, it should wheight more in the analysis than a gene which appears in just one sample.

Answer 1

As a simple approach you may want to try the "fuzzy" list implemented in EnrichR itself. There is some information about this in the Help Center under the "What is a gene set" heading.

This way you can code the number of data sets a gene appears as a value between 0 and 1. In the Help Center this is used to convert fold changes or p-values to a membership score but I don't see why this would not also work for the number of data sets a SNP is observed in.

Since it is relatively easy to implement this can give you an first impression of what to expect. It then may be combined with analyzing each data set independently and then comparing the overlap.