Extracting Marker Info ( SNPs used for calling contamination )

[shubhambinf1987]

Hi
We are using VBID2 on some adaptive sampling ONT data sets and wondering how to get the actual markers (SNPs ) used for predicting / detecting contamination ? The log file only provides

[SimplePileup] Total Number Markers: 1713
[SimplePileup] Total Number Bases: 48156

As of now we are using 1000g.phase3.100k.b38 which is significantly slowing the process ( both time and memory )

The main goal is to see if a consistent set of markers used across multiple samples , and then use that to generate custom ref.vcf file following guidelines posted here: https://github.com/Griffan/VerifyBamID, to force VBID2 to only check regions provided in custom reference file instead of all 100K SNPs.

Any other suggestions is greatly appreciated.

Thanks and looking forward.

[Griffan]

Hi @shubhambinf1987,
I see what you are look for. Have you tried option "—OutputPileup", which will output the pileup file for debugging on markers that have effective coverage for the sample. I assume you can take the intersection among these pileup markers across a few testing samples you had. And then subset the reference vcf with these common markers. Please let me know if this works.

[shubhambinf1987]

@Griffan Thanks for the reply.
I also tried your pileup method for ONT reads and while it works at contamination levels of 1,3,5,10 it fails when there is 50 % or 0% ( original sample) contamination , with all PCs = 0 and FREEMIX been 0.5. ( even though logs say success )

Pileup cmd: samtools mpileup
-B
-l 1000g.phase3.100k.b38.vcf.gz.dat.bed

vbid2 cmd:
/verifybamid2
--SVDPrefix 1000g.phase3.100k.b38.vcf.gz.dat
--Reference GCA_000001405.15_GRCh38_no_alt_analysis_set.fasta
--PileupFile
--Output pileup.out.vbid2
--Verbose --NumThread 30 --max-depth 500 --DisableSanityCheck

Is there a way to fix this ( It works with -baminput instead of -pileup but takes forever to run , so pilup input method is optimal in our case)

[Griffan]

Hi @shubhambinf1987, I didn’t mean to use samtools pileup, I meant that VB2 itself has the option "—OutputPileup" to output markers that have effective coverage.

[shubhambinf1987]

Hi @Griffan thanks for the reply.
In order to get --OutputPileup I have to run VBID2 with input BAM which is taking forever to run for ONT data. Hence I was hoping to bypass that by running SAMtools mpileup on the BAM to generate a pileup and pass it using --PileupFile to VBID2 instead of directly passing the BAM file. This is different from the question I initially asked regarding extracting useful markers for detecting contamination.

However --PileupFile method fails at 0% and 50% contamination, where PCs are all 0 and the FREEMIX is 0.5. I was wondering if this is a known bug ?

[Griffan]

How many markers are available under 0, 50%?

[shubhambinf1987]

@Griffan about 1800 SNPs from counting lines in the pileup file which was generated with VBID2 1000g.phase3.100k.b38.vcf.gz.dat
I have attached snapshot of how FREEMIX score failed at different contamination level ( provided by the leftmost col) and the line count of the pileup file provided in the right most col.

[Griffan]

Could you:

1. confirm if those runs are successful, or output "failed to converge"
2. reset different "random seed" and initial values for those failed runs
3. plot allele frequency on these markers, you should expect fixed frequency clusters for 0% and 50% runs, to confirm if there is signal at all

[shubhambinf1987]

Hi Griffan
Sorry for the delayed reply. Here are the responses. We are still seeing VBID2 failing when running with samples that are not contaminated at all ( using the pileup input method ).

1. confirm if those runs are successful, or output "failed to converge" : We ran VBID2 on multiple samples that are clean of any issues ( generating samtools mpileup and then running using --PileupFile option) that doesn't have any contamination and it failed each time, even though it says "NOTICE - Success!" in the logs, FREEMIX:0.5 and PC1/2=0, I have attached a log file here for your kind consideration.

2. reset different "random seed" and initial values for those failed runs: We tried running with different random seed number ( ex: --seed 1, --seed 123 & default ) and still received same outcome.

3. plot allele frequency on these markers, you should expect fixed frequency clusters for 0% and 50% runs, to confirm if there is signal at all : Attached are 2 plots of AF of 2 such samples that has 0 contamination. They are generated by using samtools mpileup on the sample along with 1000g.phase3.100k.b38.vcf.gz.dat.bed, and then running VARSCAN to generate the VCF ( this resulted in approx. 560 variants )

Kindly let us know any advice/guidance you may have in this regard, thanks and appreciate your help in advance.

failed_log.txt

[Griffan]

Hi @shubhambinf1987,
it looks like your mpileup file format was not recognized by VB2, could you please compare the format of the pileup file you provided follows the format described with this line:

VerifyBamID/SimplePileupViewer.cpp

Line 770 in 4b362f5

ss>>pChr>>pPos>>refAllele>>depth>>seq>>qual;

The format my code used may be an outdated single-sample pileup format. Let me know if you find any discrepancy and if it can be easily amended on your side. Thanks!

[shubhambinf1987]

HI @Griffan

I already tried with the correct format of pileup and still same result. Attached is the log as well as the pileup file.

log_github.txt

test.pileup.txt

Thanks.

Hi @shubhambinf1987, it looks like your mpileup file format was not recognized by VB2, could you please compare the format of the pileup file you provided follows the format described with this line:

VerifyBamID/SimplePileupViewer.cpp

Line 770 in 4b362f5

ss>>pChr>>pPos>>refAllele>>depth>>seq>>qual;

The format my code used may be an outdated single-sample pileup format. Let me know if you find any discrepancy and if it can be easily amended on your side. Thanks!

[Griffan]

Hi @shubhambinf1987,
Could you please try two steps:

1. enable sanity check (i.e. remove DisableSanityCheck from your cmdline) to see how many markers we actually have
2. enable the debugging msg at:

   VerifyBamID/ContaminationEstimator.h

   Line 271 in 4b362f5

   // std::cerr <<i<<"th marker\t"<<tmpBase[j]<<"\t"<<tmpQual[j]<<"\t"<<altBase<<"\tlocalAlpha:"<<localAlpha<<"\tgeno1:"<<geno1<<"\tgeno2:"<<geno2