1
$\begingroup$

I am relatively new to genome analysis and would like to compare 30+ vcf files, 1 sample per file. After filtering there are about ~20 variants per vcf file.

I would like to evaluate to mutational landscape between samples and for each sample as well.

One option that I found is a mutation heatmap where each row is a sample and each column a position in the genome. Is this the preferred way of displaying this type of information and how would I include gene information?

At this point I am interested in getting an overview first of the affected positions and then dig deeper using e.g. zygosity.

Thank you very much for helping out!

2nd edit I would also like to evaluate if there is a bias towards certain genes. For example are certain oncogenes mutated more frequently in this cohort or is the accumulation random and unique within each sample?

Best, M

Improve this question
$\endgroup$
5
  • 1
    $\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
    Feb 28 at 11:42
  • 3
    $\begingroup$ What exactly do you want to compare? Simply whether a variant was called in the different samples? Would you want to take zygosity into account (is it a match if the same variant is found heterozygous in one sample and homozygoyus in another)? Do you want to compare the genes which had variants found in them and how they segregate in your samples? What exactly does "visualize a variant" mean? $\endgroup$
    – terdon ♦
    Feb 28 at 13:57
  • $\begingroup$ Thank you very much for the comments and I updated the question. Currently I am searching for a good way to display the position of a variant in the genome for each sample and include e.g. gene annotation. @terdon Then I would like to compare the affected genes, for example are the mutations biased towards oncogenes. $\endgroup$
    – Macintosh
    Mar 2 at 14:28
  • $\begingroup$ Thanks for the extra context. I'll reiterate that you cannot assume that the lack of a variant in a PVCF file implies that the sequence is homozygous reference at that locus. That locus could be poorly sequenced and therefore be a "no call". You need to "joint call" the original GVCF files to be confident. OK, but your question: I don't work in cancer, so I'm unfamiliar with what you're trying to do. It sounds like you have ~60 genotypes. Both R and Python have libraries for reading VCFs and visualization. The R library ggplot has gemo_tile which can create heat maps. $\endgroup$
    – Daniel King
    Mar 3 at 1:01
  • $\begingroup$ It's a bit hard for me to advise about your second edit because it involves questions of science (I'm a software engineer). Do you have a particular statistical distribution of variation within genes that you think is normal? Again, these questions might be obvious to someone working in cancer, but I'm unfamiliar with standard practice there. Perhaps a two-sample Student's t-test is appropriate? With just 60 genotypes, I'd just manually copy each one from the VCF and into an R script and use the R Student's t-Test function: stat.ethz.ch/R-manual/R-devel/library/stats/html/t.test.html . $\endgroup$
    – Daniel King
    Mar 3 at 1:10

1 Answer 1

Reset to default
3
$\begingroup$

Single sample PVCF files are generally not combinable in the way you describe. PVCF stands for "project" VCF. As far as I know, PVCF isn't defined anywhere but, it refers to a multi-sample VCF that does not contain reference blocks. It only contains lines for which at least one sample has at least one alternate allele. We distinguish this from a GVCF which is a single-sample VCF that does contain reference blocks. We can't combine two PVCFs in a principled way because some variant might only appear in one PVCF. We lack the reference block information necessary to decide if a sample in the other VCF is 0/0 or "no call".

Do all the VCF files have the same set of variants? If all the VCF files have the same set of variants, you might try using Hail (disclaimer: I'm a Hail maintainer):

import hail as hl

vcf_files = ['sample1.vcf.bgz', ...]
mts = [hl.import_vcf(f) for f in vcf_files]
mt = hl.MatrixTable.union_cols(mts)  # assumes each file has the same set of variants
mt.write('combined.mt')  # save to Hail's native format so the following code is fast

mt = hl.read_matrix_table('combined.mt')
mt.show(n_rows=5, n_cols=10)  # see a table of the first 5 rows and 10 cols

mt = hl.variant_qc(mt)  # adds a bunch of variant QC annotations
mt.describe()  # inspect available annotations

You mentioned "genomic distribution". Are you referring to the allele frequency in this sample? You can see a genome-wide scatter plot of reference allele frequency like this:

from hail.ggplot import *
fig = ggplot(mt) + geom_point(aes(x=mt.locus, y=mt.variant_qc.AF[0]))
fig.show()
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for your answer. Each vcf file is derived from a different sample and I expect minimal overlap between the called mutations. I updated my post, but in brief, I would like to summarise the mutational burden between samples and evaluate if there is a bias towards certain genes. for example are certain oncogenes mutated more frequently in this cohort or is the accumulation more random? $\endgroup$
    – Macintosh
    Mar 2 at 14:36

Your Answer

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

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