dancooke
commented
4 years ago Just to clarify some terminology, as you appear to be using "allele" and "haplotype" inconsistently, when referring to variant calls:
- An allele is one of the values in the
REForALTcolumn of a single VCF record (i.e. a variant site). The site genotype (GT) is composed of alleles. - A haplotype is the sequence determined by the alleles called within some genomic region.
We can only call haplotypes if we can correctly phase alleles at heterozygous variant sites. It's therefore possible to correctly infer alleles and incorrectly infer haplotypes (i.e. phase errors), but not vice versa.
Your interpretation of the phase sets is correct. However, I don't quite understand your comment
This is already a bit confusing for the haplotype reconstruction as only one of the predicted haplotypes (green) has the last 2 SNVs (the lower 2 tracks are simulated reads for the correct alleles)
If I'm reading these pileups correctly (each track corresponds to a strain?), then it appears to me that Octopus has called this region correctly, however, it's certainly possible that the following could have called instead
gene1 62 . A G 3076.53 PASS AC=1;AN=2;DP=1137;MQ=255;NS=1 GT:GQ:DP:MQ:PS:PQ:FT 0|1:3077:1137:255:62:100:PASS
gene1 158 . G A 3076.53 PASS AC=1;AN=2;DP=1509;MQ=255;NS=1 GT:GQ:DP:MQ:PS:PQ:FT 0|1:3077:1509:255:62:100:PASS
gene1 162 . C T 3076.53 PASS AC=1;AN=2;DP=1502;MQ=255;NS=1 GT:GQ:DP:MQ:PS:PQ:FT 1|0:3077:1502:255:62:100:PASS
gene1 419 . T C 3076.53 PASS AC=1;AN=1;DP=1069;MQ=255;NS=1 GT:GQ:DP:MQ:PS:PQ:FT 1|1:15:1069:255:62:100:PASS
gene1 434 . C T 3076.53 PASS AC=1;AN=1;DP=1020;MQ=255;NS=1 GT:GQ:DP:MQ:PS:PQ:FT 1|1:15:1020:255:62:100:PASSThe reason that it wasn't called this way is due to default heuristics in the algorithm that prevent attempting to construct haplotypes where there are clearly no reads able to inform phasing. This heuristic is applied before the downstream (w.r.t the current calling window) sites have been called, so it escapes the algorithm that these sites are actually homozygous and could have been included. You can possibly achieve the calls I proposed by adding the --backtrack-level AGGRESSIVE option, as this will force the algorithm to try to join previously 'called' haplotypes.
I'm struggling to parse your second example. I don't see follow the calls Octopus makes with what I'm seeing on the pileups, plus Octopus is only calling two phase sets here, not 5 as you suggest. Did you post the right example?
Coming back to your question...
I have tested different parameter settings, but I am still not sure which ones would be best suited to improve results, as overall predictions work but are not really accurate. What parameters do I need to set to discriminate between highly similar alleles and enforce strict read phasing, --phasing-level conservative?
What exactly do you mean by "not really accurate", are you seeing false positive or negative variant calls, or do you believe that the phase calls are incorrect? I don't understand what you mean by "to discriminate between highly similar alleles and enforce strict read phasing", what do you mean by "similar" and "strict read phasing"? Octopus phases alleles using genotype posteriors, which are informed by genotype likelihoods (i.e. read information) and prior information, but the prior is only weakly informative for the polyclone model, so is unlikely to outweigh the information from the genotype likelihoods. If you simply want longer phase sets then there are numerous things to try:
- Set
--lagging-level AGGRESSIVE - Set
--backtrack-level AGGRESSIVE - Set
--extension-level AGGRESSIVE - Set
--max-haplotypes 400, or an even larger value (default is 200). - Set
--min-phase-score 5, or lower (default is 10).
But ultimately the phase length will be restricted by the information provided by the reads.
How can I extract the inferred haplotypes from the vcf or would I need to assemble reads in the bam-output with the same "hi" tag?
In what format do you want the inferred haplotypes (i.e. what do you mean by "extract")? If you simply want the sequence then you could try using something like bcftools consensus, but you'll need to manually split the genotypes up into composing haplotypes. It wouldn't make sense to try to "assemble" reads using the hi tag in this way since this is precisely what Octopus has already done. The hi tag tell you which haplotype Octopus thinks the read originated from (if it is not ambiguous).
mherold1
Hello,
thank you for your efforts in providing this tool. I have a question that I could not solve from the extensive documentation and thought maybe you could help me out.
My use case: I have mixtures of bacterial strains and I would like to discriminate the strains using sequencing data for several marker genes that have been amplified. I thought I could use this tool (with the
-C polyclone -P 1option) to identify variants on these marker genes and reconstruct which alleles are likely to be present in the mixture. The examples shown below are from a test sample with a known mix of 3 strains, where strains can share the same allele (example1) or have distinct alleles (example2) of a marker gene. This is the command I used for the examples:My questions: I have tested different parameter settings, but I am still not sure which ones would be best suited to improve results, as overall predictions work but are not really accurate. What parameters do I need to set to discriminate between highly similar alleles and enforce strict read phasing,
--phasing-level conservative? How can I extract the inferred haplotypes from the vcf or would I need to assemble reads in the bam-output with the same "hi" tag?I am also not sure if I understand the vcf output correctly. I found a previous issue that helped me a bit (https://github.com/luntergroup/octopus/issues/32) but this was not specific for ploidy 1. I assume that from the GT field I can reconstruct the genotypes as the order is consistent within the same phase group. So in the following case from the output vcf (example1):
There is phase group 62 (first 3 SNVs:62, 158, 162) and PS 419 with 2 SNVs (419, 434) that would give me 2 distinct alleles regarding these 5 variants: allele1: 1-1-0-1-1 allele2: 0-0-1-1-1 This is already a bit confusing for the haplotype reconstruction as only one of the predicted haplotypes (green) has the last 2 SNVs (the lower 2 tracks are simulated reads for the correct alleles):
A more complex example with a different gene (with 3 distinct alleles):
Here I have several phase groups and tried to "reconstruct" alleles similarly as in example1, in square brackets I noted to which correct phenotype (simulated read tracks) a combination corresponded to: first 5 SNVs in PS 53: 1-0-1-1-1 [1st, green] 1-0-1-0-0 (wrong) 0-1-0-0-0 [2nd and 3rd, blue in the beginning] 2 SNVs in PS 107: 1-0 [1st]
1-0 [1st]
0-1 [2nd and 3rd] 9 SNVs in PS 179: 1-0-1-0-0-0-0-0-0 [1st] 0-1-0-1-1-1-1-1-1 [3rd] 0-1-0-1-1-1-1-2-0 [2nd] 5 SNVs in PS 335: 0-1-0-1-0 [1st] 1-0-0-0-0 [2nd] 1-0-1-0-1 [3rd] last SNVs 2 in PS 440: 1-1 [1st] 0-0 [2nd and 3rd] 0-1 (wrong, but there are actually several reads with this combination)
Why are the same genotypes within a phase group collapsed sometimes (example1, PS:419) but not for others (example2 PS:107)?
output of
octopus --version:Thanks for your help.