A weird oncoplot

[beginner984]

Hi

I have whole exome sequencing of eight samples. I have annotated the vcf files with annovar and by maftools I see a weird number of mutations

> my_maf
An object of class  MAF 
ID summary     Mean Median
1:             NCBI_Build      NA       NA     NA
2:                 Center      NA       NA     NA
3:                Samples       8       NA     NA
4:                 nGenes     572       NA     NA
5:        Frame_Shift_Del    1680  210.000  221.0
6:        Frame_Shift_Ins     839  104.875  102.5
7:           In_Frame_Del     124   15.500   16.0
8:           In_Frame_Ins      59    7.375    5.0
9:      Missense_Mutation    6839  854.875  852.5
10:      Nonsense_Mutation     483   60.375   60.0
11:       Nonstop_Mutation      12    1.500    0.5
12:            Splice_Site     151   18.875   20.0
13: Translation_Start_Site       3    0.375    0.0
14:                  total   10190 1273.750 1274.5

You see I have a median of 1274.5 mutations which is really weird to me in compared to a maf object from maftools

> maftool_maf
An object of class  MAF 
ID          summary  Mean Median
1:        NCBI_Build               37    NA     NA
2:            Center genome.wustl.edu    NA     NA
3:           Samples              193    NA     NA
4:            nGenes             1241    NA     NA
5:   Frame_Shift_Del               52 0.269      0
6:   Frame_Shift_Ins               91 0.472      0
7:      In_Frame_Del               10 0.052      0
8:      In_Frame_Ins               42 0.218      0
9: Missense_Mutation             1342 6.953      7
10: Nonsense_Mutation              103 0.534      0
11:       Splice_Site               92 0.477      0
12:             total             1732 8.974      9

I am seeing a weird oncoplot like this

As you are seeing for every gene all samples just have a multi hit event

> oncoplot(maf = anno_maf, top = 10)
Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
> 

Have you ever had such an observation?

[ShixiangWang]

@beginner984 Hi, why your WES got so many variants? I think you should double-check your mutation caller. If possible, use two different softwares to confirm the mutation list.

I think this is not an issue directly related to ANNOVAR or maftools as they are just used for annotation and downstream visualization/analysis.

[beginner984]

@ShixiangWang thank you for answering me

So you also agree that I have a lot of mutations? Do you think if filtering can help? I have tumour only samples without matched normal, so might be most of mutations are gremlin because there is no normal sample to filter them

What you think?

[ShixiangWang]

@beginner984 I am not familiar with the mutation calling processing, but using a pool of normal from other data or annotated SNPs from database like dbSNP should help.

It's rare to see more than one somatic mutation occur in the same gene. So I think filtering would work.

[CC-DrDU]

Hi. I've got the same problem about the oncoplot, and I also only have WES data of tumor cell lines. Do you think it is ok to perform Mutect2 for somatic variant calling using matched tumor cell lines?

[ShixiangWang]

Why match tumor cell lines? It cannot be a reference I think.

GATK has tumor only mode in mutect2, it may help.

Please search the doc at:

https://gatk.broadinstitute.org/hc/en-us

[CC-DrDU]

It seems that the tumor-only mode in Mutect2 encourages to include "--panel-of-normals", but I dont have any "normal" data for creating it, so I was confused these days and finally used HaplotypeCaller to call variant for each tumor cell lines. I think that's the reason why I got such a weird oncoplot.

Do you think it is ok to use the pon frm 1000g (https://console.cloud.google.com/storage/browser/gatk-best-practices/somatic-hg38;tab=objects?prefix=&forceOnObjectsSortingFiltering=false), or should I just ignore the "--panel-of-normals"?

Thanks.

[CC-DrDU]

Cuz from the way I see it, neither of them are the best way. T-T

[ShixiangWang]

@CC-DrDU You can obtain panel-of-normal from public available data (hub) like TCGA/biobank (of course 1000 genomes). And also you have to filter out potential SNPs as they are not somatic variants.

That's all I know. I have no practice on WES mutation calling, so I haven't further deeper comments here.

[CC-DrDU]

@ShixiangWang Thank you for your help! I'll try. In fact I just found that David Benjamin from GATK4 had a response for a similar question. When only tumor samples are available, "Run Mutect2 and FilterMutectCalls with the best practices resources for germline resource and panel of normals. The Mutect2 pipeline is self-contained."

[ShixiangWang]

@CC-DrDU Thanks for letting me know.

[beginner984]

I need help please

That is too confusing

I re-called variants by Mutect2 with tumor only mode only (because I have whole exome sequencing of only tumor samples) but again the same weird oncoplot and a huge number of variants

This is my full code to use Mutect2 and then annovar and maftools

./gatk Mutect2 -R reference.fa -I sample.bam --germline-resource somatic-hg38_af-only-gnomad.hg38.vcf.gz --panel-of-normals 1000g_pon.hg38.vcf.gz --f1r2-tar-gz f1r2.tar.gz -O unfiltered.vcf

./gatk LearnReadOrientationModel -I f1r2.tar.gz -O read-orientation-model.tar.gz

./gatk GetPileupSummaries -I sample.bam -V somatic-hg38_small_exac_common_3.hg38.vcf.gz -L somatic-hg38_small_exac_common_3.hg38.vcf.gz -O getpileupsummaries.table

./gatk CalculateContamination -I getpileupsummaries.table -tumor-segmentation segments.table -O calculatecontamination.tabl

./gatk FilterMutectCalls -R reference.fa -V unfiltered.vcf --tumor-segmentation segments.table --contamination-table calculatecontamination.table --ob-priors read-orientation-model.tar.gz -O filtered.vcf


bcftools view -f PASS filtered.vcf > output.vcf

perl table_annovar.pl /data2/RNASeq/Angel/gatk/output.vcf humandb/ -buildver hg38 -out myanno -remove -protocol refGene,cytoBand,dbnsfp30a -operation g,r,f -nastring . -vcfinput

 var.annovar="myanno.hg38_multianno.txt"

> var.annovar.maf = annovarToMaf(annovar = var.annovar,refBuild = 'hg38',table = 'refGene')
-Reading annovar files
--Tumor_Sample_Barcode column not found. Creating sample IDs from filenames
-Processing Exonic variants
--Adding Variant_Classification
--Parsing aa-change
-Processing Non-exonic variants
--Adding Variant_Classification
-Adding Variant_Type
Finished in 1.459s elapsed (6.330s cpu) 
> laml = read.maf(var.annovar.maf)
-Validating
--Non MAF specific values in Variant_Classification column:
  Unknown
--Non MAF specific values in Variant_Type column:
  MNP
-Silent variants: 72743 
-Summarizing
-Processing clinical data
--Missing clinical data
-Finished in 0.284s elapsed (1.809s cpu) 

You are seeing the same weird oncolplot

Do you have any idea what should I do?

[CC-DrDU]

You only have one sample? That's not a usual thing......
It is very important to filter your vcf file(s) I guess. I am also a beginnner and I dont know how to filter my files after ANNOVAR annotation, so I just use SnpEff & SnpSift to filter my vcf files from Mutect2 pipelines, and then use ANNOVAR to annotate (although I know that SnpEff is also an annotation software LOL). After that my oncoplot seems perfect. Hope that this could help you.

[beginner984]

Thank you I have eight samples, I just showed one of them here

If you go through my code, you will see I have filtered my variants :(

[beginner984]

Sorry @CC-DrDU Please can you share your code by which you have called the mutations and filtered them using Mutect2 and SnpEff & SnpSift?

I think this extremely helps me to follow you and obtain a reasonable number of mutations

Thank you so much in advance

[CC-DrDU]

For each vcf generated from Mutect2, I do:
java -Xmx20g -jar $snpEff -v GRCh38.99 ${id}.mutect.filtered.vcf.gz > ${id}.mutect.filtered.ann.vcf.gz
cat ${id}.mutect.filtered.ann.vcf.gz | java -jar $SnpSift filter "( na FILTER ) | (FILTER = 'PASS')" > ${id}.mutect.filtered.ann.snpsift.vcf

And then I use ANNOVAR to annotate the "${id}.mutect.filtered.ann.snpsift.vcf". You can also read the document about SnpSift. It may help.

[CC-DrDU]

I dont think I'm doing in the smart way, but I am still learning about the filtering after annotation......

[beginner984]

Thank you for your help

I have done this but again the same huge number of mutations

Please look at my code

java -jar snpEff.jar -v hg38 /data2/RNASeq/Angel/gatk/mutect.filtered.vcf.gz > mutect.filtered.ann.vcf.gz

then

cat mutect.filtered.ann.vcf.gz | java -jar SnpSift.jar filter "( na FILTER ) | (FILTER = 'PASS')" > mutect.filtered.ann.snpsift.vcf

Then annovar

perl table_annovar.pl /data2/RNASeq/Angel/gatk/mutect.filtered.ann.snpsift.vcf humandb/ -buildver hg38 -out myanno -remove -protocol refGene,cytoBand,dbnsfp30a -operation g,r,f -nastring . -vcfinput

I don't know where I am doing wrong

Please can you share the code by which you obtained mutect.filtered.vcf.gz file from a bam file ??

May be I am doing something wrong with Mutect2

Thanks a lot for any help

[beginner984]

Sorry @CC-DrDU

Please can you share the code by which you called mutations with Mutect2 ??

As snpEff and annovar part I followed you so I think if I follow you on calling mutations with Mutect2, I might be solve the problem

Thanks for any help

[CC-DrDU]

I am sorry I forgot to write a reply...

I think you could change the pon you use. As you mentioned, you're using the 1000g pon from the GATK bundle, right? I've used it before and I also found a lot of multi-hits. So in the end I used my wild-type samples to create a pon (although they are also tumor cell lines and are not "normal" samples).

I think it should help. However I am not sure whether it's appropriate to create pon with samples that are nor "normal".

[beginner984]

Thank you so much for replying me

Can you share your code with me please?

The lines of codes you used for a bam file to obtain a filtered vcf file from Mutect2

Thanks for any help

[beginner984]

Thank you so much

I don't have two groups, rather I have eight tumor samples

Sorry, for this part normals=(./*.vcf.gz) which .vcf files I should provide?

[CC-DrDU]

If you dont have samples for reference (which means that all of your samples are in experimental group), the only thing you can do is to call mutations using your original code, which would always cause the oncoplot problem, because tumor cell lines have a lot of differences comparing to the reference genome and pon from normal samples.

[beginner984]

Sorry, for this part normals=(./*.vcf.gz) which .vcf files I should provide?

[CC-DrDU]

If you have reference samples, the " normals=(./*.vcf.gz)" would be the reference samples. But since you dont have reference samples, I dont think it is appropriate to do it in my way mentioned above. But you do can continue without oncoplot, the .txt file from ANNOVAR is enough for you to continue.

[CC-DrDU]

As you dont have two groups, there's nothing I can help further... T T

[arvinhm]

The easiest way to filter your VCF file (even without normal samples) is to use OpenCravat. You can choose multiple platforms to annotate your VCF file correctly. Then export the VCF file and use ANNOVAR to annotate (because the best practice for generation MAF file to use in maftools is using ANNOVAR rather than VCF2MAF, my experience). No coding requires in Opne Cravat, amazing. In addition, you can create your own MAF file by adding lots of feature and analyze your data, for example, you can only focus on the mutations resulted in protein dysfunction (hyper/hypo active).