// createGeneticAnalysis provides functions to create a genetic analysis // These are performed on one or more genome files. // They produce 3 kinds of results: Monogenic Diseases, Polygenic Diseases and Traits // They can be performed on a Person or a Couple // Couple analyses provide an analysis of the prospective offspring of the couple package createGeneticAnalysis // TODO: Some of the probabilities produced by this package are wrong // In this package, we are assuming that genetic recombination (the formation of the genetic sequences for the sperm/eggs) happens randomly for each allele locus // In reality, the recombination break points occur less often, and larger portions of each chromosome remain intact. // This effects the estimates and probabilities for all of the generated analyses // In particular, the probability of passing a defective gene does not increase as much as this package currently // estimates that it does, in the case of multiple defects existing in the same gene. // Also, based on my research, I believe that recombination break points are less likely to occur within genes, meaning they are more likely to occur at the gene boundaries (codons) // We need to remedy this problem and fix this package // Research gene linkage and recombination to understand more. // // Fixing this problem may require us to only allow will-pass-a-variant probabilities to be calculated from phased loci. // The phase of a loci becomes much more important in determining the will-pass-a-variant probability // Having multiple variants within a gene might not increase the probability of passing a variant, assuming all of those variants were on the same chromosome // Thus, we need phased loci to determine an accurate will-pass-a-variant probability // We will still be able to determine will-pass-a-variant probabilities for users who only have 1 variant on 1 base, // regardless of if their loci phase is known or not. That probability is 50%. // TODO: We want to eventually use neural nets for both trait and polygenic disease analysis // These will be trained on a set of genomes and will output a probability analysis for each trait/disease // This is only possible once we get access to the necessary training data // // This is how offspring trait prediction could work with the neural net model: // Both users will share all relevant SNPS base pairs that determine the trait on their profile. // Each location has 4 possible outcomes, so for 1000 SNPs, there are 4^1000 possible offspring outcomes for a given couple. (this // is actually too high because recombination break points do not occur at each locus, see genetic linkage) // This is too many options for us to check all of them. // Seekia will create 100 offspring that would be produced from both users, and run each offspring through the neural net. // Each offspring would be different. The allele from each parent for each SNP would be randomly chosen. // The user can choose how many prospective offspring to create in the settings. // More offspring will take longer, but will yield a more accurate trait probability. // Seekia will show the the average trait result and a chart showing the trait results for all created offspring. // TODO: We should not add any passing-a-variant probabilities that are dependent on an unknown phase to user profiles // For example, if the phase for any mutated bases are unknown, and this unknown phase effects the probability of // the user passing a variant, that probability should not be shared on the user's profile. // Seekia currently shares these probabilities on user profiles. // Only probabilities which are fully known should be shared. Otherwise, the user is sharing // information that is not fully accurate but is instead speculative. // // We must add this warning into the GUI for a user and couple genetic analyses too. // Basically, we need to make it clear that some probabilities are based completely on random genetic chance, and others // are based on our limited understanding of the user's genome (due to unphased data). // TODO: Add the ability to weight different genome files based on their reliability. // Some files are much more accurate because they record each location many times. import "seekia/resources/geneticReferences/locusMetadata" import "seekia/resources/geneticReferences/monogenicDiseases" import "seekia/resources/geneticReferences/polygenicDiseases" import "seekia/resources/geneticReferences/traits" import "seekia/internal/genetics/locusValue" import "seekia/internal/genetics/prepareRawGenomes" import "seekia/internal/helpers" import "errors" import "math" import "encoding/json" import "strings" import "slices" func verifyBasePair(inputBasePair string)bool{ base1, base2, delimiterFound := strings.Cut(inputBasePair, ";") if (delimiterFound == false){ return false } // I = Insertion // D = Deletion validBasesList := []string{"C", "A", "T", "G", "I", "D"} baseIsValid := slices.Contains(validBasesList, base1) if (baseIsValid == false){ return false } baseIsValid = slices.Contains(validBasesList, base2) if (baseIsValid == false){ return false } return true } //Inputs: // -map[string]string: Genome Identifier -> Genome Raw data string //Outputs: // -bool: Process completed (it was not stopped manually before completion) // -string: New Genetic analysis string // -error func CreatePersonGeneticAnalysis(genomesList []prepareRawGenomes.RawGenomeWithMetadata, updatePercentageCompleteFunction func(int)error, checkIfProcessIsStopped func()bool)(bool, string, error){ prepareRawGenomesUpdatePercentageCompleteFunction := func(newPercentage int)error{ newPercentageCompletion, err := helpers.ScaleNumberProportionally(true, newPercentage, 0, 100, 0, 50) if (err != nil){ return err } err = updatePercentageCompleteFunction(newPercentageCompletion) if (err != nil) { return err } return nil } genomesWithMetadataList, allRawGenomeIdentifiersList, multipleGenomesExist, onlyExcludeConflictsGenomeIdentifier, onlyIncludeSharedGenomeIdentifier, err := prepareRawGenomes.GetGenomesWithMetadataListFromRawGenomesList(genomesList, prepareRawGenomesUpdatePercentageCompleteFunction) if (err != nil) { return false, "", err } combinedGenomesExistString := helpers.ConvertBoolToYesOrNoString(multipleGenomesExist) metadataItem := map[string]string{ "ItemType": "Metadata", "AnalysisVersion": "1", "AnalysisType": "Person", "CombinedGenomesExist": combinedGenomesExistString, } if (multipleGenomesExist == true){ metadataItem["OnlyExcludeConflictsGenomeIdentifier"] = onlyExcludeConflictsGenomeIdentifier metadataItem["OnlyIncludeSharedGenomeIdentifier"] = onlyIncludeSharedGenomeIdentifier allRawGenomeIdentifiersListString := strings.Join(allRawGenomeIdentifiersList, "+") metadataItem["AllRawGenomeIdentifiersList"] = allRawGenomeIdentifiersListString } else { genomeIdentifier := allRawGenomeIdentifiersList[0] metadataItem["GenomeIdentifier"] = genomeIdentifier } processIsStopped := checkIfProcessIsStopped() if (processIsStopped == true){ return false, "", nil } geneticAnalysisMapList := []map[string]string{metadataItem} monogenicDiseasesList, err := monogenicDiseases.GetMonogenicDiseaseObjectsList() if (err != nil) { return false, "", err } for _, monogenicDiseaseObject := range monogenicDiseasesList{ diseaseName := monogenicDiseaseObject.DiseaseName diseaseAnalysisMap, variantsMapList, err := getMonogenicDiseaseAnalysis(genomesWithMetadataList, monogenicDiseaseObject) if (err != nil) { return false, "", err } diseaseAnalysisMap["ItemType"] = "MonogenicDisease" diseaseAnalysisMap["DiseaseName"] = diseaseName geneticAnalysisMapList = append(geneticAnalysisMapList, diseaseAnalysisMap) geneticAnalysisMapList = append(geneticAnalysisMapList, variantsMapList...) } polygenicDiseaseObjectsList, err := polygenicDiseases.GetPolygenicDiseaseObjectsList() if (err != nil) { return false, "", err } for _, diseaseObject := range polygenicDiseaseObjectsList{ diseaseAnalysisMap, diseaseLociMapList, err := getPolygenicDiseaseAnalysis(genomesWithMetadataList, diseaseObject) if (err != nil) { return false, "", err } diseaseName := diseaseObject.DiseaseName diseaseAnalysisMap["ItemType"] = "PolygenicDisease" diseaseAnalysisMap["DiseaseName"] = diseaseName geneticAnalysisMapList = append(geneticAnalysisMapList, diseaseAnalysisMap) geneticAnalysisMapList = append(geneticAnalysisMapList, diseaseLociMapList...) } traitObjectsList, err := traits.GetTraitObjectsList() if (err != nil) { return false, "", err } for _, traitObject := range traitObjectsList{ traitAnalysisMap, traitRulesMapList, err := getTraitAnalysis(genomesWithMetadataList, traitObject) if (err != nil) { return false, "", err } traitName := traitObject.TraitName traitAnalysisMap["ItemType"] = "Trait" traitAnalysisMap["TraitName"] = traitName geneticAnalysisMapList = append(geneticAnalysisMapList, traitAnalysisMap) geneticAnalysisMapList = append(geneticAnalysisMapList, traitRulesMapList...) } analysisBytes, err := json.MarshalIndent(geneticAnalysisMapList, "", "\t") if (err != nil) { return false, "", err } analysisString := string(analysisBytes) return true, analysisString, nil } //Outputs: // -bool: Process completed (was not manually stopped mid-way) // -string: Couple genetic analysis string // -error func CreateCoupleGeneticAnalysis(personAGenomesList []prepareRawGenomes.RawGenomeWithMetadata, personBGenomesList []prepareRawGenomes.RawGenomeWithMetadata, updatePercentageCompleteFunction func(int)error, checkIfProcessIsStopped func()bool)(bool, string, error){ personAPrepareRawGenomesUpdatePercentageCompleteFunction := func(newPercentage int)error{ newPercentageCompletion, err := helpers.ScaleNumberProportionally(true, newPercentage, 0, 100, 0, 25) if (err != nil){ return err } err = updatePercentageCompleteFunction(newPercentageCompletion) if (err != nil) { return err } return nil } personAGenomesWithMetadataList, allPersonARawGenomeIdentifiersList, personAHasMultipleGenomes, personAOnlyExcludeConflictsGenomeIdentifier, personAOnlyIncludeSharedGenomeIdentifier, err := prepareRawGenomes.GetGenomesWithMetadataListFromRawGenomesList(personAGenomesList, personAPrepareRawGenomesUpdatePercentageCompleteFunction) if (err != nil) { return false, "", err } processIsStopped := checkIfProcessIsStopped() if (processIsStopped == true){ return false, "", nil } personBPrepareRawGenomesUpdatePercentageCompleteFunction := func(newPercentage int)error{ newPercentageCompletion, err := helpers.ScaleNumberProportionally(true, newPercentage, 0, 100, 25, 50) if (err != nil){ return err } err = updatePercentageCompleteFunction(newPercentageCompletion) if (err != nil) { return err } return nil } personBGenomesWithMetadataList, allPersonBRawGenomeIdentifiersList, personBHasMultipleGenomes, personBOnlyExcludeConflictsGenomeIdentifier, personBOnlyIncludeSharedGenomeIdentifier, err := prepareRawGenomes.GetGenomesWithMetadataListFromRawGenomesList(personBGenomesList, personBPrepareRawGenomesUpdatePercentageCompleteFunction) if (err != nil) { return false, "", err } processIsStopped = checkIfProcessIsStopped() if (processIsStopped == true){ return false, "", nil } // The analysis will analyze either 1 or 2 genome pairs // The gui will display the results from each pair //Outputs: // -string: Pair 1 Person A Genome Identifier // -string: Pair 1 Person B Genome Identifier // -bool: Second pair exists // -string: Pair 2 Person A Genome Identifier // -string: Pair 2 Person B Genome Identifier // -error getGenomePairsToAnalyze := func()(string, string, bool, string, string, error){ if (personAHasMultipleGenomes == true && personBHasMultipleGenomes == true){ return personAOnlyExcludeConflictsGenomeIdentifier, personBOnlyExcludeConflictsGenomeIdentifier, true, personAOnlyIncludeSharedGenomeIdentifier, personBOnlyIncludeSharedGenomeIdentifier, nil } if (personAHasMultipleGenomes == true && personBHasMultipleGenomes == false){ personBGenomeIdentifier := allPersonBRawGenomeIdentifiersList[0] return personAOnlyExcludeConflictsGenomeIdentifier, personBGenomeIdentifier, true, personAOnlyIncludeSharedGenomeIdentifier, personBGenomeIdentifier, nil } if (personAHasMultipleGenomes == false && personBHasMultipleGenomes == true){ personAGenomeIdentifier := allPersonARawGenomeIdentifiersList[0] return personAGenomeIdentifier, personBOnlyExcludeConflictsGenomeIdentifier, true, personAGenomeIdentifier, personBOnlyIncludeSharedGenomeIdentifier, nil } //personAHasMultipleGenomes == false && personBHasMultipleGenomes == false personAGenomeIdentifier := allPersonARawGenomeIdentifiersList[0] personBGenomeIdentifier := allPersonBRawGenomeIdentifiersList[0] return personAGenomeIdentifier, personBGenomeIdentifier, false, "", "", nil } pair1PersonAGenomeIdentifier, pair1PersonBGenomeIdentifier, genomePair2Exists, pair2PersonAGenomeIdentifier, pair2PersonBGenomeIdentifier, err := getGenomePairsToAnalyze() if (err != nil){ return false, "", err } personAHasMultipleGenomesString := helpers.ConvertBoolToYesOrNoString(personAHasMultipleGenomes) personBHasMultipleGenomesString := helpers.ConvertBoolToYesOrNoString(personBHasMultipleGenomes) secondPairExistsString := helpers.ConvertBoolToYesOrNoString(genomePair2Exists) metadataItem := map[string]string{ "ItemType": "Metadata", "AnalysisVersion": "1", "AnalysisType": "Couple", "PersonAHasMultipleGenomes": personAHasMultipleGenomesString, "PersonBHasMultipleGenomes": personBHasMultipleGenomesString, "Pair1PersonAGenomeIdentifier": pair1PersonAGenomeIdentifier, "Pair1PersonBGenomeIdentifier": pair1PersonBGenomeIdentifier, "SecondPairExists": secondPairExistsString, } if (genomePair2Exists == true){ // At least 1 of the people have multiple genomes metadataItem["Pair2PersonAGenomeIdentifier"] = pair2PersonAGenomeIdentifier metadataItem["Pair2PersonBGenomeIdentifier"] = pair2PersonBGenomeIdentifier if (personAHasMultipleGenomes == true){ metadataItem["PersonAOnlyExcludeConflictsGenomeIdentifier"] = personAOnlyExcludeConflictsGenomeIdentifier metadataItem["PersonAOnlyIncludeSharedGenomeIdentifier"] = personAOnlyIncludeSharedGenomeIdentifier } if (personBHasMultipleGenomes == true){ metadataItem["PersonBOnlyExcludeConflictsGenomeIdentifier"] = personBOnlyExcludeConflictsGenomeIdentifier metadataItem["PersonBOnlyIncludeSharedGenomeIdentifier"] = personBOnlyIncludeSharedGenomeIdentifier } } newAnalysisMapList := []map[string]string{metadataItem} //First we add monogenic disease probabilities monogenicDiseasesList, err := monogenicDiseases.GetMonogenicDiseaseObjectsList() if (err != nil) { return false, "", err } for _, diseaseObject := range monogenicDiseasesList{ diseaseName := diseaseObject.DiseaseName variantsList := diseaseObject.VariantsList diseaseIsDominantOrRecessive := diseaseObject.DominantOrRecessive offspringDiseaseMap := map[string]string{ "ItemType": "MonogenicDisease", "DiseaseName": diseaseName, } personADiseaseAnalysisMap, personAVariantsMapList, err := getMonogenicDiseaseAnalysis(personAGenomesWithMetadataList, diseaseObject) if (err != nil) { return false, "", err } personBDiseaseAnalysisMap, personBVariantsMapList, err := getMonogenicDiseaseAnalysis(personBGenomesWithMetadataList, diseaseObject) if (err != nil) { return false, "", err } // This will calculate the probability of monogenic disease for the offspring from the two specified genomes // It then adds the pair entry to the offspringDiseaseMap addPairProbabilitiesToDiseaseMap := func(personAGenomeIdentifier string, personBGenomeIdentifier string)error{ //Outputs: // -bool: Probability is known // -int: Probability of passing a disease variant (value between 0 and 100) // -int: Number of variants tested // -error getPersonWillPassDiseaseVariantProbability := func(personDiseaseAnalysisMap map[string]string, genomeIdentifier string)(bool, int, int, error){ probabilityMapKey := genomeIdentifier + "_ProbabilityOfPassingADiseaseVariant" genomeDiseasePercentageProbability, exists := personDiseaseAnalysisMap[probabilityMapKey] if (exists == false) { return false, 0, 0, errors.New("Malformed person analysis map list: Missing probability of passing a variant value") } if (genomeDiseasePercentageProbability == "Unknown"){ return false, 0, 0, nil } genomeDiseaseProbabilityInt, err := helpers.ConvertStringToInt(genomeDiseasePercentageProbability) if (err != nil) { return false, 0, 0, errors.New("Malformed person analysis map list: Invalid probability of passing a variant value: " + genomeDiseasePercentageProbability) } numberOfTestedVariantsString, exists := personDiseaseAnalysisMap[genomeIdentifier + "_NumberOfVariantsTested"] if (exists == false){ return false, 0, 0, errors.New("Malformed person analysis map list: Probabilities map missing _NumberOfVariantsTested") } numberOfTestedVariants, err := helpers.ConvertStringToInt(numberOfTestedVariantsString) if (err != nil){ return false, 0, 0, errors.New("Malformed person analysis map list: Contains invalid _NumberOfVariantsTested") } return true, genomeDiseaseProbabilityInt, numberOfTestedVariants, nil } personAProbabilityIsKnown, personAWillPassVariantProbability, personANumberOfVariantsTested, err := getPersonWillPassDiseaseVariantProbability(personADiseaseAnalysisMap, personAGenomeIdentifier) if (err != nil) { return err } personBProbabilityIsKnown, personBWillPassVariantProbability, personBNumberOfVariantsTested, err := getPersonWillPassDiseaseVariantProbability(personBDiseaseAnalysisMap, personBGenomeIdentifier) if (err != nil) { return err } personANumberOfVariantsTestedString := helpers.ConvertIntToString(personANumberOfVariantsTested) personBNumberOfVariantsTestedString := helpers.ConvertIntToString(personBNumberOfVariantsTested) offspringDiseaseMap[personAGenomeIdentifier + "_NumberOfVariantsTested"] = personANumberOfVariantsTestedString offspringDiseaseMap[personBGenomeIdentifier + "_NumberOfVariantsTested"] = personBNumberOfVariantsTestedString genomePairIdentifier := personAGenomeIdentifier + "+" + personBGenomeIdentifier if (personAProbabilityIsKnown == false || personBProbabilityIsKnown == false){ offspringDiseaseMap[genomePairIdentifier + "_ProbabilityOffspringHasDisease"] = "Unknown" offspringDiseaseMap[genomePairIdentifier + "_ProbabilityOffspringHasVariant"] = "Unknown" return nil } percentageProbabilityOffspringHasDisease, percentageProbabilityOffspringHasVariant, err := GetOffspringMonogenicDiseaseProbabilities(diseaseIsDominantOrRecessive, personAWillPassVariantProbability, personBWillPassVariantProbability) if (err != nil) { return err } probabilityOffspringHasDiseaseString := helpers.ConvertIntToString(percentageProbabilityOffspringHasDisease) probabilityOffspringHasVariantString := helpers.ConvertIntToString(percentageProbabilityOffspringHasVariant) offspringDiseaseMap[genomePairIdentifier + "_ProbabilityOffspringHasDisease"] = probabilityOffspringHasDiseaseString offspringDiseaseMap[genomePairIdentifier + "_ProbabilityOffspringHasVariant"] = probabilityOffspringHasVariantString return nil } err = addPairProbabilitiesToDiseaseMap(pair1PersonAGenomeIdentifier, pair1PersonBGenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addPairProbabilitiesToDiseaseMap(pair2PersonAGenomeIdentifier, pair2PersonBGenomeIdentifier) if (err != nil) { return false, "", err } } // Now we calculate the probabilities for individual variants // We add a map for each variant offspringVariantsMapList := make([]map[string]string, 0, len(variantsList)) for _, variantObject := range variantsList{ variantIdentifier := variantObject.VariantIdentifier offspringVariantMap := map[string]string{ "ItemType": "MonogenicDiseaseVariant", "VariantIdentifier": variantIdentifier, "DiseaseName": diseaseName, } addPairProbabilitiesToVariantMap := func(personAGenomeIdentifier string, personBGenomeIdentifier string)error{ genomePairIdentifier := personAGenomeIdentifier + "+" + personBGenomeIdentifier //Outputs: // -bool: Probability is known // -float64: Probability that person will pass variant to offspring (between 0 and 1) // -error getPersonWillPassVariantProbability := func(personVariantsMapList []map[string]string, personGenomeIdentifier string)(bool, float64, error){ for _, personVariantMap := range personVariantsMapList{ itemType, exists := personVariantMap["ItemType"] if (exists == false){ return false, 0, errors.New("PersonVariantsMapList missing ItemType") } if (itemType != "MonogenicDiseaseVariant"){ return false, 0, errors.New("Person variant map is malformed: Contains non-MonogenicDiseaseVariant itemType: " + itemType) } currentVariantIdentifier, exists := personVariantMap["VariantIdentifier"] if (exists == false){ return false, 0, errors.New("Person variant map is malformed: Missing VariantIdentifier") } if (currentVariantIdentifier != variantIdentifier){ continue } personHasVariantInfo, exists := personVariantMap[personGenomeIdentifier + "_HasVariant"] if (exists == false){ return false, 0, errors.New("Person variant map is malformed: missing person genome _HasVariant value.") } if (personHasVariantInfo == "Unknown"){ return false, 0, nil } if (personHasVariantInfo == "Yes;Yes"){ return true, 1, nil } if (personHasVariantInfo == "Yes;No" || personHasVariantInfo == "No;Yes"){ return true, 0.5, nil } if (personHasVariantInfo == "No;No"){ return true, 0, nil } return false, 0, errors.New("Person variant map is malformed: Contains invalid variant info: " + personHasVariantInfo) } return false, 0, errors.New("Cannot find disease variant in personVariantsMapList: " + variantIdentifier) } // Outputs: // -bool: Probabilities are known // -int: Lower bound Percentage Probability that offspring will have 0 mutations // -int: Upper bound Percentage Probability that offspring will have 0 mutations // -int: Lower bound Percentage Probability that offspring will have 1 mutation // -int: Upper bound Percentage Probability that offspring will have 1 mutation // -int: Lower bound Percentage Probability that offspring will have 2 mutations // -int: Upper bound Percentage Probability that offspring will have 2 mutations // -error getOffspringVariantProbabilities := func()(bool, int, int, int, int, int, int, error){ //Outputs: // -int: Percentage Probability of 0 mutations // -int: Percentage Probability of 1 mutation // -int: Percentage Probability of 2 mutations // -error getOffspringProbabilities := func(personAWillPassVariantProbability float64, personBWillPassVariantProbability float64)(int, int, int, error){ // This is the probability that neither will pass variant // P = P(!A) * P(!B) probabilityOf0Mutations := (1 - personAWillPassVariantProbability) * (1 - personBWillPassVariantProbability) // This is the probability that either will pass variant, but not both // P(A XOR B) = P(A) + P(B) - (2 * P(A and B)) probabilityOf1Mutation := personAWillPassVariantProbability + personBWillPassVariantProbability - (2 * personAWillPassVariantProbability * personBWillPassVariantProbability) // This is the probability that both will pass variant // P(A and B) = P(A) * P(B) probabilityOf2Mutations := personAWillPassVariantProbability * personBWillPassVariantProbability percentageProbabilityOf0Mutations, err := helpers.FloorFloat64ToInt(probabilityOf0Mutations * 100) if (err != nil) { return 0, 0, 0, err } percentageProbabilityOf1Mutation, err := helpers.FloorFloat64ToInt(probabilityOf1Mutation * 100) if (err != nil) { return 0, 0, 0, err } percentageProbabilityOf2Mutations, err := helpers.FloorFloat64ToInt(probabilityOf2Mutations * 100) if (err != nil) { return 0, 0, 0, err } return percentageProbabilityOf0Mutations, percentageProbabilityOf1Mutation, percentageProbabilityOf2Mutations, nil } personAProbabilityIsKnown, personAWillPassVariantProbability, err := getPersonWillPassVariantProbability(personAVariantsMapList, personAGenomeIdentifier) if (err != nil) { return false, 0, 0, 0, 0, 0, 0, err } personBProbabilityIsKnown, personBWillPassVariantProbability, err := getPersonWillPassVariantProbability(personBVariantsMapList, personBGenomeIdentifier) if (err != nil) { return false, 0, 0, 0, 0, 0, 0, err } if (personAProbabilityIsKnown == false && personBProbabilityIsKnown == false){ return false, 0, 0, 0, 0, 0, 0, nil } if (personAProbabilityIsKnown == true && personBProbabilityIsKnown == false){ bestCasePersonBWillPassVariantProbability := float64(0) worstCasePersonBWillPassVariantProbability := float64(1) bestCase0MutationsProbability, bestCase1MutationProbability, bestCase2MutationsProbability, err := getOffspringProbabilities(personAWillPassVariantProbability, bestCasePersonBWillPassVariantProbability) if (err != nil) { return false, 0, 0, 0, 0, 0, 0, err } worstCase0MutationsProbability, worstCase1MutationProbability, worstCase2MutationsProbability, err := getOffspringProbabilities(personAWillPassVariantProbability, worstCasePersonBWillPassVariantProbability) if (err != nil) { return false, 0, 0, 0, 0, 0, 0, err } lowerBound1MutationProbability := min(bestCase1MutationProbability, worstCase1MutationProbability) upperBound1MutationProbability := max(bestCase1MutationProbability, worstCase1MutationProbability) return true, worstCase0MutationsProbability, bestCase0MutationsProbability, lowerBound1MutationProbability, upperBound1MutationProbability, bestCase2MutationsProbability, worstCase2MutationsProbability, nil } if (personAProbabilityIsKnown == false && personBProbabilityIsKnown == true){ bestCasePersonAWillPassVariantProbability := float64(0) worstCasePersonAWillPassVariantProbability := float64(1) bestCase0MutationsProbability, bestCase1MutationProbability, bestCase2MutationsProbability, err := getOffspringProbabilities(bestCasePersonAWillPassVariantProbability, personBWillPassVariantProbability) if (err != nil) { return false, 0, 0, 0, 0, 0, 0, err } worstCase0MutationsProbability, worstCase1MutationProbability, worstCase2MutationsProbability, err := getOffspringProbabilities(worstCasePersonAWillPassVariantProbability, personBWillPassVariantProbability) if (err != nil) { return false, 0, 0, 0, 0, 0, 0, err } lowerBound1MutationProbability := min(bestCase1MutationProbability, worstCase1MutationProbability) upperBound1MutationProbability := max(bestCase1MutationProbability, worstCase1MutationProbability) return true, worstCase0MutationsProbability, bestCase0MutationsProbability, lowerBound1MutationProbability, upperBound1MutationProbability, bestCase2MutationsProbability, worstCase2MutationsProbability, nil } offspring0MutationsProbability, offspring1MutationProbability, offspring2MutationsProbability, err := getOffspringProbabilities(personAWillPassVariantProbability, personBWillPassVariantProbability) if (err != nil) { return false, 0, 0, 0, 0, 0, 0, err } return true, offspring0MutationsProbability, offspring0MutationsProbability, offspring1MutationProbability, offspring1MutationProbability, offspring2MutationsProbability, offspring2MutationsProbability, nil } probabilitiesKnown, probabilityOf0MutationsLowerBound, probabilityOf0MutationsUpperBound, probabilityOf1MutationLowerBound, probabilityOf1MutationUpperBound, probabilityOf2MutationsLowerBound, probabilityOf2MutationsUpperBound, err := getOffspringVariantProbabilities() if (err != nil) { return err } if (probabilitiesKnown == false){ offspringVariantMap[genomePairIdentifier + "_ProbabilityOf0MutationsLowerBound"] = "Unknown" offspringVariantMap[genomePairIdentifier + "_ProbabilityOf0MutationsUpperBound"] = "Unknown" offspringVariantMap[genomePairIdentifier + "_ProbabilityOf1MutationLowerBound"] = "Unknown" offspringVariantMap[genomePairIdentifier + "_ProbabilityOf1MutationUpperBound"] = "Unknown" offspringVariantMap[genomePairIdentifier + "_ProbabilityOf2MutationsLowerBound"] = "Unknown" offspringVariantMap[genomePairIdentifier + "_ProbabilityOf2MutationsUpperBound"] = "Unknown" return nil } probabilityOf0MutationsLowerBoundString := helpers.ConvertIntToString(probabilityOf0MutationsLowerBound) probabilityOf0MutationsUpperBoundString := helpers.ConvertIntToString(probabilityOf0MutationsUpperBound) probabilityOf1MutationLowerBoundString := helpers.ConvertIntToString(probabilityOf1MutationLowerBound) probabilityOf1MutationUpperBoundString := helpers.ConvertIntToString(probabilityOf1MutationUpperBound) probabilityOf2MutationsLowerBoundString := helpers.ConvertIntToString(probabilityOf2MutationsLowerBound) probabilityOf2MutationsUpperBoundString := helpers.ConvertIntToString(probabilityOf2MutationsUpperBound) offspringVariantMap[genomePairIdentifier + "_ProbabilityOf0MutationsLowerBound"] = probabilityOf0MutationsLowerBoundString offspringVariantMap[genomePairIdentifier + "_ProbabilityOf0MutationsUpperBound"] = probabilityOf0MutationsUpperBoundString offspringVariantMap[genomePairIdentifier + "_ProbabilityOf1MutationLowerBound"] = probabilityOf1MutationLowerBoundString offspringVariantMap[genomePairIdentifier + "_ProbabilityOf1MutationUpperBound"] = probabilityOf1MutationUpperBoundString offspringVariantMap[genomePairIdentifier + "_ProbabilityOf2MutationsLowerBound"] = probabilityOf2MutationsLowerBoundString offspringVariantMap[genomePairIdentifier + "_ProbabilityOf2MutationsUpperBound"] = probabilityOf2MutationsUpperBoundString return nil } err = addPairProbabilitiesToVariantMap(pair1PersonAGenomeIdentifier, pair1PersonBGenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addPairProbabilitiesToVariantMap(pair2PersonAGenomeIdentifier, pair2PersonBGenomeIdentifier) if (err != nil) { return false, "", err } } offspringVariantsMapList = append(offspringVariantsMapList, offspringVariantMap) } // Now we check for conflicts in monogenic disease results // For couples, a conflict is when either genome pair has differing results for any disease probability/variant probability // This is different from a person having conflicts within their different genomes // Each genome pair only uses 1 genome from each person. if (genomePair2Exists == true){ // Conflicts are only possible if two genome pairs exist checkIfConflictExists := func()(bool, error){ genomePair1Identifier := pair1PersonAGenomeIdentifier + "+" + pair1PersonBGenomeIdentifier genomePair2Identifier := pair2PersonAGenomeIdentifier + "+" + pair2PersonBGenomeIdentifier allGenomePairIdentifiersList := []string{genomePair1Identifier, genomePair2Identifier} probabilityOffspringHasDisease := "" probabilityOffspringHasVariant := "" for index, genomePairIdentifier := range allGenomePairIdentifiersList{ currentProbabilityOffspringHasDisease, exists := offspringDiseaseMap[genomePairIdentifier + "_ProbabilityOffspringHasDisease"] if (exists == false) { return false, errors.New("Cannot find _ProbabilityOffspringHasDisease key when searching for conflicts.") } currentProbabilityOffspringHasVariant, exists := offspringDiseaseMap[genomePairIdentifier + "_ProbabilityOffspringHasVariant"] if (exists == false) { return false, errors.New("Cannot find _ProbabilityOffspringHasVariant key when searching for conflicts.") } if (index == 0){ probabilityOffspringHasDisease = currentProbabilityOffspringHasDisease probabilityOffspringHasVariant = currentProbabilityOffspringHasVariant continue } if (currentProbabilityOffspringHasDisease != probabilityOffspringHasDisease){ return true, nil } if (currentProbabilityOffspringHasVariant != probabilityOffspringHasVariant){ return true, nil } } for _, variantMap := range offspringVariantsMapList{ probabilityOf0MutationsLowerBound := "" probabilityOf0MutationsUpperBound := "" probabilityOf1MutationLowerBound := "" probabilityOf1MutationUpperBound := "" probabilityOf2MutationsLowerBound := "" probabilityOf2MutationsUpperBound := "" for index, genomePairIdentifier := range allGenomePairIdentifiersList{ currentProbabilityOf0MutationsLowerBound, exists := variantMap[genomePairIdentifier + "_ProbabilityOf0MutationsLowerBound"] if (exists == false){ return false, errors.New("Cannot find _ProbabilityOf0MutationsLowerBound key when searching for conflicts.") } currentProbabilityOf0MutationsUpperBound, exists := variantMap[genomePairIdentifier + "_ProbabilityOf0MutationsUpperBound"] if (exists == false){ return false, errors.New("Cannot find _ProbabilityOf0MutationsUpperBound key when searching for conflicts.") } currentProbabilityOf1MutationLowerBound, exists := variantMap[genomePairIdentifier + "_ProbabilityOf1MutationLowerBound"] if (exists == false){ return false, errors.New("Cannot find _ProbabilityOf1MutationLowerBound key when searching for conflicts.") } currentProbabilityOf1MutationUpperBound, exists := variantMap[genomePairIdentifier + "_ProbabilityOf1MutationUpperBound"] if (exists == false){ return false, errors.New("Cannot find _ProbabilityOf1MutationUpperBound key when searching for conflicts.") } currentProbabilityOf2MutationsLowerBound, exists := variantMap[genomePairIdentifier + "_ProbabilityOf2MutationsLowerBound"] if (exists == false){ return false, errors.New("Cannot find _ProbabilityOf2MutationsLowerBound key when searching for conflicts.") } currentProbabilityOf2MutationsUpperBound, exists := variantMap[genomePairIdentifier + "_ProbabilityOf2MutationsUpperBound"] if (exists == false){ return false, errors.New("Cannot find _ProbabilityOf2MutationsUpperBound key when searching for conflicts.") } if (index == 0){ probabilityOf0MutationsLowerBound = currentProbabilityOf0MutationsLowerBound probabilityOf0MutationsUpperBound = currentProbabilityOf0MutationsUpperBound probabilityOf1MutationLowerBound = currentProbabilityOf1MutationLowerBound probabilityOf1MutationUpperBound = currentProbabilityOf1MutationUpperBound probabilityOf2MutationsLowerBound = currentProbabilityOf2MutationsLowerBound probabilityOf2MutationsUpperBound = currentProbabilityOf2MutationsUpperBound continue } if (currentProbabilityOf0MutationsLowerBound != probabilityOf0MutationsLowerBound){ return true, nil } if (currentProbabilityOf0MutationsUpperBound != probabilityOf0MutationsUpperBound){ return true, nil } if (currentProbabilityOf1MutationLowerBound != probabilityOf1MutationLowerBound){ return true, nil } if (currentProbabilityOf1MutationUpperBound != probabilityOf1MutationUpperBound){ return true, nil } if (currentProbabilityOf2MutationsLowerBound != probabilityOf2MutationsLowerBound){ return true, nil } if (currentProbabilityOf2MutationsUpperBound != probabilityOf2MutationsUpperBound){ return true, nil } } } return false, nil } conflictExists, err := checkIfConflictExists() if (err != nil) { return false, "", err } conflictExistsString := helpers.ConvertBoolToYesOrNoString(conflictExists) offspringDiseaseMap["ConflictExists"] = conflictExistsString } newAnalysisMapList = append(newAnalysisMapList, offspringDiseaseMap) newAnalysisMapList = append(newAnalysisMapList, offspringVariantsMapList...) } // Step 2: Polygenic Diseases polygenicDiseaseObjectsList, err := polygenicDiseases.GetPolygenicDiseaseObjectsList() if (err != nil){ return false, "", err } for _, diseaseObject := range polygenicDiseaseObjectsList{ diseaseName := diseaseObject.DiseaseName diseaseLociList := diseaseObject.LociList _, personALociMapList, err := getPolygenicDiseaseAnalysis(personAGenomesWithMetadataList, diseaseObject) if (err != nil) { return false, "", err } _, personBLociMapList, err := getPolygenicDiseaseAnalysis(personBGenomesWithMetadataList, diseaseObject) if (err != nil) { return false, "", err } offspringLociMapList := make([]map[string]string, 0, len(diseaseLociList)) for _, locusObject := range diseaseLociList{ locusIdentifier := locusObject.LocusIdentifier locusRiskWeightsMap := locusObject.RiskWeightsMap locusOddsRatiosMap := locusObject.OddsRatiosMap offspringLocusMap := map[string]string{ "ItemType": "PolygenicDiseaseLocus", "LocusIdentifier": locusIdentifier, "DiseaseName": diseaseName, } // This will calculate the average risk weight and probability change for the offspring from the two specified genomes // It then adds the pair entry to the diseaseMap addPairDiseaseLocusInfoToLocusMap := func(personAGenomeIdentifier string, personBGenomeIdentifier string)error{ getPersonGenomeLocusBasePair := func(personGenomeIdentifier string, personLociMapList []map[string]string)(bool, string, error){ for _, locusMap := range personLociMapList{ itemType, exists := locusMap["ItemType"] if (exists == false){ return false, "", errors.New("getPolygenicDiseaseAnalysis returning lociMapList with item missing ItemType") } if (itemType != "PolygenicDiseaseLocus"){ return false, "", errors.New("getPolygenicDiseaseAnalysis returning lociMapList with non-PolygenicDiseaseLocus item: " + itemType) } currentLocusIdentifier, exists := locusMap["LocusIdentifier"] if (exists == false){ return false, "", errors.New("getPolygenicDiseaseAnalysis returning lociMapList with item missing LocusIdentifier") } if (currentLocusIdentifier != locusIdentifier){ continue } locusBasePair, exists := locusMap[personGenomeIdentifier + "_LocusBasePair"] if (exists == false){ return false, "", errors.New("getPolygenicDiseaseAnalysis returning lociMapList with item missing _LocusBasePair") } if (locusBasePair == "Unknown"){ return false, "", nil } return true, locusBasePair, nil } return false, "", errors.New("getPolygenicDiseaseAnalysis returning lociMapList missing a locus: " + locusIdentifier) } genomePairIdentifier := personAGenomeIdentifier + "+" + personBGenomeIdentifier personALocusBasePairKnown, personALocusBasePair, err := getPersonGenomeLocusBasePair(personAGenomeIdentifier, personALociMapList) if (err != nil) { return err } personBLocusBasePairKnown, personBLocusBasePair, err := getPersonGenomeLocusBasePair(personBGenomeIdentifier, personBLociMapList) if (err != nil) { return err } if (personALocusBasePairKnown == false || personBLocusBasePairKnown == false){ offspringLocusMap[genomePairIdentifier + "_OffspringRiskWeight"] = "Unknown" offspringLocusMap[genomePairIdentifier + "_OffspringOddsRatio"] = "Unknown" offspringLocusMap[genomePairIdentifier + "_OffspringUnknownOddsRatiosWeightSum"] = "Unknown" return nil } offspringAverageRiskWeight, offspringOddsRatioKnown, offspringAverageOddsRatio, averageUnknownOddsRatiosWeightSum, err := GetOffspringPolygenicDiseaseLocusInfo(locusRiskWeightsMap, locusOddsRatiosMap, personALocusBasePair, personBLocusBasePair) if (err != nil) { return err } averageRiskWeightString := helpers.ConvertIntToString(offspringAverageRiskWeight) offspringLocusMap[genomePairIdentifier + "_OffspringRiskWeight"] = averageRiskWeightString if (offspringOddsRatioKnown == false){ offspringLocusMap[genomePairIdentifier + "_OffspringOddsRatio"] = "Unknown" } else { averageOddsRatioString := helpers.ConvertFloat64ToString(offspringAverageOddsRatio) offspringLocusMap[genomePairIdentifier + "_OffspringOddsRatio"] = averageOddsRatioString } averageUnknownOddsRatiosWeightSumString := helpers.ConvertIntToString(averageUnknownOddsRatiosWeightSum) offspringLocusMap[genomePairIdentifier + "_OffspringUnknownOddsRatiosWeightSum"] = averageUnknownOddsRatiosWeightSumString return nil } err = addPairDiseaseLocusInfoToLocusMap(pair1PersonAGenomeIdentifier, pair1PersonBGenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addPairDiseaseLocusInfoToLocusMap(pair2PersonAGenomeIdentifier, pair2PersonBGenomeIdentifier) if (err != nil) { return false, "", err } } offspringLociMapList = append(offspringLociMapList, offspringLocusMap) } // Now we iterate through loci to determine genome pair disease info offspringPolygenicDiseaseMap := map[string]string{ "ItemType": "PolygenicDisease", "DiseaseName": diseaseName, } polygenicDiseaseLociMap, err := polygenicDiseases.GetPolygenicDiseaseLociMap(diseaseName) if (err != nil) { return false, "", err } addPairPolygenicDiseaseInfoToDiseaseMap := func(personAGenomeIdentifier string, personBGenomeIdentifier string)error{ genomePairIdentifier := personAGenomeIdentifier + "+" + personBGenomeIdentifier numberOfLociTested := 0 summedRiskWeights := 0 minimumPossibleRiskWeightSum := 0 maximumPossibleRiskWeightSum := 0 for _, locusMap := range offspringLociMapList{ locusIdentifier, exists := locusMap["LocusIdentifier"] if (exists == false){ return errors.New("offspringLociMapList item missing LocusIdentifier") } locusObject, exists := polygenicDiseaseLociMap[locusIdentifier] if (exists == false){ return errors.New("offspringLociMapList contains locus not found in allDiseaseLociObjectsMap") } locusMinimumRiskWeight := locusObject.MinimumRiskWeight locusMaximumRiskWeight := locusObject.MaximumRiskWeight minimumPossibleRiskWeightSum += locusMinimumRiskWeight maximumPossibleRiskWeightSum += locusMaximumRiskWeight locusRiskWeight, exists := locusMap[genomePairIdentifier + "_OffspringRiskWeight"] if (exists == false){ return errors.New("offspringLociMapList contains locusMap missing _OffspringRiskWeight") } if (locusRiskWeight == "Unknown"){ continue } numberOfLociTested += 1 locusRiskWeightInt, err := helpers.ConvertStringToInt(locusRiskWeight) if (err != nil){ return errors.New("offspringLociMapList contains locusMap with invalid _OffspringRiskWeight: " + locusRiskWeight) } summedRiskWeights += locusRiskWeightInt } if (numberOfLociTested == 0){ // No loci were tested offspringPolygenicDiseaseMap[genomePairIdentifier + "_NumberOfLociTested"] = "0" offspringPolygenicDiseaseMap[genomePairIdentifier + "_OffspringRiskScore"] = "Unknown" return nil } numberOfLociTestedString := helpers.ConvertIntToString(numberOfLociTested) offspringPolygenicDiseaseMap[genomePairIdentifier + "_NumberOfLociTested"] = numberOfLociTestedString pairDiseaseRiskScore, err := helpers.ScaleNumberProportionally(true, summedRiskWeights, minimumPossibleRiskWeightSum, maximumPossibleRiskWeightSum, 0, 10) if (err != nil) { return err } pairDiseaseRiskScoreString := helpers.ConvertIntToString(pairDiseaseRiskScore) offspringPolygenicDiseaseMap[genomePairIdentifier + "_OffspringRiskScore"] = pairDiseaseRiskScoreString return nil } addPairPolygenicDiseaseInfoToDiseaseMap(pair1PersonAGenomeIdentifier, pair1PersonBGenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addPairPolygenicDiseaseInfoToDiseaseMap(pair2PersonAGenomeIdentifier, pair2PersonBGenomeIdentifier) if (err != nil) { return false, "", err } } if (genomePair2Exists == true){ // We check for conflicts // Conflicts are only possible if two genome pairs exist checkIfConflictExists := func()(bool, error){ genomePair1Identifier := pair1PersonAGenomeIdentifier + "+" + pair1PersonBGenomeIdentifier genomePair2Identifier := pair2PersonAGenomeIdentifier + "+" + pair2PersonBGenomeIdentifier allGenomePairIdentifiersList := []string{genomePair1Identifier, genomePair2Identifier} for _, locusMap := range offspringLociMapList{ offspringRiskWeight := "" offspringOddsRatio := "" offspringUnknownOddsRatiosWeightSum := "" for index, genomePairIdentifier := range allGenomePairIdentifiersList{ currentOffspringRiskWeight, exists := locusMap[genomePairIdentifier + "_OffspringRiskWeight"] if (exists == false){ return false, errors.New("Cannot find _OffspringRiskWeight key when searching for conflicts.") } currentOffspringOddsRatio, exists := locusMap[genomePairIdentifier + "_OffspringOddsRatio"] if (exists == false){ return false, errors.New("Cannot find _OffspringOddsRatio key when searching for conflicts.") } currentOffspringUnknownOddsRatiosWeightSum, exists := locusMap[genomePairIdentifier + "_OffspringUnknownOddsRatiosWeightSum"] if (exists == false){ return false, errors.New("Cannot find _OffspringUnknownOddsRatiosWeightSum key when searching for conflicts.") } if (index == 0){ offspringRiskWeight = currentOffspringRiskWeight offspringOddsRatio = currentOffspringOddsRatio offspringUnknownOddsRatiosWeightSum = currentOffspringUnknownOddsRatiosWeightSum continue } if (currentOffspringRiskWeight != offspringRiskWeight){ return true, nil } if (currentOffspringOddsRatio != offspringOddsRatio){ return true, nil } if (currentOffspringUnknownOddsRatiosWeightSum != offspringUnknownOddsRatiosWeightSum){ return true, nil } } } return false, nil } conflictExists, err := checkIfConflictExists() if (err != nil) { return false, "", err } conflictExistsString := helpers.ConvertBoolToYesOrNoString(conflictExists) offspringPolygenicDiseaseMap["ConflictExists"] = conflictExistsString } newAnalysisMapList = append(newAnalysisMapList, offspringPolygenicDiseaseMap) newAnalysisMapList = append(newAnalysisMapList, offspringLociMapList...) } // Step 3: Traits traitObjectsList, err := traits.GetTraitObjectsList() if (err != nil) { return false, "", err } for _, traitObject := range traitObjectsList{ traitName := traitObject.TraitName traitRulesList := traitObject.RulesList _, personATraitRulesMapList, err := getTraitAnalysis(personAGenomesWithMetadataList, traitObject) if (err != nil) { return false, "", err } _, personBTraitRulesMapList, err := getTraitAnalysis(personBGenomesWithMetadataList, traitObject) if (err != nil) { return false, "", err } offspringRulesMapList := make([]map[string]string, 0, len(traitRulesList)) for _, ruleObject := range traitRulesList{ ruleIdentifier := ruleObject.RuleIdentifier // This is a list that describes the locus rsids and their values that must be fulfilled to pass the rule ruleLocusObjectsList := ruleObject.LociList offspringRuleMap := map[string]string{ "ItemType": "TraitRule", "RuleIdentifier": ruleIdentifier, "TraitName": traitName, } // This will calculate the number of offspring that will pass the rule for the offspring from the two specified genomes // It then adds the pair entry to the ruleMap addPairTraitRuleInfoToRuleMap := func(personAGenomeIdentifier string, personBGenomeIdentifier string)error{ //Outputs: // -bool: All rule loci are known // -map[string]string: RSID -> Locus base pair value // -error getPersonGenomeRuleLociValuesMap := func(personGenomeIdentifier string, personTraitRulesMapList []map[string]string)(bool, map[string]string, error){ for _, ruleMap := range personTraitRulesMapList{ itemType, exists := ruleMap["ItemType"] if (exists == false){ return false, nil, errors.New("getTraitAnalysis returning rulesMapList with item missing ItemType") } if (itemType != "TraitRule"){ return false, nil, errors.New("getTraitAnalysis returning rulesMapList with non-TraitRule item: " + itemType) } currentRuleIdentifier, exists := ruleMap["RuleIdentifier"] if (exists == false){ return false, nil, errors.New("getTraitAnalysis returning rulesMapList with item missing RuleIdentifier") } if (currentRuleIdentifier != ruleIdentifier){ continue } // Map structure: Locus identifier -> Locus base pair value personTraitLociBasePairsMap := make(map[string]string) for _, locusObject := range ruleLocusObjectsList{ locusIdentifier := locusObject.LocusIdentifier personLocusBasePair, exists := ruleMap[personGenomeIdentifier + "_RuleLocusBasePair_" + locusIdentifier] if (exists == false){ return false, nil, errors.New("_LocusBasePair_ value not found in person trait analysis ruleMap") } if (personLocusBasePair == "Unknown"){ // Not all locus values are known, thus, we cannot determine if the genome passes the rule return false, nil, nil } personTraitLociBasePairsMap[locusIdentifier] = personLocusBasePair } return true, personTraitLociBasePairsMap, nil } return false, nil, errors.New("getTraitAnalysis returning rulesMapList missing a rule: " + ruleIdentifier) } allPersonALociKnown, personAGenomeRuleLociValuesMap, err := getPersonGenomeRuleLociValuesMap(personAGenomeIdentifier, personATraitRulesMapList) if (err != nil) { return err } allPersonBLociKnown, personBGenomeRuleLociValuesMap, err := getPersonGenomeRuleLociValuesMap(personBGenomeIdentifier, personBTraitRulesMapList) if (err != nil) { return err } genomePairIdentifier := personAGenomeIdentifier + "+" + personBGenomeIdentifier if (allPersonALociKnown == false || allPersonBLociKnown == false){ // We only know how many of the 4 prospective offspring pass the rule if all loci are known for both people's genomes offspringRuleMap[genomePairIdentifier + "_OffspringProbabilityOfPassingRule"] = "Unknown" return nil } // This is a probability between 0 and 1 offspringProbabilityOfPassingRule := float64(1) for _, ruleLocusObject := range ruleLocusObjectsList{ locusIdentifier := ruleLocusObject.LocusIdentifier personALocusBasePair, exists := personAGenomeRuleLociValuesMap[locusIdentifier] if (exists == false){ return errors.New("personAGenomeRuleLociValuesMap missing locusIdentifier") } personBLocusBasePair, exists := personBGenomeRuleLociValuesMap[locusIdentifier] if (exists == false){ return errors.New("personBGenomeRuleLociValuesMap missing locusIdentifier") } locusRequiredBasePairsList := ruleLocusObject.BasePairsList offspringProbabilityOfPassingRuleLocus, err := GetOffspringTraitRuleLocusInfo(locusRequiredBasePairsList, personALocusBasePair, personBLocusBasePair) if (err != nil) { return err } offspringProbabilityOfPassingRule *= offspringProbabilityOfPassingRuleLocus } offspringPercentageProbabilityOfPassingRule := offspringProbabilityOfPassingRule * 100 offspringPercentageProbabilityOfPassingRuleString := helpers.ConvertFloat64ToStringRounded(offspringPercentageProbabilityOfPassingRule, 0) offspringRuleMap[genomePairIdentifier + "_OffspringProbabilityOfPassingRule"] = offspringPercentageProbabilityOfPassingRuleString return nil } err = addPairTraitRuleInfoToRuleMap(pair1PersonAGenomeIdentifier, pair1PersonBGenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addPairTraitRuleInfoToRuleMap(pair2PersonAGenomeIdentifier, pair2PersonBGenomeIdentifier) if (err != nil) { return false, "", err } } offspringRulesMapList = append(offspringRulesMapList, offspringRuleMap) } // Now we iterate through rules to determine genome pair trait info offspringTraitMap := map[string]string{ "ItemType": "Trait", "TraitName": traitName, } traitRulesMap, err := traits.GetTraitRulesMap(traitName) if (err != nil) { return false, "", err } addPairTraitInfoToTraitMap := func(personAGenomeIdentifier string, personBGenomeIdentifier string)error{ genomePairIdentifier := personAGenomeIdentifier + "+" + personBGenomeIdentifier if (len(traitRulesList) == 0){ // This trait is not yet able to be analyzed. // This feature will be added soon // These traits will use neural networks instead of rules offspringTraitMap[genomePairIdentifier + "_NumberOfRulesTested"] = "0" return nil } numberOfRulesTested := 0 // Map Structure: Outcome -> Average score from all rules averageOutcomeScoresMap := make(map[string]float64) for _, ruleMap := range offspringRulesMapList{ ruleIdentifier, exists := ruleMap["RuleIdentifier"] if (exists == false){ return errors.New("offspringRulesMapList item missing RuleIdentifier") } offspringPercentageProbabilityOfPassingRule, exists := ruleMap[genomePairIdentifier + "_OffspringProbabilityOfPassingRule"] if (exists == false){ return errors.New("offspringRulesMapList contains ruleMap missing _OffspringProbabilityOfPassingRule") } if (offspringPercentageProbabilityOfPassingRule == "Unknown"){ continue } numberOfRulesTested += 1 offspringPercentageProbabilityOfPassingRuleFloat64, err := helpers.ConvertStringToFloat64(offspringPercentageProbabilityOfPassingRule) if (err != nil){ return errors.New("offspringRulesMapList contains ruleMap with invalid _OffspringProbabilityOfPassingRule: " + offspringPercentageProbabilityOfPassingRule) } // This is the 0 - 1 probability value offspringProbabilityOfPassingRule := offspringPercentageProbabilityOfPassingRuleFloat64/100 ruleObject, exists := traitRulesMap[ruleIdentifier] if (exists == false){ return errors.New("offspringRulesMapList contains rule not found in traitRulesMap: " + ruleIdentifier) } ruleOutcomePointsMap := ruleObject.OutcomePointsMap for outcomeName, outcomePointsEffect := range ruleOutcomePointsMap{ pointsToAdd := float64(outcomePointsEffect) * offspringProbabilityOfPassingRule averageOutcomeScoresMap[outcomeName] += pointsToAdd } } numberOfRulesTestedString := helpers.ConvertIntToString(numberOfRulesTested) offspringTraitMap[genomePairIdentifier + "_NumberOfRulesTested"] = numberOfRulesTestedString if (numberOfRulesTested == 0){ return nil } traitOutcomesList := traitObject.OutcomesList for _, outcomeName := range traitOutcomesList{ getOffspringAverageOutcomeScore := func()float64{ averageOutcomeScore, exists := averageOutcomeScoresMap[outcomeName] if (exists == false){ // No rules effected this outcome. return 0 } return averageOutcomeScore } offspringAverageOutcomeScore := getOffspringAverageOutcomeScore() offspringAverageOutcomeScoreString := helpers.ConvertFloat64ToStringRounded(offspringAverageOutcomeScore, 2) offspringTraitMap[genomePairIdentifier + "_OffspringAverageOutcomeScore_" + outcomeName] = offspringAverageOutcomeScoreString } return nil } addPairTraitInfoToTraitMap(pair1PersonAGenomeIdentifier, pair1PersonBGenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addPairTraitInfoToTraitMap(pair2PersonAGenomeIdentifier, pair2PersonBGenomeIdentifier) if (err != nil) { return false, "", err } } if (genomePair2Exists == true){ // We check for conflicts // Conflicts are only possible if two genome pairs exist checkIfConflictExists := func()(bool, error){ genomePair1Identifier := pair1PersonAGenomeIdentifier + "+" + pair1PersonBGenomeIdentifier genomePair2Identifier := pair2PersonAGenomeIdentifier + "+" + pair2PersonBGenomeIdentifier allGenomePairIdentifiersList := []string{genomePair1Identifier, genomePair2Identifier} for _, ruleMap := range offspringRulesMapList{ offspringProbabilityOfPassingRule := "" for index, genomePairIdentifier := range allGenomePairIdentifiersList{ currentOffspringProbabilityOfPassingRule, exists := ruleMap[genomePairIdentifier + "_OffspringProbabilityOfPassingRule"] if (exists == false){ return false, errors.New("Cannot find _OffspringProbabilityOfPassingRule key when searching for conflicts.") } if (index == 0){ offspringProbabilityOfPassingRule = currentOffspringProbabilityOfPassingRule continue } if (currentOffspringProbabilityOfPassingRule != offspringProbabilityOfPassingRule){ return true, nil } } } return false, nil } conflictExists, err := checkIfConflictExists() if (err != nil) { return false, "", err } conflictExistsString := helpers.ConvertBoolToYesOrNoString(conflictExists) offspringTraitMap["ConflictExists"] = conflictExistsString } newAnalysisMapList = append(newAnalysisMapList, offspringTraitMap) newAnalysisMapList = append(newAnalysisMapList, offspringRulesMapList...) } analysisBytes, err := json.MarshalIndent(newAnalysisMapList, "", "\t") if (err != nil) { return false, "", err } analysisString := string(analysisBytes) return true, analysisString, nil } // We also use this function when calculating offspring probabilities between users in viewProfileGui.go //Outputs: // -int: Percentage probability offspring has disease (0-100) // -int: Percentage probability offspring has variant (0-100) // -error func GetOffspringMonogenicDiseaseProbabilities(dominantOrRecessive string, personAWillPassVariantPercentageProbability int, personBWillPassVariantPercentageProbability int)(int, int, error){ if (dominantOrRecessive != "Dominant" && dominantOrRecessive != "Recessive"){ return 0, 0, errors.New("GetOffspringMonogenicDiseaseProbabilities called with invalid dominantOrRecessive: " + dominantOrRecessive) } personAWillPassVariantProbability := float64(personAWillPassVariantPercentageProbability)/100 personBWillPassVariantProbability := float64(personBWillPassVariantPercentageProbability)/100 if (personAWillPassVariantProbability < 0 || personAWillPassVariantProbability > 1){ return 0, 0, errors.New("GetOffspringMonogenicDiseaseProbabilities called with invalid personAWillPassVariantProbability") } if (personBWillPassVariantProbability < 0 || personBWillPassVariantProbability > 1){ return 0, 0, errors.New("GetOffspringMonogenicDiseaseProbabilities called with invalid personBWillPassVariantProbability") } // The probability offspring has a variant = the probability that either parent passes a variant (inclusive or) // We find the probability of the offspring having a monogenic disease variant as follows: // P(A U B) = P(A) + P(B) - P(A ∩ B) // (Probability of person A passing a variant) + (Probability of person B passing a variant) - (Probability of offspring having disease) // A person with a variant may have the disease, or just be a carrier. probabilityOffspringHasVariant := personAWillPassVariantProbability + personBWillPassVariantProbability - (personAWillPassVariantProbability * personBWillPassVariantProbability) if (dominantOrRecessive == "Dominant"){ // The probability of having the monogenic disease is the same as the probability of having a variant percentageProbabilityOffspringHasVariant := int(probabilityOffspringHasVariant * 100) return percentageProbabilityOffspringHasVariant, percentageProbabilityOffspringHasVariant, nil } // We find the probability of the offspring having the mongenic disease as follows: // P(A and B) = P(A) * P(B) // (Probability of person A Passing a variant) * (Probability of person B passing a variant) probabilityOffspringHasDisease := personAWillPassVariantProbability * personBWillPassVariantProbability percentageProbabilityOffspringHasDisease := probabilityOffspringHasDisease * 100 percentageProbabilityOffspringHasVariant := probabilityOffspringHasVariant * 100 // This conversion remove any digits after the radix point // This will not result in any false 0% values, an example being 0.9% becoming 0% // This is because the lowest non-zero probability a person can have for passing a variant is 50% // Thus, the lowest non-zero probability of an offspring having a disease is 25% percentageProbabilityOffspringHasDiseaseInt := int(percentageProbabilityOffspringHasDisease) percentageProbabilityOffspringHasVariantInt := int(percentageProbabilityOffspringHasVariant) return percentageProbabilityOffspringHasDiseaseInt, percentageProbabilityOffspringHasVariantInt, nil } //Outputs: // -int: Offspring average risk weight // -bool: Odds ratio is known // -float64: Offspring average odds ratio // -int: Offspring unknown odds ratios weight sum // -error func GetOffspringPolygenicDiseaseLocusInfo(locusRiskWeightsMap map[string]int, locusOddsRatiosMap map[string]float64, personALocusBasePair string, personBLocusBasePair string)(int, bool, float64, int, error){ personABase1, personABase2, delimiterFound := strings.Cut(personALocusBasePair, ";") if (delimiterFound == false){ return 0, false, 0, 0, errors.New("GetOffspringPolygenicDiseaseLocusInfo called with invalid person A locus base pair: " + personALocusBasePair) } personBBase1, personBBase2, delimiterFound := strings.Cut(personBLocusBasePair, ";") if (delimiterFound == false){ return 0, false, 0, 0, errors.New("GetOffspringPolygenicDiseaseLocusInfo called with invalid person B locus base pair: " + personBLocusBasePair) } // We create the 4 options for the offspring's bases at this locus offspringBasePairOutcomeA := personABase1 + ";" + personBBase1 offspringBasePairOutcomeB := personABase2 + ";" + personBBase2 offspringBasePairOutcomeC := personABase1 + ";" + personBBase2 offspringBasePairOutcomeD := personABase2 + ";" + personBBase1 baseOutcomesList := []string{offspringBasePairOutcomeA, offspringBasePairOutcomeB, offspringBasePairOutcomeC, offspringBasePairOutcomeD} summedRiskWeight := 0 numberOfSummedOddsRatios := 0 summedOddsRatios := float64(0) numberOfSummedUnknownOddsRatioWeights := 0 summedUnknownOddsRatioWeights := 0 for _, outcomeBasePair := range baseOutcomesList{ isValid := verifyBasePair(outcomeBasePair) if (isValid == false){ return 0, false, 0, 0, errors.New("GetOffspringPolygenicDiseaseLocusInfo called with invalid locus base pair: " + outcomeBasePair) } offspringOutcomeRiskWeight, exists := locusRiskWeightsMap[outcomeBasePair] if (exists == false){ // We do not know the risk weight for this base pair // We treat this as a 0 risk for both weight and odds ratio summedOddsRatios += 1 numberOfSummedOddsRatios += 1 continue } summedRiskWeight += offspringOutcomeRiskWeight offspringOutcomeOddsRatio, exists := locusOddsRatiosMap[outcomeBasePair] if (exists == false){ // This particular outcome has no known odds ratio // We add it to the unknown odds ratio weights sum summedUnknownOddsRatioWeights += offspringOutcomeRiskWeight numberOfSummedUnknownOddsRatioWeights += 1 } else { summedOddsRatios += offspringOutcomeOddsRatio numberOfSummedOddsRatios += 1 } } averageRiskWeight := summedRiskWeight/4 getAverageUnknownOddsRatiosWeightSum := func()int{ if (numberOfSummedUnknownOddsRatioWeights == 0){ return 0 } averageUnknownOddsRatiosWeightSum := summedUnknownOddsRatioWeights/numberOfSummedUnknownOddsRatioWeights return averageUnknownOddsRatiosWeightSum } averageUnknownOddsRatiosWeightSum := getAverageUnknownOddsRatiosWeightSum() if (numberOfSummedOddsRatios == 0){ return averageRiskWeight, false, 0, averageUnknownOddsRatiosWeightSum, nil } averageOddsRatio := summedOddsRatios/float64(numberOfSummedOddsRatios) return averageRiskWeight, true, averageOddsRatio, averageUnknownOddsRatiosWeightSum, nil } //Outputs: // -float64: Probability of offspring passing rule (0-1) // -error func GetOffspringTraitRuleLocusInfo(locusRequiredBasePairsList []string, personALocusBasePair string, personBLocusBasePair string)(float64, error){ personABase1, personABase2, delimiterFound := strings.Cut(personALocusBasePair, ";") if (delimiterFound == false){ return 0, errors.New("GetOffspringTraitRuleLocusInfo called with invalid personA locus base pair: " + personALocusBasePair) } personBBase1, personBBase2, delimiterFound := strings.Cut(personBLocusBasePair, ";") if (delimiterFound == false){ return 0, errors.New("GetOffspringTraitRuleLocusInfo called with invalid personB locus base pair: " + personBLocusBasePair) } // We create the 4 options for the offspring's bases at this locus offspringBasePairOutcomeA := personABase1 + ";" + personBBase1 offspringBasePairOutcomeB := personABase2 + ";" + personBBase2 offspringBasePairOutcomeC := personABase1 + ";" + personBBase2 offspringBasePairOutcomeD := personABase2 + ";" + personBBase1 baseOutcomesList := []string{offspringBasePairOutcomeA, offspringBasePairOutcomeB, offspringBasePairOutcomeC, offspringBasePairOutcomeD} numberOfOffspringOutcomesWhomPassRuleLocus := 0 for _, outcomeBasePair := range baseOutcomesList{ isValid := verifyBasePair(outcomeBasePair) if (isValid == false){ return 0, errors.New("GetOffspringTraitRuleLocusInfo called with invalid locus base pair: " + outcomeBasePair) } outcomePassesRuleLocus := slices.Contains(locusRequiredBasePairsList, outcomeBasePair) if (outcomePassesRuleLocus == true){ numberOfOffspringOutcomesWhomPassRuleLocus += 1 } } offspringProbabilityOfPassingRuleLocus := float64(numberOfOffspringOutcomesWhomPassRuleLocus)/float64(4) return offspringProbabilityOfPassingRuleLocus, nil } // This function will retrieve the base pair of the locus from the input genome map // We need this because each rsID has aliases, so we must sometimes check those aliases to find locus values // // Outputs: // -bool: Valid base pair value found // -string: Base 1 Value (Nucleotide base for the SNP) // -string: Base 2 Value (Nucleotide base for the SNP) // -bool: Base pairs are phased // -error func getGenomeLocusBasePair(inputGenomeMap map[int64]locusValue.LocusValue, locusRSID int64)(bool, string, string, bool, error){ // Outputs: // -bool: Locus value found // -locusValue.LocusValue // -error getLocusValue := func()(bool, locusValue.LocusValue, error){ currentLocusValue, exists := inputGenomeMap[locusRSID] if (exists == true){ return true, currentLocusValue, nil } // We check for aliases anyAliasesExist, rsidAliasesList, err := locusMetadata.GetRSIDAliases(locusRSID) if (err != nil) { return false, locusValue.LocusValue{}, err } if (anyAliasesExist == false){ return false, locusValue.LocusValue{}, nil } for _, rsidAlias := range rsidAliasesList{ currentLocusValue, exists := inputGenomeMap[rsidAlias] if (exists == true){ return true, currentLocusValue, nil } } return false, locusValue.LocusValue{}, nil } locusValueFound, locusValueObject, err := getLocusValue() if (err != nil) { return false, "", "", false, err } if (locusValueFound == false){ return false, "", "", false, nil } base1Value := locusValueObject.Base1Value base2Value := locusValueObject.Base2Value locusIsPhased := locusValueObject.LocusIsPhased return true, base1Value, base2Value, locusIsPhased, nil } //Outputs: // -map[string]string: Monogenic disease analysis map (contains probabilities for each genomeIdentifier and number of tested variants) // -[]map[string]string: Variants map list (contains variant info for each genomeIdentifier) // -error func getMonogenicDiseaseAnalysis(inputGenomesWithMetadataList []prepareRawGenomes.GenomeWithMetadata, diseaseObject monogenicDiseases.MonogenicDisease)(map[string]string, []map[string]string, error){ diseaseName := diseaseObject.DiseaseName dominantOrRecessive := diseaseObject.DominantOrRecessive variantsList := diseaseObject.VariantsList // We use this map to keep track of which RSIDs corresponds to each variant // Map Structure: Variant Identifier -> []rsID variantRSIDsMap := make(map[string][]int64) variantsMapList := make([]map[string]string, 0, len(variantsList)) for _, variantObject := range variantsList{ variantIdentifier := variantObject.VariantIdentifier // Map Structure: // -GenomeIdentifier -> ("Yes"/"No" + ";" + "Yes"/"No") or "Unknown" variantMap := map[string]string{ "ItemType": "MonogenicDiseaseVariant", "VariantIdentifier": variantIdentifier, "DiseaseName": diseaseName, } variantRSID := variantObject.VariantRSID variantDefectiveBase := variantObject.DefectiveBase variantRSIDsList := []int64{variantRSID} // We add aliases to variantRSIDsList anyAliasesExist, rsidAliasesList, err := locusMetadata.GetRSIDAliases(variantRSID) if (err != nil) { return nil, nil, err } if (anyAliasesExist == true){ variantRSIDsList = append(variantRSIDsList, rsidAliasesList...) } variantRSIDsMap[variantIdentifier] = variantRSIDsList for _, genomeWithMetadataObject := range inputGenomesWithMetadataList{ genomeIdentifier := genomeWithMetadataObject.GenomeIdentifier genomeMap := genomeWithMetadataObject.GenomeMap basePairValueFound, base1Value, base2Value, locusIsPhased, err := getGenomeLocusBasePair(genomeMap, variantRSID) if (err != nil) { return nil, nil, err } if (basePairValueFound == false){ variantMap[genomeIdentifier + "_HasVariant"] = "Unknown" continue } getBaseIsVariantMutationBool := func(inputBase string)bool{ if (inputBase == variantDefectiveBase){ return true } // Base could be mutated to a different unhealthy base // That mutation could be a neutral/healthier change // We only care about this specific variant return false } base1IsDefective := getBaseIsVariantMutationBool(base1Value) base2IsDefective := getBaseIsVariantMutationBool(base2Value) base1IsDefectiveString := helpers.ConvertBoolToYesOrNoString(base1IsDefective) base2IsDefectiveString := helpers.ConvertBoolToYesOrNoString(base2IsDefective) genomeHasVariantMapValue := base1IsDefectiveString + ";" + base2IsDefectiveString locusIsPhasedString := helpers.ConvertBoolToYesOrNoString(locusIsPhased) variantMap[genomeIdentifier + "_HasVariant"] = genomeHasVariantMapValue variantMap[genomeIdentifier + "_LocusIsPhased"] = locusIsPhasedString //TODO: Add LocusIsPhased to readGeneticAnalysis package } variantsMapList = append(variantsMapList, variantMap) } // Now we determine probability that user will pass disease to offspring, and probability that user has disease // We compute this for each genome diseaseAnalysisMap := make(map[string]string) // We will use this list when checking for conflicts allGenomeIdentifiersList := make([]string, 0, len(inputGenomesWithMetadataList)) for _, genomeWithMetadataObject := range inputGenomesWithMetadataList{ genomeIdentifier := genomeWithMetadataObject.GenomeIdentifier allGenomeIdentifiersList = append(allGenomeIdentifiersList, genomeIdentifier) numberOfVariantsTested := 0 for _, variantMap := range variantsMapList{ genomeHasVariantValue, exists := variantMap[genomeIdentifier + "_HasVariant"] if (exists == false){ return nil, nil, errors.New("variantMap malformed: Map missing a _HasVariant for genome") } if (genomeHasVariantValue != "Unknown"){ numberOfVariantsTested += 1 } } numberOfVariantsTestedString := helpers.ConvertIntToString(numberOfVariantsTested) diseaseAnalysisMap[genomeIdentifier + "_NumberOfVariantsTested"] = numberOfVariantsTestedString // Outputs: // -bool: Probability is known (will be false if the genome has no tested variants) // -float64: Probability Person has disease // -float64: Probability Person will pass a defect (variant) to offspring // -error getPersonDiseaseInfo := func()(bool, float64, float64, error){ anyProbabilityKnown := false // These variables are used to count the number of defective variants that exist on each chromosome numberOfVariants_Chromosome1 := 0 numberOfVariants_Chromosome2 := 0 numberOfVariants_UnknownChromosome := 0 // We use this map to keep track of how many mutations exist for each rsid // This allows us to know if 2 different variant mutations exist for a single RSID // For example, base1 is a different deleterious mutation than base2 // If this ever happens, we know that the user has the disease, // because both copies of the gene locus are defective. rsidMutationsMap := make(map[int64]int) for _, variantMap := range variantsMapList{ personHasVariantStatus, exists := variantMap[genomeIdentifier + "_HasVariant"] if (exists == false){ return false, 0, 0, errors.New("variantMap malformed: Map missing a _HasVariant for genome.") } if (personHasVariantStatus == "Unknown"){ // We don't know if the genome has the variant. Skip to next variant continue } anyProbabilityKnown = true locusIsPhasedStatus, exists := variantMap[genomeIdentifier + "_LocusIsPhased"] if (exists == false){ return false, 0, 0, errors.New("variantMap malformed: Map missing _LocusIsPhased for genome.") } base1HasVariant, base2HasVariant, foundSemicolon := strings.Cut(personHasVariantStatus, ";") if (foundSemicolon == false){ return false, 0, 0, errors.New("variantMap is malformed: Contains invalid personHasVariantValue: " + personHasVariantStatus) } if (base1HasVariant == "No" && base2HasVariant == "No"){ // Neither chromosome contains the variant mutation. continue } if (base1HasVariant == "Yes" && base2HasVariant == "Yes"){ // Both chromosomes contain the same variant mutation. // Person has the disease. // Person will definitely pass disease variant to offspring. return true, 1, 1, nil } // We know that this variant exists on 1 of the bases, but not both. variantIdentifier, exists := variantMap["VariantIdentifier"] if (exists == false){ return false, 0, 0, errors.New("VariantMap missing VariantIdentifier.") } variantRSIDsList, exists := variantRSIDsMap[variantIdentifier] if (exists == false){ return false, 0, 0, errors.New("variantRSIDMap missing variantIdentifier.") } for _, rsid := range variantRSIDsList{ rsidMutationsMap[rsid] += 1 } if (locusIsPhasedStatus == "Yes"){ if (base1HasVariant == "Yes"){ numberOfVariants_Chromosome1 += 1 } if (base2HasVariant == "Yes"){ numberOfVariants_Chromosome2 += 1 } } else { if (base1HasVariant == "Yes" || base2HasVariant == "Yes"){ numberOfVariants_UnknownChromosome += 1 } } } if (anyProbabilityKnown == false){ return false, 0, 0, nil } totalNumberOfVariants := numberOfVariants_Chromosome1 + numberOfVariants_Chromosome2 + numberOfVariants_UnknownChromosome if (totalNumberOfVariants == 0){ // Person does not have any disease variants. // They do not have the disease, and have no chance of passing a disease variant return true, 0, 0, nil } // Now we check to see if there are any loci which have 2 different variants, one for each base for _, numberOfMutations := range rsidMutationsMap{ if (numberOfMutations >= 2){ // Person has 2 mutations on the same location // They must have the disease, and will definitely pass a variant to their offspring return true, 1, 1, nil } } // At this point, we know that there are no homozygous variant mutations // All variant mutations are heterozygous, meaning the other chromosome base is healthy // Probability is expressed as a float between 0 - 1 getProbabilityPersonHasDisease := func()float64{ if (dominantOrRecessive == "Dominant"){ // Only 1 variant is needed for the person to have the disease // We know they have at least 1 variant return 1 } // dominantOrRecessive == "Recessive" if (totalNumberOfVariants == 1){ // There is only 1 variant in total. // This single variant cannot exist on both chromosomes. // The person does not have the disease return 0 } if (numberOfVariants_Chromosome1 >= 1 && numberOfVariants_Chromosome2 >= 1){ // We know there is at least 1 variant mutation on each chromosome // Therefore, the person has the disease return 1 } if (numberOfVariants_UnknownChromosome == 0){ // We know that variants do not exist on both chromosomes, only on 1. // Thus, the person does not have the disease return 0 } if (numberOfVariants_Chromosome1 == 0 && numberOfVariants_Chromosome2 == 0){ // All variants have an unknown phase, and we know there are multiple of them. // The probability the person does not have the disease is the probability that all mutations are on the same chromosome // We calculate the probability that all of the mutations are on the same chromosome // If they are all on the same chromosome, the person does not have the disease // If at least 1 variant exists on both chromosomes, the person has the disease // We calculate the probability that all variants are all on the same chromosome // Probability of n variants existing on the same chromosome = 1/(2^n) // P(X) = Probability all variants existing on chromosome 1 = 1/(2^n) // P(Y) = Probability all variants existing on chromosome 2 = 1/(2^n) // P(A U B) = P(A) + P(B) (If A and B are mutually exclusive) // P(X) and P(Y) are mutually exclusive // Probability of n variants existing on either chromosome 1 or 2 = P(X) + P(Y) // Probability person has disease = !P(X U Y) probabilityAllVariantsAreOnOneChromosome := 1/(math.Pow(2, float64(numberOfVariants_UnknownChromosome))) probabilityPersonHasDisease := 1 - (probabilityAllVariantsAreOnOneChromosome * 2) return probabilityPersonHasDisease } // We know that there is at least 1 variant whose phase is known // We know that there are no variants whose phase is known which exist on both chromosomes // We know there are at least some variants whose phase is unknown // The probability that the person has the disease is // the probability that the unknown-phase variants exist on the same chromosome as the ones which we know exist do // This probability is 50% for each unknown-phase variant probabilityAllVariantsAreOnOneChromosome := 1/(math.Pow(2, float64(numberOfVariants_UnknownChromosome))) probabilityPersonHasDisease := 1 - probabilityAllVariantsAreOnOneChromosome return probabilityPersonHasDisease } probabilityPersonHasDisease := getProbabilityPersonHasDisease() // We know all variants are heterozygous // Probability person will not pass any of n variants to their offspring: 1/(2^n) // Probability person will pass at least 1 of n variants to their offspring: 1 - (1/(2^n)) probabilityPersonWillPassAnyVariant := 1 - (1/(math.Pow(2, float64(totalNumberOfVariants)))) return true, probabilityPersonHasDisease, probabilityPersonWillPassAnyVariant, nil } probabilityKnown, probabilityPersonHasDisease, probabilityPersonWillPassAnyVariant, err := getPersonDiseaseInfo() if (err != nil) { return nil, nil, err } if (probabilityKnown == false){ diseaseAnalysisMap[genomeIdentifier + "_ProbabilityOfHavingDisease"] = "Unknown" diseaseAnalysisMap[genomeIdentifier + "_ProbabilityOfPassingADiseaseVariant"] = "Unknown" continue } percentageProbabilityPersonHasDisease := probabilityPersonHasDisease * 100 percentageProbabilityPersonWillPassADiseaseVariant := probabilityPersonWillPassAnyVariant * 100 probabilityOfHavingDiseaseString := helpers.ConvertFloat64ToStringRounded(percentageProbabilityPersonHasDisease, 0) probabilityOfPassingADiseaseVariantString := helpers.ConvertFloat64ToStringRounded(percentageProbabilityPersonWillPassADiseaseVariant, 0) diseaseAnalysisMap[genomeIdentifier + "_ProbabilityOfHavingDisease"] = probabilityOfHavingDiseaseString diseaseAnalysisMap[genomeIdentifier + "_ProbabilityOfPassingADiseaseVariant"] = probabilityOfPassingADiseaseVariantString } if (len(allGenomeIdentifiersList) == 1){ // We do not need to check for conflicts // Nothing left to do. Analysis is complete. return diseaseAnalysisMap, variantsMapList, nil } // We check for conflicts getConflictExistsBool := func()(bool, error){ // We start with disease analysis map probabilityOfHavingDisease := "" probabilityOfPassingAVariant := "" for index, genomeIdentifier := range allGenomeIdentifiersList{ currentProbabilityOfHavingDisease, exists := diseaseAnalysisMap[genomeIdentifier + "_ProbabilityOfHavingDisease"] if (exists == false){ return false, errors.New("Cannot create analysis: diseaseAnalysisMap missing _ProbabilityOfHavingDisease") } currentProbabilityOfPassingAVariant, exists := diseaseAnalysisMap[genomeIdentifier + "_ProbabilityOfPassingADiseaseVariant"] if (exists == false){ return false, errors.New("Cannot create analysis: diseaseAnalysisMap missing _ProbabilityOfPassingADiseaseVariant") } if (index == 0){ probabilityOfHavingDisease = currentProbabilityOfHavingDisease probabilityOfPassingAVariant = currentProbabilityOfPassingAVariant continue } if (currentProbabilityOfHavingDisease != probabilityOfHavingDisease){ return true, nil } if (currentProbabilityOfPassingAVariant != probabilityOfPassingAVariant){ return true, nil } } // Now we test variants for conflicts for _, variantMap := range variantsMapList{ // Each variant value is either "Yes;No", "No;Yes", "No;No", or "Yes;Yes" // Two different genomes have "Yes;No" and "No;Yes", it will not count as a conflict // If the locus is unphased, then there is no difference between "Yes;No" and "No;Yes" // If the locus is phased, then this flip is only meaningful if it effects the probability of disease/passing a variant // We already checked those probabilities for conflicts earlier // Therefore, any flip is not considered a conflict genomeHasVariantStatus := "" for index, genomeIdentifier := range allGenomeIdentifiersList{ currentGenomeHasVariantStatus, exists := variantMap[genomeIdentifier + "_HasVariant"] if (exists == false){ return false, errors.New("variantMap missing genomeIdentifier key when checking for conflicts") } if (index == 0){ genomeHasVariantStatus = currentGenomeHasVariantStatus continue } if (currentGenomeHasVariantStatus != genomeHasVariantStatus){ if (currentGenomeHasVariantStatus == "Yes;No" && genomeHasVariantStatus == "No;Yes"){ continue } if (currentGenomeHasVariantStatus == "No;Yes" && genomeHasVariantStatus == "Yes;No"){ continue } return true, nil } } } return false, nil } conflictExists, err := getConflictExistsBool() if (err != nil) { return nil, nil, err } conflictExistsString := helpers.ConvertBoolToYesOrNoString(conflictExists) diseaseAnalysisMap["ConflictExists"] = conflictExistsString return diseaseAnalysisMap, variantsMapList, nil } //Outputs: // -map[string]string: Polygenic Disease analysis map (contains weight/probability for each genomeIdentifier and number of loci tested) // -[]map[string]string: Loci map list (contains locus info for each genomeIdentifier) // -error func getPolygenicDiseaseAnalysis(inputGenomesWithMetadataList []prepareRawGenomes.GenomeWithMetadata, diseaseObject polygenicDiseases.PolygenicDisease)(map[string]string, []map[string]string, error){ diseaseName := diseaseObject.DiseaseName diseaseLociList := diseaseObject.LociList lociMapList := make([]map[string]string, 0, len(diseaseLociList)) // Map Structure: Locus identifier -> Disease Locus object locusObjectsMap := make(map[string]polygenicDiseases.DiseaseLocus) for _, locusObject := range diseaseLociList{ locusIdentifier := locusObject.LocusIdentifier locusObjectsMap[locusIdentifier] = locusObject locusMap := map[string]string{ "ItemType": "PolygenicDiseaseLocus", "LocusIdentifier": locusIdentifier, "DiseaseName": diseaseName, } locusRSID := locusObject.LocusRSID locusRiskWeightsMap := locusObject.RiskWeightsMap locusOddsRatiosMap := locusObject.OddsRatiosMap for _, genomeWithMetadataObject := range inputGenomesWithMetadataList{ genomeIdentifier := genomeWithMetadataObject.GenomeIdentifier genomeMap := genomeWithMetadataObject.GenomeMap //Outputs: // -bool: Disease locus info is known // -string: Locus Bases // -int: Genome disease locus risk weight // -bool: Genome disease locus odds ratio known // -float64: Genome disease locus odds ratio // -error getGenomeDiseaseLocusInfo := func()(bool, string, int, bool, float64, error){ basePairValueFound, base1Value, base2Value, _, err := getGenomeLocusBasePair(genomeMap, locusRSID) if (err != nil) { return false, "", 0, false, 0, err } if (basePairValueFound == false){ return false, "", 0, false, 0, nil } locusBasePair := base1Value + ";" + base2Value riskWeight, exists := locusRiskWeightsMap[locusBasePair] if (exists == false){ // This is an unknown base combination // We will treat it as a 0 risk weight return true, locusBasePair, 0, true, 1, nil } if (riskWeight == 0){ return true, locusBasePair, 0, true, 1, nil } oddsRatio, exists := locusOddsRatiosMap[locusBasePair] if (exists == false){ return true, locusBasePair, 0, false, 0, nil } return true, locusBasePair, riskWeight, true, oddsRatio, nil } diseaseLocusInfoKnown, locusBasePair, locusRiskWeight, locusOddsRatioKnown, locusOddsRatio, err := getGenomeDiseaseLocusInfo() if (err != nil) { return nil, nil, err } if (diseaseLocusInfoKnown == false){ locusMap[genomeIdentifier + "_LocusBasePair"] = "Unknown" locusMap[genomeIdentifier + "_RiskWeight"] = "Unknown" locusMap[genomeIdentifier + "_OddsRatio"] = "Unknown" } else { locusMap[genomeIdentifier + "_LocusBasePair"] = locusBasePair riskWeightString := helpers.ConvertIntToString(locusRiskWeight) locusMap[genomeIdentifier + "_RiskWeight"] = riskWeightString if (locusOddsRatioKnown == false){ locusMap[genomeIdentifier + "_OddsRatio"] = "Unknown" } else { locusOddsRatioString := helpers.ConvertFloat64ToString(locusOddsRatio) locusMap[genomeIdentifier + "_OddsRatio"] = locusOddsRatioString } } } lociMapList = append(lociMapList, locusMap) } // Now we construct polygenic disease probability map for each genome diseaseAnalysisMap := make(map[string]string) // We use this list to check for conflicts later allGenomeIdentifiersList := make([]string, 0, len(inputGenomesWithMetadataList)) for _, genomeWithMetadataObject := range inputGenomesWithMetadataList{ genomeIdentifier := genomeWithMetadataObject.GenomeIdentifier allGenomeIdentifiersList = append(allGenomeIdentifiersList, genomeIdentifier) numberOfLociTested := 0 minimumPossibleRiskWeightSum := 0 maximumPossibleRiskWeightSum := 0 summedDiseaseRiskWeight := 0 for _, locusMap := range lociMapList{ genomeRiskWeight, exists := locusMap[genomeIdentifier + "_RiskWeight"] if (exists == false){ return nil, nil, errors.New("locusMap malformed: Map missing a genomeIdentifier riskWeight") } if (genomeRiskWeight == "Unknown"){ // The genome does not have this locus // We cannot test for risk weight continue } numberOfLociTested += 1 locusIdentifier, exists := locusMap["LocusIdentifier"] if (exists == false) { return nil, nil, errors.New("locusMap malformed: Map missing LocusIdentifier") } locusObject, exists := locusObjectsMap[locusIdentifier] if (exists == false){ return nil, nil, errors.New("LocusObjectsMap missing locus: " + locusIdentifier) } locusMinimumWeight := locusObject.MinimumRiskWeight locusMaximumWeight := locusObject.MaximumRiskWeight minimumPossibleRiskWeightSum += locusMinimumWeight maximumPossibleRiskWeightSum += locusMaximumWeight genomeRiskWeightInt, err := helpers.ConvertStringToInt(genomeRiskWeight) if (err != nil) { return nil, nil, errors.New("Locus map malformed: Contains invalid _RiskWeight: " + genomeRiskWeight) } summedDiseaseRiskWeight += genomeRiskWeightInt } numberOfLociTestedString := helpers.ConvertIntToString(numberOfLociTested) diseaseAnalysisMap[genomeIdentifier + "_NumberOfLociTested"] = numberOfLociTestedString if (numberOfLociTested == 0){ diseaseAnalysisMap[genomeIdentifier + "_RiskScore"] = "Unknown" continue } diseaseRiskScore, err := helpers.ScaleNumberProportionally(true, summedDiseaseRiskWeight, minimumPossibleRiskWeightSum, maximumPossibleRiskWeightSum, 0, 10) if (err != nil) { return nil, nil, err } diseaseRiskScoreString := helpers.ConvertIntToString(diseaseRiskScore) diseaseAnalysisMap[genomeIdentifier + "_RiskScore"] = diseaseRiskScoreString } if (len(allGenomeIdentifiersList) == 1){ // We do not need to check for conflicts // Nothing left to do. Analysis is complete. return diseaseAnalysisMap, lociMapList, nil } // We check for conflicts getConflictExistsBool := func()(bool, error){ for _, locusMap := range lociMapList{ locusBasePair := "" for index, genomeIdentifier := range allGenomeIdentifiersList{ currentLocusBasePair, exists := locusMap[genomeIdentifier + "_LocusBasePair"] if (exists == false){ return false, errors.New("Cannot find _LocusBasePair key when searching for conflicts") } if (index == 0){ locusBasePair = currentLocusBasePair continue } if (currentLocusBasePair != locusBasePair){ return true, nil } } } return false, nil } conflictExists, err := getConflictExistsBool() if (err != nil) { return nil, nil, err } conflictExistsString := helpers.ConvertBoolToYesOrNoString(conflictExists) diseaseAnalysisMap["ConflictExists"] = conflictExistsString return diseaseAnalysisMap, lociMapList, nil } //Outputs: // -map[string]string: Trait analysis map (contains weight/probability for each genomeIdentifier and number of rules tested) // -[]map[string]string: Loci map list (contains locus info for each genomeIdentifier) // -error func getTraitAnalysis(inputGenomesWithMetadataList []prepareRawGenomes.GenomeWithMetadata, traitObject traits.Trait)(map[string]string, []map[string]string, error){ traitName := traitObject.TraitName traitLociList := traitObject.LociList traitRulesList := traitObject.RulesList // We first add loci values to trait analysis map traitAnalysisMap := make(map[string]string) // We use this list to check for conflicts later allGenomeIdentifiersList := make([]string, 0, len(inputGenomesWithMetadataList)) for _, genomeWithMetadataObject := range inputGenomesWithMetadataList{ genomeIdentifier := genomeWithMetadataObject.GenomeIdentifier allGenomeIdentifiersList = append(allGenomeIdentifiersList, genomeIdentifier) genomeMap := genomeWithMetadataObject.GenomeMap for _, locusRSID := range traitLociList{ locusRSIDString := helpers.ConvertInt64ToString(locusRSID) locusBasePairKnown, locusBase1, locusBase2, locusIsPhased, err := getGenomeLocusBasePair(genomeMap, locusRSID) if (err != nil) { return nil, nil, err } if (locusBasePairKnown == false){ traitAnalysisMap[genomeIdentifier + "_LocusValue_rs" + locusRSIDString] = "Unknown" } else { locusBasePair := locusBase1 + ";" + locusBase2 locusIsPhasedString := helpers.ConvertBoolToYesOrNoString(locusIsPhased) traitAnalysisMap[genomeIdentifier + "_LocusValue_rs" + locusRSIDString] = locusBasePair traitAnalysisMap[genomeIdentifier + "_LocusIsPhased_rs" + locusRSIDString] = locusIsPhasedString } } } // This map stores the TraitRule maps rulesMapList := make([]map[string]string, 0, len(traitRulesList)) // Map Structure: Rule identifier -> Rule object ruleObjectsMap := make(map[string]traits.TraitRule) if (len(traitRulesList) != 0){ // This trait contains at least 1 rule for _, ruleObject := range traitRulesList{ ruleIdentifier := ruleObject.RuleIdentifier ruleObjectsMap[ruleIdentifier] = ruleObject ruleMap := map[string]string{ "ItemType": "TraitRule", "RuleIdentifier": ruleIdentifier, "TraitName": traitName, } ruleLociList := ruleObject.LociList for _, genomeWithMetadataObject := range inputGenomesWithMetadataList{ genomeIdentifier := genomeWithMetadataObject.GenomeIdentifier genomeMap := genomeWithMetadataObject.GenomeMap // We add locus base pairs to ruleMap // We also check to see if genome passes all rule loci // We only consider a rule Known if the genome either passes all loci, or fails to pass 1 locus allRuleLociKnown := true anyLocusFailureKnown := false //This is true if at least 1 locus is known to not pass for _, locusObject := range ruleLociList{ locusIdentifier := locusObject.LocusIdentifier locusRSID := locusObject.LocusRSID locusBasePairKnown, locusBase1, locusBase2, _, err := getGenomeLocusBasePair(genomeMap, locusRSID) if (err != nil) { return nil, nil, err } if (locusBasePairKnown == false){ ruleMap[genomeIdentifier + "_RuleLocusBasePair_" + locusIdentifier] = "Unknown" // The genome has failed to pass a single rule locus, thus, the rule is not passed // We still continue so we can add all locus base pairs to the ruleMap allRuleLociKnown = false continue } locusBasePair := locusBase1 + ";" + locusBase2 ruleMap[genomeIdentifier + "_RuleLocusBasePair_" + locusIdentifier] = locusBasePair if (anyLocusFailureKnown == true){ // We don't have to check if the rule locus is passed, because we already know the rule is not passed continue } locusBasePairsList := locusObject.BasePairsList genomePassesRuleLocus := slices.Contains(locusBasePairsList, locusBasePair) if (genomePassesRuleLocus == false){ // We know the rule is not passed // We still continue so we can add all locus base pairs to the ruleMap anyLocusFailureKnown = true continue } } if (anyLocusFailureKnown == true){ ruleMap[genomeIdentifier + "_PassesRule"] = "No" } else if (allRuleLociKnown == false){ ruleMap[genomeIdentifier + "_PassesRule"] = "Unknown" } else { ruleMap[genomeIdentifier + "_PassesRule"] = "Yes" } } rulesMapList = append(rulesMapList, ruleMap) } // Now we construct trait outcome scores map for each genome traitOutcomesList := traitObject.OutcomesList for _, genomeWithMetadataObject := range inputGenomesWithMetadataList{ genomeIdentifier := genomeWithMetadataObject.GenomeIdentifier // A rule is considered tested if at least 1 locus within the rule is known to fail, or all loci within the rule pass numberOfRulesTested := 0 // Map Structure: Trait outcome -> Number of points traitOutcomeScoresMap := make(map[string]int) for _, ruleMap := range rulesMapList{ genomePassesRule, exists := ruleMap[genomeIdentifier + "_PassesRule"] if (exists == false){ return nil, nil, errors.New("ruleMap malformed: Map missing a genomeIdentifier _PassesRule") } if (genomePassesRule == "Unknown"){ continue } numberOfRulesTested += 1 if (genomePassesRule == "No"){ continue } ruleIdentifier, exists := ruleMap["RuleIdentifier"] if (exists == false) { return nil, nil, errors.New("ruleMap malformed: Map missing RuleIdentifier") } ruleObject, exists := ruleObjectsMap[ruleIdentifier] if (exists == false){ return nil, nil, errors.New("RuleObjectsMap missing rule: " + ruleIdentifier) } ruleOutcomePointsMap := ruleObject.OutcomePointsMap for traitOutcome, pointsChange := range ruleOutcomePointsMap{ traitOutcomeScoresMap[traitOutcome] += pointsChange } } numberOfRulesTestedString := helpers.ConvertIntToString(numberOfRulesTested) traitAnalysisMap[genomeIdentifier + "_NumberOfRulesTested"] = numberOfRulesTestedString if (numberOfRulesTested == 0){ // No rules were tested. The number of points in each outcome are all unknown. continue } // We add all outcomes for which there were no points for _, traitOutcome := range traitOutcomesList{ _, exists := traitOutcomeScoresMap[traitOutcome] if (exists == false){ traitOutcomeScoresMap[traitOutcome] = 0 } } for traitOutcome, outcomeScore := range traitOutcomeScoresMap{ outcomeScoreString := helpers.ConvertIntToString(outcomeScore) traitAnalysisMap[genomeIdentifier + "_OutcomeScore_" + traitOutcome] = outcomeScoreString } } } else { // No rules exist for this trait for _, genomeIdentifier := range allGenomeIdentifiersList{ traitAnalysisMap[genomeIdentifier + "_NumberOfRulesTested"] = "0" } } if (len(allGenomeIdentifiersList) == 1){ // We do not need to check for conflicts // Nothing left to do. Analysis is complete. return traitAnalysisMap, rulesMapList, nil } // We check for conflicts getConflictExistsBool := func()(bool, error){ //TODO: Check for locus value conflicts once locus values are used in neural network prediction. // We only have to check each rule result, because they will determine the overall person results for _, ruleMap := range rulesMapList{ passesRule := "" for index, genomeIdentifier := range allGenomeIdentifiersList{ currentPassesRule, exists := ruleMap[genomeIdentifier + "_PassesRule"] if (exists == false){ return false, errors.New("Cannot find _PassesRule key when searching for conflicts") } if (index == 0){ passesRule = currentPassesRule continue } if (currentPassesRule != passesRule){ return true, nil } } } return false, nil } conflictExists, err := getConflictExistsBool() if (err != nil) { return nil, nil, err } conflictExistsString := helpers.ConvertBoolToYesOrNoString(conflictExists) traitAnalysisMap["ConflictExists"] = conflictExistsString return traitAnalysisMap, rulesMapList, nil }