Answers.Body: <p>Single sample PVCF files are generally not combinable in the way you describe. PVCF stands for &quot;project&quot; 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 &quot;no call&quot;.</p>
  <p>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 <a href="https://hail.is/" rel="nofollow noreferrer">Hail</a> (disclaimer: I'm a Hail maintainer):</p>
  <pre class="lang-py prettyprint-override"><code>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
  </code></pre>
  <p>You mentioned &quot;genomic distribution&quot;. Are you referring to the allele frequency in this sample? You can see a genome-wide scatter plot of reference allele frequency like this:</p>
  <pre class="lang-py prettyprint-override"><code>from hail.ggplot import *
  fig = ggplot(mt) + geom_point(aes(x=mt.locus, y=mt.variant_qc.AF[0]))
  fig.show()
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-03-02T14:36:21.293Z
Answers.Text: 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?
Answers.UserDisplayName: Macintosh
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-02-28T17:12:49.200Z
Answers.LastActivityDate: 2023-03-01T18:45:57.363Z
Answers.LastEditDate: 2023-03-01T18:45:57.363Z
Answers.UserDisplayName: Daniel King
Body: <p>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.</p>
  <p>I would like to evaluate to mutational landscape between samples and for each sample as well.</p>
  <p>One option that I found is a <em><strong>mutation heatmap</strong></em> 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?</p>
  <p>At this point I am interested in getting an overview first of the affected positions and then dig deeper using e.g. zygosity.</p>
  <p>Thank you very much for helping out!</p>
  <p><strong>2nd edit</strong>
  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?</p>
  <p>Best,
  M</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-02-28T11:42:31.833Z
Comments.Text: 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.
Comments.UserDisplayName: Community
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-02-28T13:57:54.337Z
Comments.Text: 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?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-02T14:28:53.673Z
Comments.Text: 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.
Comments.UserDisplayName: Macintosh
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-03T01:01:40.873Z
Comments.Text: 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.
Comments.UserDisplayName: Daniel King
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-03T01:10:47.490Z
Comments.Text: 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: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/t.test.html .
Comments.UserDisplayName: Daniel King
Tags: <vcf><variant-calling><visualization>
Title: variant visualisation of VCF files for 30+ samples
ThreadText: # variant visualisation of VCF files for 30+ samples
  
  (Asked by: username_0 on 2023-02-28)
  
  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
  
  **Tags:** vcf, variant-calling, visualization
  
  > **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 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?  
  > **username_0**: 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. @username_2 Then I would like to compare the affected genes, for example are the mutations biased towards oncogenes.  
  > **username_3**: 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.  
  > **username_3**: 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: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/t.test.html .  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_3 on 2023-02-28. Score: 3
  
  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](https://hail.is/) (disclaimer: I'm a Hail maintainer):
  
  ```py
  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:
  
  ```py
  from hail.ggplot import *
  fig = ggplot(mt) + geom_point(aes(x=mt.locus, y=mt.variant_qc.AF[0]))
  fig.show()
  
  ```
  
  > **Comments:**
  > **username_0**: 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?  
  
  ---
  
  