// createCoupleGeneticAnalysis provides functions to create a Couple genetic analysis // Couple analyses provide an analysis of the prospective offspring of a pair of people // These analyses are performed on one or more genome files. // Analyses contain 3 categories of results: Monogenic Diseases, Polygenic Diseases and Traits // Use createPersonGeneticAnalysis.go to create Person analyses package createCoupleGeneticAnalysis // Disclaimer: I am a novice in the ways of genetics. This package could be flawed in numerous ways. // TODO: We want to eventually use neural nets for polygenic disease analysis (see geneticPrediction.go) // This is only possible once we get access to the necessary training data import "seekia/resources/trainedPredictionModels" import "seekia/resources/geneticReferences/locusMetadata" import "seekia/resources/geneticReferences/monogenicDiseases" import "seekia/resources/geneticReferences/polygenicDiseases" import "seekia/resources/geneticReferences/traits" import "seekia/internal/encoding" import "seekia/internal/genetics/createPersonGeneticAnalysis" import "seekia/internal/genetics/geneticAnalysis" import "seekia/internal/genetics/locusValue" import "seekia/internal/genetics/prepareRawGenomes" import "seekia/internal/helpers" import "errors" import mathRand "math/rand/v2" import "maps" import "reflect" //Outputs: // -bool: Process completed (was not manually stopped mid-way) // -string: Couple genetic analysis string (encoded in MessagePack) // -error func CreateCoupleGeneticAnalysis(person1GenomesList []prepareRawGenomes.RawGenomeWithMetadata, person2GenomesList []prepareRawGenomes.RawGenomeWithMetadata, updatePercentageCompleteFunction func(int)error, checkIfProcessIsStopped func()bool)(bool, string, error){ person1PrepareRawGenomesUpdatePercentageCompleteFunction := func(newPercentage int)error{ newPercentageCompletion, err := helpers.ScaleIntProportionally(true, newPercentage, 0, 100, 0, 25) if (err != nil) { return err } err = updatePercentageCompleteFunction(newPercentageCompletion) if (err != nil) { return err } return nil } anyUsefulLocationsExist, person1GenomesWithMetadataList, allPerson1RawGenomeIdentifiersList, person1HasMultipleGenomes, person1OnlyExcludeConflictsGenomeIdentifier, person1OnlyIncludeSharedGenomeIdentifier, err := prepareRawGenomes.GetGenomesWithMetadataListFromRawGenomesList(person1GenomesList, person1PrepareRawGenomesUpdatePercentageCompleteFunction) if (err != nil) { return false, "", err } if (anyUsefulLocationsExist == false){ // We should have checked for this when genomes were first imported. return false, "", errors.New("CreateCoupleGeneticAnalysis called with person1GenomesList that does not contain any useful genomes") } if (len(person1GenomesList) > 1 && (len(person1GenomesList) != (len(person1GenomesWithMetadataList)-2)) ){ // If there is more than 1 genome, 2 combined genomes are created // We are checking to make sure that none of the input genomes were dropped due to not having any locations // We should have checked to make sure each input genome has useful locations when each genome was first imported. return false, "", errors.New("CreateCoupleGeneticAnalysis called with person1GenomesList containing at least 1 genome without useful locations.") } processIsStopped := checkIfProcessIsStopped() if (processIsStopped == true){ return false, "", nil } person2PrepareRawGenomesUpdatePercentageCompleteFunction := func(newPercentage int)error{ newPercentageCompletion, err := helpers.ScaleIntProportionally(true, newPercentage, 0, 100, 25, 50) if (err != nil){ return err } err = updatePercentageCompleteFunction(newPercentageCompletion) if (err != nil) { return err } return nil } anyUsefulLocationsExist, person2GenomesWithMetadataList, allPerson2RawGenomeIdentifiersList, person2HasMultipleGenomes, person2OnlyExcludeConflictsGenomeIdentifier, person2OnlyIncludeSharedGenomeIdentifier, err := prepareRawGenomes.GetGenomesWithMetadataListFromRawGenomesList(person2GenomesList, person2PrepareRawGenomesUpdatePercentageCompleteFunction) if (err != nil) { return false, "", err } if (anyUsefulLocationsExist == false){ // We should have checked for this when genomes were first imported. return false, "", errors.New("CreateCoupleGeneticAnalysis called with person2GenomesList that does not contain any useful genomes") } if (len(person2GenomesList) > 1 && (len(person2GenomesList) != (len(person2GenomesWithMetadataList)-2)) ){ // If there is more than 1 genome, 2 combined genomes are created // We are checking to make sure that none of the input genomes were dropped due to not having any locations // We should have checked to make sure each input genome has useful locations when each genome was first imported. return false, "", errors.New("CreateCoupleGeneticAnalysis called with person2GenomesList containing at least 1 genome without useful locations.") } processIsStopped = checkIfProcessIsStopped() if (processIsStopped == true){ return false, "", nil } // This map stores each genome's locus values // Map Structure: Genome Identifier -> Genome locus values map (rsID -> Locus Value) person1GenomesMap := make(map[[16]byte]map[int64]locusValue.LocusValue) for _, genomeWithMetadata := range person1GenomesWithMetadataList{ genomeIdentifier := genomeWithMetadata.GenomeIdentifier genomeMap := genomeWithMetadata.GenomeMap person1GenomesMap[genomeIdentifier] = genomeMap } // This map stores each genome's locus values // Map Structure: Genome Identifier -> Genome locus values map (rsID -> Locus Value) person2GenomesMap := make(map[[16]byte]map[int64]locusValue.LocusValue) for _, genomeWithMetadata := range person2GenomesWithMetadataList{ genomeIdentifier := genomeWithMetadata.GenomeIdentifier genomeMap := genomeWithMetadata.GenomeMap person2GenomesMap[genomeIdentifier] = genomeMap } // The analysis will analyze either 1 or 2 genome pairs // The gui will display the results from each pair //Outputs: // -[16]byte: Pair 1 Person 1 Genome Identifier // -[16]byte: Pair 1 Person 2 Genome Identifier // -bool: Second pair exists // -[16]byte: Pair 2 Person 1 Genome Identifier // -[16]byte: Pair 2 Person 2 Genome Identifier // -error getGenomePairsToAnalyze := func()([16]byte, [16]byte, bool, [16]byte, [16]byte, error){ if (person1HasMultipleGenomes == true && person2HasMultipleGenomes == true){ return person1OnlyExcludeConflictsGenomeIdentifier, person2OnlyExcludeConflictsGenomeIdentifier, true, person1OnlyIncludeSharedGenomeIdentifier, person2OnlyIncludeSharedGenomeIdentifier, nil } if (person1HasMultipleGenomes == true && person2HasMultipleGenomes == false){ person2GenomeIdentifier := allPerson2RawGenomeIdentifiersList[0] return person1OnlyExcludeConflictsGenomeIdentifier, person2GenomeIdentifier, true, person1OnlyIncludeSharedGenomeIdentifier, person2GenomeIdentifier, nil } if (person1HasMultipleGenomes == false && person2HasMultipleGenomes == true){ person1GenomeIdentifier := allPerson1RawGenomeIdentifiersList[0] return person1GenomeIdentifier, person2OnlyExcludeConflictsGenomeIdentifier, true, person1GenomeIdentifier, person2OnlyIncludeSharedGenomeIdentifier, nil } //person1HasMultipleGenomes == false && person2HasMultipleGenomes == false person1GenomeIdentifier := allPerson1RawGenomeIdentifiersList[0] person2GenomeIdentifier := allPerson2RawGenomeIdentifiersList[0] return person1GenomeIdentifier, person2GenomeIdentifier, false, [16]byte{}, [16]byte{}, nil } pair1Person1GenomeIdentifier, pair1Person2GenomeIdentifier, genomePair2Exists, pair2Person1GenomeIdentifier, pair2Person2GenomeIdentifier, err := getGenomePairsToAnalyze() if (err != nil){ return false, "", err } newCoupleAnalysis := geneticAnalysis.CoupleAnalysis{ AnalysisVersion: 1, Pair1Person1GenomeIdentifier: pair1Person1GenomeIdentifier, Pair1Person2GenomeIdentifier: pair1Person2GenomeIdentifier, SecondPairExists: genomePair2Exists, Person1HasMultipleGenomes: person1HasMultipleGenomes, Person2HasMultipleGenomes: person2HasMultipleGenomes, } if (genomePair2Exists == true){ // At least 1 of the people have multiple genomes newCoupleAnalysis.Pair2Person1GenomeIdentifier = pair2Person1GenomeIdentifier newCoupleAnalysis.Pair2Person2GenomeIdentifier = pair2Person2GenomeIdentifier if (person1HasMultipleGenomes == true){ newCoupleAnalysis.Person1OnlyExcludeConflictsGenomeIdentifier = person1OnlyExcludeConflictsGenomeIdentifier newCoupleAnalysis.Person1OnlyIncludeSharedGenomeIdentifier = person1OnlyIncludeSharedGenomeIdentifier } if (person2HasMultipleGenomes == true){ newCoupleAnalysis.Person2OnlyExcludeConflictsGenomeIdentifier = person2OnlyExcludeConflictsGenomeIdentifier newCoupleAnalysis.Person2OnlyIncludeSharedGenomeIdentifier = person2OnlyIncludeSharedGenomeIdentifier } } // We compute and add monogenic disease probabilities monogenicDiseasesList, err := monogenicDiseases.GetMonogenicDiseaseObjectsList() if (err != nil) { return false, "", err } // Map Structure: Disease Name -> OffspringMonogenicDiseaseInfo offspringMonogenicDiseasesMap := make(map[string]geneticAnalysis.OffspringMonogenicDiseaseInfo) for _, diseaseObject := range monogenicDiseasesList{ diseaseName := diseaseObject.DiseaseName variantsList := diseaseObject.VariantsList diseaseIsDominantOrRecessive := diseaseObject.DominantOrRecessive person1DiseaseAnalysisObject, err := createPersonGeneticAnalysis.GetPersonMonogenicDiseaseAnalysis(person1GenomesWithMetadataList, diseaseObject) if (err != nil) { return false, "", err } person2DiseaseAnalysisObject, err := createPersonGeneticAnalysis.GetPersonMonogenicDiseaseAnalysis(person2GenomesWithMetadataList, diseaseObject) if (err != nil) { return false, "", err } // This map stores the quantity of variants tested in each person's genome // Map Structure: Genome Identifier -> Number of variants tested quantityOfVariantsTestedMap := make(map[[16]byte]int) // This map stores the offspring disease probabilities for each genome pair. // A genome pair is a concatenation of two genome identifiers // If a map entry doesn't exist, the probabilities are unknown for that genome pair // Map Structure: Genome Pair Identifier -> OffspringMonogenicDiseaseProbabilities offspringMonogenicDiseaseInfoMap := make(map[[32]byte]geneticAnalysis.OffspringGenomePairMonogenicDiseaseInfo) // This will calculate the probability of monogenic disease for the offspring from the two specified genomes // It also calculates the probabilities for each monogenic disease variant for the offspring // It then adds the genome pair disease information to the offspringMonogenicDiseaseInfoMap and quantityOfVariantsTestedMap addGenomePairInfoToDiseaseMaps := func(person1GenomeIdentifier [16]byte, person2GenomeIdentifier [16]byte)error{ //Outputs: // -bool: Probability is known // -int: Probability of passing a disease variant (value between 0 and 100) // -int: Number of variants tested getPersonWillPassDiseaseVariantProbability := func(personDiseaseAnalysisObject geneticAnalysis.PersonMonogenicDiseaseInfo, genomeIdentifier [16]byte)(bool, int, int){ personDiseaseInfoMap := personDiseaseAnalysisObject.MonogenicDiseaseInfoMap personGenomeDiseaseInfoObject, exists := personDiseaseInfoMap[genomeIdentifier] if (exists == false){ return false, 0, 0 } personGenomeProbabilityOfPassingADiseaseVariant := personGenomeDiseaseInfoObject.ProbabilityOfPassingADiseaseVariant personGenomeQuantityOfVariantsTested := personGenomeDiseaseInfoObject.QuantityOfVariantsTested return true, personGenomeProbabilityOfPassingADiseaseVariant, personGenomeQuantityOfVariantsTested } person1ProbabilityIsKnown, person1WillPassVariantProbability, person1NumberOfVariantsTested := getPersonWillPassDiseaseVariantProbability(person1DiseaseAnalysisObject, person1GenomeIdentifier) person2ProbabilityIsKnown, person2WillPassVariantProbability, person2NumberOfVariantsTested := getPersonWillPassDiseaseVariantProbability(person2DiseaseAnalysisObject, person2GenomeIdentifier) offspringHasDiseaseProbabilityIsKnown, percentageProbabilityOffspringHasDisease, offspringHasAVariantProbabilityIsKnown, percentageProbabilityOffspringHasVariant, err := GetOffspringMonogenicDiseaseProbabilities(diseaseIsDominantOrRecessive, person1ProbabilityIsKnown, person1WillPassVariantProbability, person2ProbabilityIsKnown, person2WillPassVariantProbability) if (err != nil) { return err } // Now we calculate the probabilities for individual variants offspringVariantsInfoMap := make(map[[3]byte]geneticAnalysis.OffspringMonogenicDiseaseVariantInfo) for _, variantObject := range variantsList{ variantIdentifierHex := variantObject.VariantIdentifier variantIdentifier, err := encoding.DecodeHexStringTo3ByteArray(variantIdentifierHex) if (err != nil) { return err } //Outputs: // -bool: Probability is known // -float64: Probability that person will pass variant to offspring (between 0 and 1) // -error getPersonWillPassVariantProbability := func(personDiseaseAnalysisObject geneticAnalysis.PersonMonogenicDiseaseInfo, personGenomeIdentifier [16]byte)(bool, float64, error){ personDiseaseInfoMap := personDiseaseAnalysisObject.MonogenicDiseaseInfoMap personGenomeDiseaseInfoObject, exists := personDiseaseInfoMap[personGenomeIdentifier] if (exists == false){ return false, 0, nil } personGenomeDiseaseVariantsMap := personGenomeDiseaseInfoObject.VariantsInfoMap personVariantInfoObject, exists := personGenomeDiseaseVariantsMap[variantIdentifier] if (exists == false){ // The genome does not have information about this variant return false, 0, nil } base1HasVariant := personVariantInfoObject.Base1HasVariant base2HasVariant := personVariantInfoObject.Base2HasVariant if (base1HasVariant == true && base2HasVariant == true){ return true, 1, nil } if (base1HasVariant == true && base2HasVariant == false){ return true, 0.5, nil } if (base1HasVariant == false && base2HasVariant == true){ return true, 0.5, nil } // Neither base has a variant return true, 0, nil } person1VariantProbabilityIsKnown, person1WillPassVariantProbability, err := getPersonWillPassVariantProbability(person1DiseaseAnalysisObject, person1GenomeIdentifier) if (err != nil) { return err } person2VariantProbabilityIsKnown, person2WillPassVariantProbability, err := getPersonWillPassVariantProbability(person2DiseaseAnalysisObject, person2GenomeIdentifier) if (err != nil) { return err } if (person1VariantProbabilityIsKnown == false && person2VariantProbabilityIsKnown == false){ continue } // Outputs: // -float64: Best Case Person 1 will pass variant probability (0-1) // -float64: Worst Case Person 1 will pass variant probability (0-1) // -float64: Best Case Person 2 will pass variant probability (0-1) // -float64: Worst Case Person 2 will pass variant probability (0-1) getBestAndWorstCaseProbabilities := func()(float64, float64, float64, float64){ if (person1VariantProbabilityIsKnown == true && person2VariantProbabilityIsKnown == false){ return person1WillPassVariantProbability, person1WillPassVariantProbability, float64(0), float64(1) } if (person1VariantProbabilityIsKnown == false && person2VariantProbabilityIsKnown == true){ return float64(0), float64(1), person2WillPassVariantProbability, person2WillPassVariantProbability } // Both people's probabilities are known return person1WillPassVariantProbability, person1WillPassVariantProbability, person2WillPassVariantProbability, person2WillPassVariantProbability } bestCasePerson1WillPassVariantProbability, worstCasePerson1WillPassVariantProbability, bestCasePerson2WillPassVariantProbability, worstCasePerson2WillPassVariantProbability := getBestAndWorstCaseProbabilities() //Outputs: // -int: Percentage Probability of 0 mutations // -int: Percentage Probability of 1 mutation // -int: Percentage Probability of 2 mutations // -error getOffspringVariantMutationProbabilities := func(person1WillPassVariantProbability float64, person2WillPassVariantProbability float64)(int, int, int, error){ // This is the probability that neither person will pass the variant // P = P(!A) * P(!B) probabilityOf0Mutations := (1 - person1WillPassVariantProbability) * (1 - person2WillPassVariantProbability) // This is the probability that either person will pass the variant, but not both // P(A XOR B) = P(A) + P(B) - (2 * P(A and B)) probabilityOf1Mutation := person1WillPassVariantProbability + person2WillPassVariantProbability - (2 * person1WillPassVariantProbability * person2WillPassVariantProbability) // This is the probability that both people will pass the variant // P(A and B) = P(A) * P(B) probabilityOf2Mutations := person1WillPassVariantProbability * person2WillPassVariantProbability 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 } bestCase0MutationsProbability, bestCase1MutationProbability, bestCase2MutationsProbability, err := getOffspringVariantMutationProbabilities(bestCasePerson1WillPassVariantProbability, bestCasePerson2WillPassVariantProbability) if (err != nil) { return err } worstCase0MutationsProbability, worstCase1MutationProbability, worstCase2MutationsProbability, err := getOffspringVariantMutationProbabilities(worstCasePerson1WillPassVariantProbability, worstCasePerson2WillPassVariantProbability) if (err != nil) { return err } // We have to figure out which 1-mutation-probability is lower // The best case probabilities can actually result in a higher probability for 1 mutation // than the worst case probabilities // This is because in a worst-case-scenario, the probability of 1 mutation might be 0 because // the probability of 2 mutations is 100 // This is not needed for the 0 and 2 mutation probabilities because TODO lowerBound1MutationProbability := min(bestCase1MutationProbability, worstCase1MutationProbability) upperBound1MutationProbability := max(bestCase1MutationProbability, worstCase1MutationProbability) newOffspringMonogenicDiseaseVariantInfoObject := geneticAnalysis.OffspringMonogenicDiseaseVariantInfo{ ProbabilityOf0MutationsLowerBound: worstCase0MutationsProbability, ProbabilityOf0MutationsUpperBound: bestCase0MutationsProbability, ProbabilityOf1MutationLowerBound: lowerBound1MutationProbability, ProbabilityOf1MutationUpperBound: upperBound1MutationProbability, ProbabilityOf2MutationsLowerBound: bestCase2MutationsProbability, ProbabilityOf2MutationsUpperBound: worstCase2MutationsProbability, } offspringVariantsInfoMap[variantIdentifier] = newOffspringMonogenicDiseaseVariantInfoObject } if (person1ProbabilityIsKnown == false && person2ProbabilityIsKnown == false && len(offspringVariantsInfoMap) == 0){ // We have no information about this genome pair's disease risk // We don't add this genome pair's disease info to the diseaseInfoMap return nil } quantityOfVariantsTestedMap[person1GenomeIdentifier] = person1NumberOfVariantsTested quantityOfVariantsTestedMap[person2GenomeIdentifier] = person2NumberOfVariantsTested newOffspringGenomePairMonogenicDiseaseInfoObject := geneticAnalysis.OffspringGenomePairMonogenicDiseaseInfo{ ProbabilityOffspringHasDiseaseIsKnown: offspringHasDiseaseProbabilityIsKnown, ProbabilityOffspringHasDisease: percentageProbabilityOffspringHasDisease, ProbabilityOffspringHasVariantIsKnown: offspringHasAVariantProbabilityIsKnown, ProbabilityOffspringHasVariant: percentageProbabilityOffspringHasVariant, VariantsInfoMap: offspringVariantsInfoMap, } genomePairIdentifier := helpers.JoinTwo16ByteArrays(person1GenomeIdentifier, person2GenomeIdentifier) offspringMonogenicDiseaseInfoMap[genomePairIdentifier] = newOffspringGenomePairMonogenicDiseaseInfoObject return nil } err = addGenomePairInfoToDiseaseMaps(pair1Person1GenomeIdentifier, pair1Person2GenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addGenomePairInfoToDiseaseMaps(pair2Person1GenomeIdentifier, pair2Person2GenomeIdentifier) if (err != nil) { return false, "", err } } newOffspringMonogenicDiseaseInfoObject := geneticAnalysis.OffspringMonogenicDiseaseInfo{ QuantityOfVariantsTestedMap: quantityOfVariantsTestedMap, MonogenicDiseaseInfoMap: offspringMonogenicDiseaseInfoMap, } // 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 (len(offspringMonogenicDiseaseInfoMap) >= 2){ // Conflicts are only possible if two genome pairs with information about this disease exist checkIfConflictExists := func()(bool, error){ probabilityOffspringHasDiseaseIsKnown := false probabilityOffspringHasDisease := 0 probabilityOffspringHasVariantIsKnown := false probabilityOffspringHasVariant := 0 offspringVariantsInfoMap := make(map[[3]byte]geneticAnalysis.OffspringMonogenicDiseaseVariantInfo) firstItemReached := false for _, offspringMonogenicDiseaseInfoObject := range offspringMonogenicDiseaseInfoMap{ currentProbabilityOffspringHasDiseaseIsKnown := offspringMonogenicDiseaseInfoObject.ProbabilityOffspringHasDiseaseIsKnown currentProbabilityOffspringHasDisease := offspringMonogenicDiseaseInfoObject.ProbabilityOffspringHasDisease currentProbabilityOffspringHasVariantIsKnown := offspringMonogenicDiseaseInfoObject.ProbabilityOffspringHasVariantIsKnown currentProbabilityOffspringHasVariant := offspringMonogenicDiseaseInfoObject.ProbabilityOffspringHasVariant currentOffspringVariantsInfoMap := offspringMonogenicDiseaseInfoObject.VariantsInfoMap if (firstItemReached == false){ probabilityOffspringHasDiseaseIsKnown = currentProbabilityOffspringHasDiseaseIsKnown probabilityOffspringHasDisease = currentProbabilityOffspringHasDisease probabilityOffspringHasVariantIsKnown = currentProbabilityOffspringHasVariantIsKnown probabilityOffspringHasVariant = currentProbabilityOffspringHasVariant offspringVariantsInfoMap = currentOffspringVariantsInfoMap firstItemReached = true continue } if (currentProbabilityOffspringHasDiseaseIsKnown != probabilityOffspringHasDiseaseIsKnown){ return true, nil } if (currentProbabilityOffspringHasDisease != probabilityOffspringHasDisease){ return true, nil } if (currentProbabilityOffspringHasVariantIsKnown != probabilityOffspringHasVariantIsKnown){ return true, nil } if (currentProbabilityOffspringHasVariant != probabilityOffspringHasVariant){ return true, nil } areEqual := maps.Equal(currentOffspringVariantsInfoMap, offspringVariantsInfoMap) if (areEqual == false){ return true, nil } } return false, nil } conflictExists, err := checkIfConflictExists() if (err != nil) { return false, "", err } newOffspringMonogenicDiseaseInfoObject.ConflictExists = conflictExists } offspringMonogenicDiseasesMap[diseaseName] = newOffspringMonogenicDiseaseInfoObject } newCoupleAnalysis.MonogenicDiseasesMap = offspringMonogenicDiseasesMap // Step 2: Polygenic Diseases polygenicDiseaseObjectsList, err := polygenicDiseases.GetPolygenicDiseaseObjectsList() if (err != nil){ return false, "", err } // Map Structure: Disease Name -> OffspringPolygenicDiseaseInfo offspringPolygenicDiseasesMap := make(map[string]geneticAnalysis.OffspringPolygenicDiseaseInfo) for _, diseaseObject := range polygenicDiseaseObjectsList{ diseaseName := diseaseObject.DiseaseName // This map stores the polygenic disease info for each genome pair // Map Structure: Genome Pair Identifier -> OffspringGenomePairPolygenicDiseaseInfo offspringPolygenicDiseaseInfoMap := make(map[[32]byte]geneticAnalysis.OffspringGenomePairPolygenicDiseaseInfo) // This will calculate the disease info for the offspring from the two specified genomes to the diseases map // It then adds the pair entry to the offspringPolygenicDiseaseInfoMap addGenomePairDiseaseInfoToDiseaseMap := func(person1GenomeIdentifier [16]byte, person2GenomeIdentifier [16]byte)error{ person1LocusValuesMap, exists := person1GenomesMap[person1GenomeIdentifier] if (exists == false){ return errors.New("addGenomePairDiseaseInfoToDiseaseMap called with unknown person1GenomeIdentifier.") } person2LocusValuesMap, exists := person2GenomesMap[person2GenomeIdentifier] if (exists == false){ return errors.New("addGenomePairDiseaseInfoToDiseaseMap called with unknown person2GenomeIdentifier.") } neuralNetworkExists, anyOffspringLocusKnown, offspringAverageRiskScore, accuracyRangesMap, predictedRiskScoresList, quantityOfLociKnown, quantityOfParentalPhasedLoci, err := GetOffspringPolygenicDiseaseAnalysis(diseaseObject, person1LocusValuesMap, person2LocusValuesMap) if (err != nil) { return err } if (neuralNetworkExists == false){ // We cannot analyze this disease return nil } if (anyOffspringLocusKnown == false){ // We have no information about this genome pair's disease risk // We don't add this genome pair's disease info to the diseaseInfoMap return nil } newOffspringGenomePairPolygenicDiseaseInfo := geneticAnalysis.OffspringGenomePairPolygenicDiseaseInfo{ OffspringAverageRiskScore: offspringAverageRiskScore, PredictionConfidenceRangesMap: accuracyRangesMap, QuantityOfLociKnown: quantityOfLociKnown, QuantityOfParentalPhasedLoci: quantityOfParentalPhasedLoci, SampleOffspringRiskScoresList: predictedRiskScoresList, } genomePairIdentifier := helpers.JoinTwo16ByteArrays(person1GenomeIdentifier, person2GenomeIdentifier) offspringPolygenicDiseaseInfoMap[genomePairIdentifier] = newOffspringGenomePairPolygenicDiseaseInfo return nil } err = addGenomePairDiseaseInfoToDiseaseMap(pair1Person1GenomeIdentifier, pair1Person2GenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addGenomePairDiseaseInfoToDiseaseMap(pair2Person1GenomeIdentifier, pair2Person2GenomeIdentifier) if (err != nil) { return false, "", err } } if (len(offspringPolygenicDiseaseInfoMap) == 0){ // No disease analysis was performed continue } newOffspringPolygenicDiseaseInfoObject := geneticAnalysis.OffspringPolygenicDiseaseInfo{ PolygenicDiseaseInfoMap: offspringPolygenicDiseaseInfoMap, } if (len(offspringPolygenicDiseaseInfoMap) >= 2){ // We check for conflicts // Conflicts are only possible if two genome pairs with disease info for this disease exist checkIfConflictExists := func()(bool, error){ currentGenomePairPolygenicDiseaseInfo := geneticAnalysis.OffspringGenomePairPolygenicDiseaseInfo{} firstItemReached := false for _, genomePairDiseaseInfoObject := range offspringPolygenicDiseaseInfoMap{ if (firstItemReached == false){ currentGenomePairPolygenicDiseaseInfo = genomePairDiseaseInfoObject firstItemReached = true continue } areEqual := reflect.DeepEqual(genomePairDiseaseInfoObject, currentGenomePairPolygenicDiseaseInfo) if (areEqual == false){ return true, nil } } return false, nil } conflictExists, err := checkIfConflictExists() if (err != nil) { return false, "", err } newOffspringPolygenicDiseaseInfoObject.ConflictExists = conflictExists } offspringPolygenicDiseasesMap[diseaseName] = newOffspringPolygenicDiseaseInfoObject } newCoupleAnalysis.PolygenicDiseasesMap = offspringPolygenicDiseasesMap // Step 3: Traits traitObjectsList, err := traits.GetTraitObjectsList() if (err != nil) { return false, "", err } // Map Structure: Trait Name -> Trait Info Object offspringDiscreteTraitsMap := make(map[string]geneticAnalysis.OffspringDiscreteTraitInfo) // Map Structure: Trait Name -> Trait Info Object offspringNumericTraitsMap := make(map[string]geneticAnalysis.OffspringNumericTraitInfo) for _, traitObject := range traitObjectsList{ traitName := traitObject.TraitName traitIsDiscreteOrNumeric := traitObject.DiscreteOrNumeric if (traitIsDiscreteOrNumeric == "Discrete"){ // This map stores the trait info for each genome pair // Map Structure: Genome Pair Identifier -> OffspringGenomePairDiscreteTraitInfo offspringTraitInfoMap := make(map[[32]byte]geneticAnalysis.OffspringGenomePairDiscreteTraitInfo) // This will add the offspring trait information for the provided genome pair to the offspringTraitInfoMap addGenomePairTraitInfoToOffspringMap := func(person1GenomeIdentifier [16]byte, person2GenomeIdentifier [16]byte)error{ person1LocusValuesMap, exists := person1GenomesMap[person1GenomeIdentifier] if (exists == false){ return errors.New("addGenomePairTraitInfoToOffspringMap called with unknown person1GenomeIdentifier.") } person2LocusValuesMap, exists := person2GenomesMap[person2GenomeIdentifier] if (exists == false){ return errors.New("addGenomePairTraitInfoToOffspringMap called with unknown person2GenomeIdentifier.") } newOffspringGenomePairTraitInfo := geneticAnalysis.OffspringGenomePairDiscreteTraitInfo{} neuralNetworkExists, neuralNetworkAnalysisExists, outcomeProbabilitiesMap, averagePredictionConfidence, quantityOfLociTested, quantityOfParentalPhasedLoci, err := GetOffspringDiscreteTraitAnalysis_NeuralNetwork(traitObject, person1LocusValuesMap, person2LocusValuesMap) if (err != nil) { return err } if (neuralNetworkExists == true){ newOffspringGenomePairTraitInfo.NeuralNetworkExists = true if (neuralNetworkAnalysisExists == true){ newOffspringGenomePairTraitInfo.NeuralNetworkAnalysisExists = true newOffspringGenomePairTraitInfo_NeuralNetwork := geneticAnalysis.OffspringGenomePairDiscreteTraitInfo_NeuralNetwork{ OffspringOutcomeProbabilitiesMap: outcomeProbabilitiesMap, AverageConfidence: averagePredictionConfidence, QuantityOfLociKnown: quantityOfLociTested, QuantityOfParentalPhasedLoci: quantityOfParentalPhasedLoci, } newOffspringGenomePairTraitInfo.NeuralNetworkAnalysis = newOffspringGenomePairTraitInfo_NeuralNetwork } } anyRulesExist, rulesAnalysisExists, quantityOfRulesTested, quantityOfLociKnown, offspringProbabilityOfPassingRulesMap, offspringOutcomeProbabilitiesMap, err := GetOffspringDiscreteTraitAnalysis_Rules(traitObject, person1LocusValuesMap, person2LocusValuesMap) if (err != nil) { return err } if (anyRulesExist == true){ newOffspringGenomePairTraitInfo.RulesExist = true if (rulesAnalysisExists == true){ newOffspringGenomePairTraitInfo.RulesAnalysisExists = true newOffspringGenomePairTraitInfo_Rules := geneticAnalysis.OffspringGenomePairDiscreteTraitInfo_Rules{ QuantityOfRulesTested: quantityOfRulesTested, QuantityOfLociKnown: quantityOfLociKnown, ProbabilityOfPassingRulesMap: offspringProbabilityOfPassingRulesMap, OffspringOutcomeProbabilitiesMap: offspringOutcomeProbabilitiesMap, } newOffspringGenomePairTraitInfo.RulesAnalysis = newOffspringGenomePairTraitInfo_Rules } } genomePairIdentifier := helpers.JoinTwo16ByteArrays(person1GenomeIdentifier, person2GenomeIdentifier) offspringTraitInfoMap[genomePairIdentifier] = newOffspringGenomePairTraitInfo return nil } err = addGenomePairTraitInfoToOffspringMap(pair1Person1GenomeIdentifier, pair1Person2GenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addGenomePairTraitInfoToOffspringMap(pair2Person1GenomeIdentifier, pair2Person2GenomeIdentifier) if (err != nil) { return false, "", err } } newOffspringTraitInfoObject := geneticAnalysis.OffspringDiscreteTraitInfo{ TraitInfoMap: offspringTraitInfoMap, } if (len(offspringTraitInfoMap) >= 2){ // We check for conflicts // Conflicts are only possible if two genome pairs exist with information about the trait checkIfConflictExists := func()(bool, error){ // We check for conflicts between each genome pair's outcome scores and trait rules maps genomePairTraitInfoObject := geneticAnalysis.OffspringGenomePairDiscreteTraitInfo{} firstItemReached := false for _, currentGenomePairTraitInfoObject := range offspringTraitInfoMap{ if (firstItemReached == false){ genomePairTraitInfoObject = currentGenomePairTraitInfoObject firstItemReached = true continue } areEqual := reflect.DeepEqual(genomePairTraitInfoObject, currentGenomePairTraitInfoObject) if (areEqual == false){ return true, nil } } return false, nil } conflictExists, err := checkIfConflictExists() if (err != nil) { return false, "", err } newOffspringTraitInfoObject.ConflictExists = conflictExists } offspringDiscreteTraitsMap[traitName] = newOffspringTraitInfoObject } else { // traitIsDiscreteOrNumeric == "Numeric" // This map stores the trait info for each genome pair // Map Structure: Genome Pair Identifier -> OffspringGenomePairNumericTraitInfo offspringTraitInfoMap := make(map[[32]byte]geneticAnalysis.OffspringGenomePairNumericTraitInfo) // This will add the offspring trait information for the provided genome pair to the offspringTraitInfoMap addGenomePairTraitInfoToOffspringMap := func(person1GenomeIdentifier [16]byte, person2GenomeIdentifier [16]byte)error{ person1LocusValuesMap, exists := person1GenomesMap[person1GenomeIdentifier] if (exists == false){ return errors.New("addGenomePairTraitInfoToOffspringMap called with unknown person1GenomeIdentifier.") } person2LocusValuesMap, exists := person2GenomesMap[person2GenomeIdentifier] if (exists == false){ return errors.New("addGenomePairTraitInfoToOffspringMap called with unknown person2GenomeIdentifier.") } neuralNetworkExists, neuralNetworkAnalysisExists, averageOutcome, predictionConfidenceRangesMap, sampleOffspringOutcomesList, quantityOfLociTested, quantityOfParentalPhasedLoci, err := GetOffspringNumericTraitAnalysis(traitObject, person1LocusValuesMap, person2LocusValuesMap) if (err != nil) { return err } if (neuralNetworkExists == false){ // Predictions are not possible for this trait return nil } if (neuralNetworkAnalysisExists == false){ // No locations for this trait exists in which both user's genomes contain information return nil } newOffspringGenomePairTraitInfo := geneticAnalysis.OffspringGenomePairNumericTraitInfo{ OffspringAverageOutcome: averageOutcome, PredictionConfidenceRangesMap: predictionConfidenceRangesMap, QuantityOfParentalPhasedLoci: quantityOfParentalPhasedLoci, QuantityOfLociKnown: quantityOfLociTested, SampleOffspringOutcomesList: sampleOffspringOutcomesList, } genomePairIdentifier := helpers.JoinTwo16ByteArrays(person1GenomeIdentifier, person2GenomeIdentifier) offspringTraitInfoMap[genomePairIdentifier] = newOffspringGenomePairTraitInfo return nil } err = addGenomePairTraitInfoToOffspringMap(pair1Person1GenomeIdentifier, pair1Person2GenomeIdentifier) if (err != nil) { return false, "", err } if (genomePair2Exists == true){ err := addGenomePairTraitInfoToOffspringMap(pair2Person1GenomeIdentifier, pair2Person2GenomeIdentifier) if (err != nil) { return false, "", err } } newOffspringTraitInfoObject := geneticAnalysis.OffspringNumericTraitInfo{ TraitInfoMap: offspringTraitInfoMap, } if (len(offspringTraitInfoMap) >= 2){ // We check for conflicts // Conflicts are only possible if two genome pairs exist with information about the trait checkIfConflictExists := func()(bool, error){ // We check for conflicts between each genome pair's outcome scores and trait rules maps genomePairTraitInfoObject := geneticAnalysis.OffspringGenomePairNumericTraitInfo{} firstItemReached := false for _, currentGenomePairTraitInfoObject := range offspringTraitInfoMap{ if (firstItemReached == false){ genomePairTraitInfoObject = currentGenomePairTraitInfoObject firstItemReached = true continue } areEqual := reflect.DeepEqual(genomePairTraitInfoObject, currentGenomePairTraitInfoObject) if (areEqual == false){ return true, nil } } return false, nil } conflictExists, err := checkIfConflictExists() if (err != nil) { return false, "", err } newOffspringTraitInfoObject.ConflictExists = conflictExists } offspringNumericTraitsMap[traitName] = newOffspringTraitInfoObject } } newCoupleAnalysis.DiscreteTraitsMap = offspringDiscreteTraitsMap newCoupleAnalysis.NumericTraitsMap = offspringNumericTraitsMap analysisBytes, err := encoding.EncodeMessagePackBytes(newCoupleAnalysis) 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: // -bool: Probability offspring has disease is known // -int: Percentage probability offspring has disease (0-100) // -bool: Probability offspring has variant is known // -int: Percentage probability offspring has variant (0-100) // -error func GetOffspringMonogenicDiseaseProbabilities(dominantOrRecessive string, person1ProbabilityIsKnown bool, person1WillPassVariantPercentageProbability int, person2ProbabilityIsKnown bool, person2WillPassVariantPercentageProbability int)(bool, int, bool, int, error){ if (dominantOrRecessive != "Dominant" && dominantOrRecessive != "Recessive"){ return false, 0, false, 0, errors.New("GetOffspringMonogenicDiseaseProbabilities called with invalid dominantOrRecessive: " + dominantOrRecessive) } if (person1ProbabilityIsKnown == false && person2ProbabilityIsKnown == false){ return false, 0, false, 0, nil } if (person1ProbabilityIsKnown == true){ if (person1WillPassVariantPercentageProbability < 0 || person1WillPassVariantPercentageProbability > 100){ return false, 0, false, 0, errors.New("GetOffspringMonogenicDiseaseProbabilities called with invalid person1WillPassVariantProbability") } } if (person2ProbabilityIsKnown == true){ if (person2WillPassVariantPercentageProbability < 0 || person2WillPassVariantPercentageProbability > 100){ return false, 0, false, 0, errors.New("GetOffspringMonogenicDiseaseProbabilities called with invalid person2WillPassVariantProbability") } } if (person1ProbabilityIsKnown == false || person2ProbabilityIsKnown == false){ // Only 1 of the person's probabilities are known getPersonWhoHasVariantProbabilityOfPassingIt := func()int{ if (person1ProbabilityIsKnown == true){ return person1WillPassVariantPercentageProbability } return person2WillPassVariantPercentageProbability } personWhoHasVariantProbabilityOfPassingIt := getPersonWhoHasVariantProbabilityOfPassingIt() if (personWhoHasVariantProbabilityOfPassingIt == 100){ if (dominantOrRecessive == "Dominant"){ // We know the offspring will have the disease and will have a variant return true, 100, true, 100, nil } //dominantOrRecessive == "Recessive" // We don't know if the offspring will have the disease, but we know they will have a variant return false, 0, true, 100, nil } if (personWhoHasVariantProbabilityOfPassingIt == 0){ if (dominantOrRecessive == "Recessive"){ // We know the offspring will not have the disease, but we don't know if they will have a variant return true, 0, false, 0, nil } } // We don't know the probabilities that the offspring will have the disease or if they will have a variant return false, 0, false, 0, nil } person1WillPassVariantProbability := float64(person1WillPassVariantPercentageProbability)/100 person2WillPassVariantProbability := float64(person2WillPassVariantPercentageProbability)/100 // 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 1 passing a variant) + (Probability of person 2 passing a variant) - (Probability of offspring having disease) // A person with a variant may have the disease, or just be a carrier. probabilityOffspringHasVariant := person1WillPassVariantProbability + person2WillPassVariantProbability - (person1WillPassVariantProbability * person2WillPassVariantProbability) 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 true, percentageProbabilityOffspringHasVariant, true, 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 1 Passing a variant) * (Probability of person 2 passing a variant) probabilityOffspringHasDisease := person1WillPassVariantProbability * person2WillPassVariantProbability 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 true, percentageProbabilityOffspringHasDiseaseInt, true, percentageProbabilityOffspringHasVariantInt, nil } //Outputs: // -bool: A neural network exists for this trait // -bool: Any loci tested (if false, no offspring polygenic disease analysis is known) // -int: Offspring Average Risk Score (Value between 0-10) // -map[int]float64: Prediction accuracy ranges map // -Map Structure: Probability prediction is accurate (X) -> Distance from prediction that must be travelled in both directions to // create a range in which the true value will fall into, X% of the time // -[]int: Sample offspring risks scores list // -int: Quantity of loci known // -int: Quantity of parental phased loci // -error func GetOffspringPolygenicDiseaseAnalysis(diseaseObject polygenicDiseases.PolygenicDisease, person1LocusValuesMap map[int64]locusValue.LocusValue, person2LocusValuesMap map[int64]locusValue.LocusValue)(bool, bool, int, map[int]float64, []int, int, int, error){ diseaseName := diseaseObject.DiseaseName modelExists := trainedPredictionModels.CheckIfAttributeNeuralNetworkExists(diseaseName) if (modelExists == false){ // Prediction is not possible for this trait return false, false, 0, nil, nil, 0, 0, nil } if (len(person1LocusValuesMap) == 0){ return true, false, 0, nil, nil, 0, 0, nil } if (len(person2LocusValuesMap) == 0){ return true, false, 0, nil, nil, 0, 0, nil } diseaseLociList := diseaseObject.LociList // First we count up the quantity of parental phased loci // We only count the quantity of phased loci for loci which are known for both parents quantityOfParentalPhasedLoci := 0 for _, rsID := range diseaseLociList{ person1LocusValue, exists := person1LocusValuesMap[rsID] if (exists == false){ continue } person2LocusValue, exists := person2LocusValuesMap[rsID] if (exists == false){ continue } person1LocusIsPhased := person1LocusValue.LocusIsPhased if (person1LocusIsPhased == true){ quantityOfParentalPhasedLoci += 1 } person2LocusIsPhased := person2LocusValue.LocusIsPhased if (person2LocusIsPhased == true){ quantityOfParentalPhasedLoci += 1 } } // We create 100 prospective offspring genomes. anyLocusValueExists, prospectiveOffspringGenomesList, err := getProspectiveOffspringGenomesList(diseaseLociList, person1LocusValuesMap, person2LocusValuesMap) if (err != nil) { return false, false, 0, nil, nil, 0, 0, err } if (anyLocusValueExists == false){ return true, false, 0, nil, nil, 0, 0, nil } // A list of predicted risk scores for each offspring predictedRiskScoresList := make([]int, 0) accuracyRangesMap := make(map[int]float64) quantityOfLociTested := 0 for index, offspringGenomeMap := range prospectiveOffspringGenomesList{ neuralNetworkExists, predictionIsKnown, predictedRiskScore, predictionAccuracyRangesMap, currentQuantityOfLociTested, _, err := createPersonGeneticAnalysis.GetPersonGenomePolygenicDiseaseAnalysis(diseaseObject, offspringGenomeMap, false) if (err != nil){ return false, false, 0, nil, nil, 0, 0, err } if (neuralNetworkExists == false){ return false, false, 0, nil, nil, 0, 0, errors.New("GetGenomeNumericTraitAnalysis claiming that neural network doesn't exist when we already checked.") } if (predictionIsKnown == false){ return false, false, 0, nil, nil, 0, 0, errors.New("GetGenomeNumericTraitAnalysis claiming that prediction is impossible when we already know at least 1 locus value exists for trait.") } predictedRiskScoresList = append(predictedRiskScoresList, predictedRiskScore) if (index == 0){ // These values should be the same for each predicted offspring accuracyRangesMap = predictionAccuracyRangesMap quantityOfLociTested = currentQuantityOfLociTested } } // We calculate the average predicted risk score outcomesSum := 0 for _, predictedRiskScore := range predictedRiskScoresList{ outcomesSum += predictedRiskScore } averageRiskScore := outcomesSum/100 return true, true, averageRiskScore, accuracyRangesMap, predictedRiskScoresList, quantityOfLociTested, quantityOfParentalPhasedLoci, nil } //Outputs: // -bool: A neural network exists for this trait // -bool: Analysis exists (at least 1 locus exists for this analysis from both people's genomes // -map[string]int: Outcome probabilities map // Map Structure: Outcome Name -> Offspring probability of outcome // -int: Average prediction confidence (the average prediction confidence for all prospective offspring) // -int: Quantity of loci tested // -int: Quantity of parental phased loci // -error func GetOffspringDiscreteTraitAnalysis_NeuralNetwork(traitObject traits.Trait, person1LocusValuesMap map[int64]locusValue.LocusValue, person2LocusValuesMap map[int64]locusValue.LocusValue)(bool, bool, map[string]int, int, int, int, error){ traitName := traitObject.TraitName traitIsDiscreteOrNumeric := traitObject.DiscreteOrNumeric if (traitIsDiscreteOrNumeric != "Discrete"){ return false, false, nil, 0, 0, 0, errors.New("GetOffspringDiscreteTraitAnalysis_NeuralNetwork called with non-discrete trait.") } modelExists := trainedPredictionModels.CheckIfAttributeNeuralNetworkExists(traitName) if (modelExists == false){ // Neural network prediction is not possible for this trait return false, false, nil, 0, 0, 0, nil } traitLociList := traitObject.LociList // First we count up the quantity of parental phased loci // We only count the quantity of phased loci for loci which are known for both parents quantityOfParentalPhasedLoci := 0 for _, rsID := range traitLociList{ person1LocusValue, exists := person1LocusValuesMap[rsID] if (exists == false){ continue } person2LocusValue, exists := person2LocusValuesMap[rsID] if (exists == false){ continue } person1LocusIsPhased := person1LocusValue.LocusIsPhased if (person1LocusIsPhased == true){ quantityOfParentalPhasedLoci += 1 } person2LocusIsPhased := person2LocusValue.LocusIsPhased if (person2LocusIsPhased == true){ quantityOfParentalPhasedLoci += 1 } } // Next, we create 100 prospective offspring genomes. anyLocusValueExists, prospectiveOffspringGenomesList, err := getProspectiveOffspringGenomesList(traitLociList, person1LocusValuesMap, person2LocusValuesMap) if (err != nil) { return false, false, nil, 0, 0, 0, err } if (anyLocusValueExists == false){ return true, false, nil, 0, 0, 0, nil } // Map Structure: Outcome Name -> Probability of outcome coming true // Because we are summing from 100 offspring, the count of outcomes is the same as the probability of an offspring having the outcome outcomeCountsMap := make(map[string]int) // This is a sum of each prediction's confidence predictionConfidencesSum := 0 quantityOfLociTested := 0 for index, offspringGenomeMap := range prospectiveOffspringGenomesList{ neuralNetworkExists, predictionIsKnown, predictedOutcome, predictionConfidence, currentQuantityOfLociTested, _, err := createPersonGeneticAnalysis.GetGenomeDiscreteTraitAnalysis_NeuralNetwork(traitObject, offspringGenomeMap, false) if (err != nil){ return false, false, nil, 0, 0, 0, err } if (neuralNetworkExists == false){ return false, false, nil, 0, 0, 0, errors.New("GetGenomeTraitAnalysis_NeuralNetwork claiming that neural network doesn't exist when we already checked.") } if (predictionIsKnown == false){ return false, false, nil, 0, 0, 0, errors.New("GetGenomeTraitAnalysis_NeuralNetwork claiming that prediction is impossible when we already know at least 1 locus value exists for trait.") } outcomeCountsMap[predictedOutcome] += 1 predictionConfidencesSum += predictionConfidence if (index == 0){ // This value should be the same for each predicted offspring quantityOfLociTested = currentQuantityOfLociTested } } averagePredictionConfidence := predictionConfidencesSum/100 return true, true, outcomeCountsMap, averagePredictionConfidence, quantityOfLociTested, quantityOfParentalPhasedLoci, nil } //Outputs: // -bool: Any rules exist (if false, rule-based prediction is not possible for this trait) // -bool: Rule-based analysis exists (if false, no offspring trait information is known, or there is an outcome tie for one of the offspring) // -int: Quantity of rules tested // -int: Quantity of loci known // -map[[3]byte]int: Offspring probability of passing rules map // Map Structure: Rule identifier -> Offspring probability of passing rule (1-100) // If a rule entry doesn't exist, we don't know the passes-rule probability for any of the offspring // -map[string]int: Offspring outcome probabilities map // Map Structure: Outcome Name -> Offspring probability of outcome (0-100) // -error func GetOffspringDiscreteTraitAnalysis_Rules(traitObject traits.Trait, person1LocusValuesMap map[int64]locusValue.LocusValue, person2LocusValuesMap map[int64]locusValue.LocusValue)(bool, bool, int, int, map[[3]byte]int, map[string]int, error){ traitRulesList := traitObject.RulesList if (len(traitRulesList) == 0){ return false, false, 0, 0, nil, nil, nil } if (len(person1LocusValuesMap) == 0){ return true, false, 0, 0, nil, nil, nil } if (len(person2LocusValuesMap) == 0){ return true, false, 0, 0, nil, nil, nil } // First, we create 100 prospective offspring genomes. traitLociList_Rules := traitObject.LociList_Rules anyLocusValueExists, prospectiveOffspringGenomesList, err := getProspectiveOffspringGenomesList(traitLociList_Rules, person1LocusValuesMap, person2LocusValuesMap) if (err != nil) { return false, false, 0, 0, nil, nil, err } if (anyLocusValueExists == false){ return true, false, 0, 0, nil, nil, nil } // Map Structure: Rule Identifier -> Number of offspring who pass the rule (out of 100 prospective offspring) // Because there are 100 offspring, this also represents the percentage probability that an offspring will pass the rule offspringPassesRulesCountMap := make(map[[3]byte]int) // This map stores the quantity of offspring who have each outcome // The probability an offspring will have this outcome is the same as the // quantity of offspring who have this outcome in our set of 100 randomly generated offspring // Map structure: Outcome name -> quantity of offspring who have this outcome outcomeCountsMap := make(map[string]int) quantityOfLociKnown := 0 for index, offspringGenomeMap := range prospectiveOffspringGenomesList{ // Now we get outcome prediction for prospective offspring anyRulesExist, quantityOfRulesTested, currentQuantityOfLociKnown, offspringPassesRulesMap, predictionOutcomeIsKnown, predictedOutcome, err := createPersonGeneticAnalysis.GetGenomeDiscreteTraitAnalysis_Rules(traitObject, offspringGenomeMap, false) if (err != nil) { return false, false, 0, 0, nil, nil, err } if (anyRulesExist == false){ return false, false, 0, 0, nil, nil, errors.New("GetGenomeTraitAnalysis_Rules returning noRulesExists when we already checked and trait rules do in-fact exist.") } if (quantityOfRulesTested == 0){ // This will be the same for each of the 100 generated offspring // No analysis is possible. return true, false, 0, currentQuantityOfLociKnown, nil, nil, nil } if (index == 0){ // currentQuantityOfLociKnown will be the same for each prospective offspring quantityOfLociKnown = currentQuantityOfLociKnown } for ruleIdentifier, genomePassesRule := range offspringPassesRulesMap{ if (genomePassesRule == true){ offspringPassesRulesCountMap[ruleIdentifier] += 1 } } if (predictionOutcomeIsKnown == false){ // There was a tie between outcomes for this offspring // We can't predict anything about this trait for this couple using rules // This is why we need to create rules which make it unlikely for a tie between outcomes to occur. return true, false, 0, 0, nil, nil, nil } outcomeCountsMap[predictedOutcome] += 1 } quantityOfRulesTested := len(offspringPassesRulesCountMap) return true, true, quantityOfRulesTested, quantityOfLociKnown, offspringPassesRulesCountMap, outcomeCountsMap, nil } //Outputs: // -bool: A neural network exists for this trait // -bool: Analysis exists (at least 1 locus exists for this analysis from both people's genomes // -float64: Average outcome for offspring // -map[int]float64: Prediction accuracy ranges map // -Map Structure: Probability prediction is accurate (X) -> Distance from predictoin that must be travelled in both directions to // create a range in which the true value will fall into, X% of the time // -[]float64: A list of 100 offspring outcomes // -int: Quantity of loci known // -int: Quantity of parental phased loci // -error func GetOffspringNumericTraitAnalysis(traitObject traits.Trait, person1LocusValuesMap map[int64]locusValue.LocusValue, person2LocusValuesMap map[int64]locusValue.LocusValue)(bool, bool, float64, map[int]float64, []float64, int, int, error){ traitName := traitObject.TraitName traitIsDiscreteOrNumeric := traitObject.DiscreteOrNumeric if (traitIsDiscreteOrNumeric != "Numeric"){ return false, false, 0, nil, nil, 0, 0, errors.New("GetOffspringNumericTraitAnalysis called with non-numeric trait.") } modelExists := trainedPredictionModels.CheckIfAttributeNeuralNetworkExists(traitName) if (modelExists == false){ // Prediction is not possible for this trait return false, false, 0, nil, nil, 0, 0, nil } traitLociList := traitObject.LociList // First we count up the quantity of parental phased loci // We only count the quantity of phased loci for loci which are known for both parents quantityOfParentalPhasedLoci := 0 for _, rsID := range traitLociList{ person1LocusValue, exists := person1LocusValuesMap[rsID] if (exists == false){ continue } person2LocusValue, exists := person2LocusValuesMap[rsID] if (exists == false){ continue } person1LocusIsPhased := person1LocusValue.LocusIsPhased if (person1LocusIsPhased == true){ quantityOfParentalPhasedLoci += 1 } person2LocusIsPhased := person2LocusValue.LocusIsPhased if (person2LocusIsPhased == true){ quantityOfParentalPhasedLoci += 1 } } // Next, we create 100 prospective offspring genomes. anyLocusValueExists, prospectiveOffspringGenomesList, err := getProspectiveOffspringGenomesList(traitLociList, person1LocusValuesMap, person2LocusValuesMap) if (err != nil) { return false, false, 0, nil, nil, 0, 0, err } if (anyLocusValueExists == false){ return true, false, 0, nil, nil, 0, 0, nil } // A list of predicted outcomes for each offspring predictedOutcomesList := make([]float64, 0) accuracyRangesMap := make(map[int]float64) quantityOfLociTested := 0 for index, offspringGenomeMap := range prospectiveOffspringGenomesList{ neuralNetworkExists, predictionIsKnown, predictedOutcome, predictionAccuracyRangesMap, currentQuantityOfLociTested, _, err := createPersonGeneticAnalysis.GetGenomeNumericTraitAnalysis(traitObject, offspringGenomeMap, false) if (err != nil){ return false, false, 0, nil, nil, 0, 0, err } if (neuralNetworkExists == false){ return false, false, 0, nil, nil, 0, 0, errors.New("GetGenomeNumericTraitAnalysis claiming that neural network doesn't exist when we already checked.") } if (predictionIsKnown == false){ return false, false, 0, nil, nil, 0, 0, errors.New("GetGenomeNumericTraitAnalysis claiming that prediction is impossible when we already know at least 1 locus value exists for trait.") } predictedOutcomesList = append(predictedOutcomesList, predictedOutcome) if (index == 0){ // These values should be the same for each predicted offspring accuracyRangesMap = predictionAccuracyRangesMap quantityOfLociTested = currentQuantityOfLociTested } } // We calculate the average predicted outcome outcomesSum := float64(0) for _, predictedOutcome := range predictedOutcomesList{ outcomesSum += predictedOutcome } averageOutcome := outcomesSum/100 return true, true, averageOutcome, accuracyRangesMap, predictedOutcomesList, quantityOfLociTested, quantityOfParentalPhasedLoci, nil } // This function will return a list of 100 prospective offspring genomes // Each genome represents an equal-probability offspring genome from both people's genomes // This function takes into account the effects of genetic linkage // Any locations which do not exist in both people's genomes will not be included // // TODO: The user should be able to choose how many prospective offspring to create in the settings. // More offspring will take longer, but will yield a more accurate trait analysis. // //Outputs: // -bool: Any locus value exists between both users // -[]map[int64]locusValue.LocusValue // -error func getProspectiveOffspringGenomesList(lociList []int64, person1LociMap map[int64]locusValue.LocusValue, person2LociMap map[int64]locusValue.LocusValue)(bool, []map[int64]locusValue.LocusValue, error){ // -We use randomness to generate the offspring genomes // -We want the results to be the same for each pair of people each time, so we have to seed our randomness generator // -This is necessary so that two people's analysis results do not change every time // -Instead, the same 2 people will produce the exact same result every time pseudorandomNumberGenerator := mathRand.New(mathRand.NewPCG(1, 2)) //Outputs: // -[]int64: A list of random breakpoints for this chromosome that are statistically accurate // -error getRandomChromosomeBreakpoints := func(chromosome int)([]int64, error){ getChromosomeLength := func()(int64, error){ // Approximate number of base pairs in each chromosome taken from: https://www.ncbi.nlm.nih.gov/books/NBK557784/ switch chromosome{ case 1:{ return 249000000, nil } case 2:{ return 243000000, nil } case 3:{ return 200000000, nil } case 4:{ return 192000000, nil } case 5:{ return 181000000, nil } case 6:{ return 170000000, nil } case 7:{ return 158000000, nil } case 8:{ return 146000000, nil } case 9:{ return 140000000, nil } case 10:{ return 135000000, nil } case 11:{ return 135000000, nil } case 12:{ return 132000000, nil } case 13:{ return 114000000, nil } case 14:{ return 106000000, nil } case 15:{ return 100000000, nil } case 16:{ return 89000000, nil } case 17:{ return 79000000, nil } case 18:{ return 76000000, nil } case 19:{ return 64000000, nil } case 20:{ return 62000000, nil } case 21:{ return 47000000, nil } case 22:{ return 50000000, nil } } chromosomeString := helpers.ConvertIntToString(chromosome) return 0, errors.New("getRandomChromosomeBreakpoints called with invalid chromosome: " + chromosomeString) } chromosomeLength, err := getChromosomeLength() if (err != nil) { return nil, err } listOfRandomBreakpoints := make([]int64, 0) // TODO: Take into account different recombination rate for each chromosome // TODO: There are also breakpoint hotspots which we need to account for // TODO: I read somewhere that recombination break points are less likely to occur within genes, // meaning they are more likely to occur at the gene boundaries (codons) // We step by 1,000,000 each time // It would be more realistic if we did it in 1 integer increments, but it would be slower for position := int64(0); position < chromosomeLength; position += 1_000_000{ //From Wikipedia: // A centimorgan (abbreviated cM) is a unit for measuring genetic linkage. // It is defined as the distance between chromosome positions (loci) for which the expected // average number of intervening chromosomal crossovers in a single generation is 0.01. // One centimorgan corresponds to about 1 million base pairs in humans on average // // A chromosomal crossover == recombination breakpoint // // For every 1,000,000 base pairs, there is a 0.01 probability that there is a breakpoint randomFloat := pseudorandomNumberGenerator.Float64() if (randomFloat <= 0.01){ // This has a 0.01, or 1% probability of being true listOfRandomBreakpoints = append(listOfRandomBreakpoints, position) } } return listOfRandomBreakpoints, nil } // Map Structure: rsID -> Locus Value offspringGenomesList := make([]map[int64]locusValue.LocusValue, 0) for i:=0; i < 100; i++{ // This map stores the chromosome breakpoints for person1 // Map Structure: Chromosome -> List of breakpoints person1ChromosomeBreakpointsMap := make(map[int][]int64) // This map stores the chromosome breakpoints for person2 // Map Structure: Chromosome -> List of breakpoints person2ChromosomeBreakpointsMap := make(map[int][]int64) // This stores the locus values for this prospective offspring // Map Structure: rsID -> Locus Value prospectiveOffspringGenome := make(map[int64]locusValue.LocusValue) for _, rsID := range lociList{ //Outputs: // -bool: Allele is known // -string: Locus base // -error getPersonAllele := func(personLociMap map[int64]locusValue.LocusValue, personBreakpointsMap map[int][]int64)(bool, string, error){ personLocusValue, exists := personLociMap[rsID] if (exists == false){ return false, "", nil } personLocusBase1 := personLocusValue.Base1Value personLocusBase2 := personLocusValue.Base2Value personLocusIsPhased := personLocusValue.LocusIsPhased if (personLocusBase1 == personLocusBase2){ // Phase doesn't matter return true, personLocusBase1, nil } if (personLocusIsPhased == false){ // Breakpoints are unnecessary // We either choose base 1 or 2 randomInt := pseudorandomNumberGenerator.IntN(2) if (randomInt == 1){ return true, personLocusBase1, nil } return true, personLocusBase2, nil } // We have a phased locus // We figure out which allele to use by seeing which allele gets inherited from our random breakpoints list // We figure out the chromosome and position of this locus locusMetadataExists, locusMetadataObject, err := locusMetadata.GetLocusMetadata(rsID) if (err != nil) { return false, "", err } if (locusMetadataExists == false){ rsIDString := helpers.ConvertInt64ToString(rsID) return false, "", errors.New("getProspectiveOffspringGenomesList called with unknown rsID: " + rsIDString) } locusPosition := locusMetadataObject.Position locusChromosome := locusMetadataObject.Chromosome getPersonChromosomeBreakpointsList := func()([]int64, error){ breakpointsList, exists := personBreakpointsMap[locusChromosome] if (exists == true){ return breakpointsList, nil } // We have to create a new breakpoints list newBreakpointsList, err := getRandomChromosomeBreakpoints(locusChromosome) if (err != nil) { return nil, err } personBreakpointsMap[locusChromosome] = newBreakpointsList return newBreakpointsList, nil } personBreakpointsList, err := getPersonChromosomeBreakpointsList() if (err != nil) { return false, "", err } getLocusListIndex := func()int{ for index, breakpointPosition := range personBreakpointsList{ if (int64(locusPosition) <= breakpointPosition){ return index } } index := len(personBreakpointsList) // This is reached if the final breakpoint in the list is less than the locus's position, or if there were no breakpoints return index } locusListIndex := getLocusListIndex() if (locusListIndex%2 == 0){ return true, personLocusBase1, nil } return true, personLocusBase2, nil } person1AlleleIsKnown, person1Allele, err := getPersonAllele(person1LociMap, person1ChromosomeBreakpointsMap) if (err != nil) { return false, nil, err } if (person1AlleleIsKnown == false){ continue } person2AlleleIsKnown, person2Allele, err := getPersonAllele(person2LociMap, person2ChromosomeBreakpointsMap) if (err != nil) { return false, nil, err } if (person2AlleleIsKnown == false){ continue } offspringLocusValue := locusValue.LocusValue{ Base1Value: person1Allele, Base2Value: person2Allele, LocusIsPhased: true, } prospectiveOffspringGenome[rsID] = offspringLocusValue } if (len(prospectiveOffspringGenome) == 0){ // We don't have any locations at which both people's genomes contain a locus value. return false, nil, nil } offspringGenomesList = append(offspringGenomesList, prospectiveOffspringGenome) } return true, offspringGenomesList, nil }