Replies: 5 comments 31 replies
|
Some random thoughts about the future data structure:
|
|
Gene names are present in both the TMB trace files and the Nirvana output, but we might need to ensure their standardization. The Local App Nirvana output includes Ensembl transcript and gene IDs, but the Dragen Nirvana output does not.. I think the molecular biologists prefer RefSeq gene IDs to Ensembl gene IDs anyway (and Nirvana output from both the Local App and Dragen does provide these). |
|
Thanks for adding the info, what I've meant is that gene names are not present in the vcf files - as that file format is not supposed to keep that information. |
|
One could have gene ID information in the INFO VCF column, even potentially utilizing various tags for various gene IDs. Dragen and Local App utilize different tags for gene IDs in the CNV VCFs. But the small variant VCFs produced by the Local App and Dragen don't refer to genes. Utilizing the Nirvana output early would mean that we can assign gene IDs to all variant calls (even those not captured by the TMB file) before we run any own annotation. |
|
I think for different types of gene ids, we should primarly not consider what we need in the end for the molecular biologists but rather something that enables us to clearly identify a gene. The mapping to whatever is necessary for the review by the molecular biologists can then be done from the identifier we agree on - if that is the same or not is only secondary. |
|
I am not currently aware of the need to generate small variant VCF files with gene IDs. That being said, we should have some gene IDs available after the planned Nirvana output parsing, so we could create a function that prints our internal small variant structure into a VCF file if desirable at some point. |
|
As for generating standard vcf files, that should heavily rely on the imported vcf parsing package, in my opinion. However, the package suggested in the guidelines of the develop branch does not seem to be supporting vcf 4.1 or 4.2 (used by localapp and dragen, respectively) - it states working with vcf 4.0, so i can imagine there may be a better package to use, maybe pysam? I will comment on the library use there as well. I cannot find changelog between vcf v4.0 and 4.1, newer versions are addressed here: https://samtools.github.io/hts-specs/VCFv4.3.pdf |
|
Good point, Tina! We should for sure find the package that supports the vcf versions we will work with. Would you be willing to test the different packages with different versions of the vcf, have a look at usability and performance? The candidates I see for now are: |
|
This is a separate issue now. |
The file |
|
The
The difference between the Local App and Dragen formats isn't just in formatting, but also in what kind of details are being provided, which applies especially to the utilized filters. At this point, I don't have enough data to determine which Dragen filters can be trusted. :-/ I don't know what should ideally be done when it comes to 'rescuing' filtered Dragen variants. Unless we can do some benchmarking ahead of our implementation, the safe approach would be to 'rescue too many variants' and let the downstream annotation filter uninteresting variants. I am thus not sure whether we can go for easy conversions here. But I can imagine us introducing new Illumina variant classes (just like TSOPPI introduced the The 'key information' should cover most of our immediate needs (with the exception of variant prioritization and filtering) and we can find it in both the gVCF files and the Nirvana output. The question is whether we should try to integrate info from both sources at the very beginning of the small variant processing, making much of the annotation available early on (TSOPPI utilizes Nirvana output only when constructing the small variant tables in order to look at non-primary transcripts). The gVCFs could in principle provide some additional technical information via the INFO and FORMAT tags, but I don't believe that we would utilize that extra data early on (one could develop own technical filters/indicator based on those values, but I think our filtering will be mostly functional-based to start with). |
|
This might be a recurring theme with other input types as well: Do we now implement the bare minimum that we know we will need, or do we also want to parse out additional information that might (or might not) be useful down the line? |
|
I think parsing all information is always better. The only issue I see is if we want to standardize between LocalApp and Dragen we would need to account for differing or missing information. I think we would need a spec to describe that. |
|
Do we agree on starting with the modest necessary minimum (the key information parsing as described above) for now, in order to implement parsing of what we already know will be required? There are many additional things that might be potentially interesting, or we might never end up using them.. the problem is that determination of what of the extra information (which might be Dragen- or LocalApp-specific) would be truly useful requires research on test data that we don't even have at this point. :-/ |
|
we can agree on the minimum now, for sure and extend when we are ready to. Ok, then I think what you suggest makes perfect sense 👍 |
|
Regarding the dragen |
|
The advantages of using specialized libraries! :-) The parsing part should be easy then, we just need to make sure that we access the desired data correctly for the output of both pipelines. |
|
Thanks to your efforts, we should now have an easier time parsing the JSON and VCF input. :-D How do we approach the variant data structure design? We need to track the key information for a lot of variants (and we might want to parse additional information from the VCF, JSON and/or the TMB files for all the variants after obtaining more experience with the Dragen output, though that is uncertain as of now). Downstream in the post-processing, I expect we will want to track a lot more information for a much smaller variant sub-set (we will then need various additional annotations for implementing our own custom filtering). Do we use the same data structure for everything? We could wait with deciding about the data structure for annotated variant candidates until a later point (and the data structure used for tracking the "key information" could simply be used to initiate a more complex data structure later in the post-processing, if needed). @tinavisnovska @marrip Tina, since the variant parsing is a prerequisite for the predisposition table creation, would you like to pause working on the latter while trying to implement the former? I could perhaps try to focus on the CANCAN/PureCN bit with Martin meanwhile (?). |
|
The whole VCF is relevant for variant recurrence tracking and for rescue of potentially interesting variants that the Illumina pipelines filter out. The rescue operation basically means iterating through the VCF and picking out qualifying variants (with no need to keep the whole VCF in memory). The recurrence tracking was done for convenience reasons with the whole recurrence table in memory, but the input VCF can be simply iterated through as well (1. read in input recurrence table -> 2. update the recurrence table with new info from input VCF files that are being iterated through -> 3. save the updated recurrence table). That being said, when merging data from the VCF, the JSON file and the TMB table, keeping everything in memory would be convenient for speed reasons, and I guess the JSON files are larger anyway and we intend to load them into memory as well? (The JSON file should include records for all the VCF variants and provide more information for them than the VCF.) |
Do you want to spend more time on the testing or do you think that using a vcf object to store all the merged info is a viable option here? |
|
As we might want to write a vcf again it would be good to keep the vcf object? I am also not super concerned with size here unless we discover they vcfs are bigger than 1 GB. |
|
Let's conclude on how to proceed with the merging, as it's a prerequisite for further small variant post-processing development. :-) I wonder about the following if we do go for the VCF extension:
Regarding point 2), I assume we'd have to write some custom functions that operate on the VCF objects (?). One could make them general, for different types of INFO fields (there are Integer/Float/Flag/Character/String fields, with one or with multiple values). Sometimes we'd like get all variants with a given float INFO field value under a specific threshold (e.g., when looking at variant frequency in a population), other times all variants with a specific flag INFO field value present (e.g., all variants with mutation hotspot status), etc. These are the INFO field items we would need from the other variant files:
In addition, some of the gVCF information could be re-formatted and/or post-processed into new INFO fields:
|
|
I think what you suggest sounds reasonable. I would be willing to give it a try and implement those methods but I would need some more detailed instructions I think. For example, do we need to add a header for each new INFO field we add. I will start on a new branch with my work and we can keep a tight feedback loop? |
Uh oh!
There was an error while loading. Please reload this page.
Uh oh!
There was an error while loading. Please reload this page.
Inputs considered so far:
Pair_IDandSample_IDvalues were used when generating the output referenced in the text belowPair_IDvalue choice) can impact the analysis output structurePair_IDandSample_IDvalues were used when generating the output referenced in the text below5. Small variants
5.1. Input files
5.1.1. Genomic VCFs - records of small variants (or the lack thereof) at all targeted positions
Dragen TSO500 2.6.2:
panel_of_normals: "Seen in at least one sample in the panel of normals vcf"multiallelic: "Site filtered because more than two alt alleles pass tumor LOD"excluded_regions: "Site overlaps with vc excluded regions bed" (Dragen "Block List" - HLA A/B/C, KMT2B/C/D, chrY, any position with VAF > 0.01 in at least 6 of 60 baseline samples)Local App 2.2.0.12:
Blacklist(see the "Block List" and theexcluded_regionsFILTER in case of Dragen)5.1.2. TMB trace files - small variants that made it through all Illumina filtering steps
Dragen TSO500 2.6.2:
Status(Germline_DB/Germline_Proxi/Somatic, corresponding to TSOPPI'sGL_DB/GL_P/SOM)ClinVarIDs,ConsequenceLocal App 2.2.0.12:
GermlineFilterDatabaseandGermlineFilterProxiGermlineFilterDatabase== True) -> the class should beGermline_DB(GL_DB)GermlineFilterProxi== True) -> the class should beGermline_Proxi(GL_P)Somatic(SOM)5.1.3. Nirvana annotation output - more extensive small variant annotation
Dragen TSO500 2.6.2:
Local App 2.2.0.12:
5.2. Small variant data usage in tsoppy modules
5.2.1. Key variant information that needs to be extracted:
5.2.2. Modules that would utilize small variant data
Modules that can do with the key information listed above:
Modules that will ideally require additional data/pre-processing:
5.2.3. Suggestions for implementation
GL_DB,GL_PandSOMclasses are specified in the TMB trace filesBL, can be created based on the tolerated corresponding FILTERsWhat internal format do we use for storing the small variant information inside tsoppy? The annotations might require tracking multiple instances of certain values (e.g., anything tied to individual transcripts), we would like to be able to add flags and filters of our own, and to efficiently select variant sub-sets based on these. Lastly, if would be nice if we could generate a standard VCF at any point from the internal format.
Suggested sources of the "key variant information":
1. genomic position (chromosome + position) and genomic change (reference allele + alternative allele):
:as the separator of the contributing elements (e.g.,chr7:140453136:A:T)Primary source: the
CHROM,POS,REFandALTcolumns of the gVCF filesALTcolumn can include multiple items and while neither the Local App nor the Dragen report multiple variant calls within a single gVCF line/entry, the Dragen output includes an artificial<NON_REF>'variant' on each line, potentially in addition to an actual variant (the<NON_REF>item should be removed from the variant ID string)ALTallele is equal to simply<NON_REF>(in case of Dragen) or.(in case of the Local App), no variant is reported at given position (there won't be any matching JSON or TMB trace file record)<NON_REF>allele placeholder - we could consider using this input insteadFILTERfield value not beingPASS) won't have a matching TMB trace file record, but they will have a JSON record and might be of interestChromosome,Position,RefCallandAltCallcolumns should be used to construct the IDchromosome,position,refAlleleandaltAllelesshould be used to construct the ID (thealtAllelesbeing a list with a single value)positionsarraypositionsarray entry per variant reported in the gVCF (the gVCF has additional entries for targeted positions with no variants)variantsitem (located on the same level aschromosome,position,refAlleleandaltAlleles) can be confusing, as it might for some variant types show different values for its ownchromosome,begin,refAlleleandaltAllelesub-items due to utilizing a different notation.. I recommend ignoring the confusing notation undervariantsSuggested INFO field header for the variant ID:
##INFO=<ID=variant_ID,Number=1,Type=String,Description="Variant identifier, composed of CHROM, POS, REF and ALT fields">
(This assumes that all merging output records include exactly one variant. Ideally, all
<NON_REF>alleles and non-variant records are removed/ignored - this could be made easier by working with Dragen VCF files instead of gVCF files - as mentioned above, but the Local App output doesn't seem to contain a corresponding alternative input option.)2. variant type
Possible sources:
2.1. The
variantTypeitem in the Nirvana JSON file (positions->variants->variantType); available for all called variantsMNV,SNV,insertionanddeletionMNV,SNV,insertion,deletionandindel2.2. The
VariantTypecolumn in the TMB trace file; only available for variants passing all filtersdeletion,insertion,MNV,SNVDELETION,INSERTION,RECOMPOSED,SNVThe meaning of the
RECOMPOSEDvariant type in the TMB trace files is unclear. I suggest that we utilize the Nirvana output, which seems more suitable. Dragen'sindeltype seems to be used for genomic changes that include both an insertion and a deletion.Suggested INFO field header for the variant type:
##INFO=<ID=variant_type,Number=1,Type=String,Description="Type of genomic change: "MNV", "SNV", "insertion", "deletion" or "indel"">
3. variant class
A distinction should be made between variants that passed all the filters and made it into the TMB trace file and the remaining (failed) variants.
3.1. The gVCFs provide information about failing pipeline filters for all called variants in their
FILTERcolumn: Variants that did not fail any filter are assigned valuePASS, the remaining variants are assigned a list of applied/failed filters (all of which should be documented in the "##FILTER" header items)The Nirvana JSON output should provide the same filtering information as the gVCF file via its
positions->filterslist item (thefilterslist is on the same level as thechromosome,position,refAlleleandaltAllelesitems.This information doesn't need to be duplicated in the INFO field, as it would be straight copy/paste.
3.2 Variants passing all filters should have an entry in the TMB trace file, which provides further classification.
Status, with possible values ofGermline_DB,Germline_ProxiandSomaticGermlineFilterDatabasecolumn has valueTrue, then the class should beGermline_DBGermlineFilterProxicolumn has valueTrue, then the class should beGermline_ProxiSomaticSuggested INFO field header for the variant class:
##INFO=<ID=Illumina_variant_class,Number=1,Type=String,Description="Variant classification provided by initial Illumina pipeline analysis: "Germline_DB", "Germline_Proxi", "Somatic", "Blacklist", "VCF_filtered"">
Variants that fail one of more VCF filters will not have an entry in the TMB trace file (and thus will not have a
Germline_DB,Germline_ProxiorSomaticclass). TSOPPI previously utilized a class "BL" ("Blacklist"/"Blocklist") for variants that failed only the Local AppBlacklistfilter (these were "rescued" and reported on). In case of Dragen output, this would correspond to a variant failing only theexcluded_regionsfilter. I suggest we continue utilizing the "Blacklist" class until/unless we find a reason not to do so. All variants not rescued from filtering in any way could be assigned classVCF_filtered.4. read depth, allelic depth and allelic fraction
There are two options for the primary source:
4.1. The gVCF files, namely their genotype fields for a given sample (the last column in the data section of the file). The preceding
FORMATfield will define the relative position of the relevant keywords.AD,DPandVFfor allelic depth, read depth and variant allelic fraction, respectively. For each called variant, there is a single read depth value, one allelic depth value per allele (reference, variant) and a single variant allelic fraction value.AD,DPandAFfor allelic depth, read depth and allelic fraction, respectively. For each called variant, there is a single read depth value, one allelic depth value per allele (reference, variant, <NON_REF>) and one allelic fraction value per non-reference allele (variant, <NON_REF>).Please note that the number and the order of keywords can vary from line to line (so the
FORMATfield should in principle be parsed for each variant separately).4.2. The Nirvana JSON output. Under positions -> samples, the following items define the allelic depth, read depth and variant allelic fraction values:
alleleDepths,totalDepthandvariantFrequenciesPlease note that the <NON_REF> allele is not taken into account in the Nirvana output of Dragen (unlike in its gVCF format).
No INFO field if we parse the data from the standard FORMAT/Genotype fields instead. Please note the "VF" vs. "AF" FORMAT/Genotype field naming inconsistency and the variable <NON_REF> allele usage.
5. overlapping genes
Primary source: The Nirvana JSON output, namely items
positions->variants->transcripts->hgnc.Suggested INFO field header for overlapping genes:
##INFO=<ID=overlapping_genes,Number=.,Type=String,Description="List of genes overlapping the variant site, as reported by Nirvana in form of HGNC gene symbols">
All reactions