Handling of applied variants with and without modifications

mzLib/Omics/BioPolymer/IHasSequenceVariants.cs

@@ -47,13 +47,19 @@ public interface IHasSequenceVariants
     /// </summary>
     List<TruncationProduct> TruncationProducts { get; }
 
-    
+
     /// <summary>
-    /// Used to construct a new variant of the same type as the original and is called in <see cref="VariantApplication"/>
+    /// Constructs a new variant biopolymer instance from the original, applying the specified sequence variants and modifications.
+    /// The method sets both the full set of database sequence variants and the subset of applied variants, ensuring that only unapplied variants remain in the SequenceVariations list,
+    /// while applied variants are tracked in AppliedSequenceVariations. This enables accurate representation of both the original and variant-specific annotations.
     /// </summary>
-    /// <remarks>The generic structure enables proteins to produce proteins and RNA to produce RNA</remarks>
-    /// <returns></returns>
-    TBioPolymerType CreateVariant<TBioPolymerType>(string variantBaseSequence, TBioPolymerType original, IEnumerable<SequenceVariation> appliedSequenceVariants,
-        IEnumerable<TruncationProduct> applicableProteolysisProducts, IDictionary<int, List<Modification>> oneBasedModifications, string sampleNameForVariants)
-        where TBioPolymerType : IHasSequenceVariants;
+    TBioPolymerType CreateVariant<TBioPolymerType>(
+        string variantBaseSequence,
+        TBioPolymerType original,
+        IEnumerable<SequenceVariation> appliedSequenceVariants,
+        IEnumerable<TruncationProduct> applicableProteolysisProducts,
+        IDictionary<int,
+            List<Modification>> oneBasedModifications,
+        string sampleNameForVariants)
+    where TBioPolymerType : IHasSequenceVariants;
 }

mzLib/Omics/BioPolymer/VariantApplication.cs

@@ -1,4 +1,5 @@
-using MzLibUtil;
+using System.Linq;
+using MzLibUtil;
 using Omics.BioPolymer;
 using Omics.Modifications;
 
@@ -94,24 +95,36 @@ public static int RestoreModificationIndex(IHasSequenceVariants protein, int var
         /// <summary>
         /// Applies multiple variant changes to a protein sequence
         /// </summary>
-        public static List<TBioPolymerType> ApplyVariants<TBioPolymerType>(TBioPolymerType protein, IEnumerable<SequenceVariation> sequenceVariations, int maxAllowedVariantsForCombinitorics, int minAlleleDepth)
+        public static List<TBioPolymerType> ApplyVariants<TBioPolymerType>(TBioPolymerType protein, IEnumerable<SequenceVariation> sequenceVariations, int maxAllowedVariantsForCombinatorics, int minAlleleDepth)
             where TBioPolymerType : IHasSequenceVariants
         {
             List<SequenceVariation> uniqueEffectsToApply = sequenceVariations
                 .GroupBy(v => v.SimpleString())
                 .Select(x => x.First())
                 .Where(v => v.Description.Genotypes.Count > 0) // this is a VCF line
-                .OrderByDescending(v => v.OneBasedBeginPosition) // apply variants at the end of the protein sequence first
+                .OrderBy(v => v.OneBasedBeginPosition)
+                .ThenBy(v => v.OneBasedEndPosition)
+                .ThenBy(v => v.OriginalSequence ?? string.Empty)
+                .ThenBy(v => v.VariantSequence ?? string.Empty) // apply variants at the end of the protein sequence first
                 .ToList();
 
-            TBioPolymerType proteinCopy = protein.CreateVariant(protein.BaseSequence, protein, null, protein.TruncationProducts, protein.OneBasedPossibleLocalizedModifications, null);
-
-            // If there aren't any variants to apply, just return the base protein
+            // If there aren't any variants to apply, just return the original as-is
             if (uniqueEffectsToApply.Count == 0)
             {
-                return new List<TBioPolymerType> { proteinCopy };
+                return new List<TBioPolymerType> { protein };
             }
 
+            // Start from a normalized copy whose “unapplied” list is the DB list.
+            TBioPolymerType proteinCopy = protein.CreateVariant(
+                protein.BaseSequence,
+                protein,
+                null,                                  // none applied yet
+                protein.TruncationProducts,
+                protein.OneBasedPossibleLocalizedModifications,
+                null);
+
+
+
             HashSet<string> individuals = new HashSet<string>(uniqueEffectsToApply.SelectMany(v => v.Description.Genotypes.Keys));
             List<TBioPolymerType> variantProteins = new();
             List<TBioPolymerType> newVariantProteins = new();
@@ -121,7 +134,7 @@ public static List<TBioPolymerType> ApplyVariants<TBioPolymerType>(TBioPolymerTy
                 newVariantProteins.Clear();
                 newVariantProteins.Add(proteinCopy);
 
-                bool tooManyHeterozygousVariants = uniqueEffectsToApply.Count(v => v.Description.Heterozygous[individual]) > maxAllowedVariantsForCombinitorics;
+                bool tooManyHeterozygousVariants = uniqueEffectsToApply.Count(v => v.Description.Heterozygous[individual]) > maxAllowedVariantsForCombinatorics;
                 foreach (var variant in uniqueEffectsToApply)
                 {
                     bool variantAlleleIsInTheGenotype = variant.Description.Genotypes[individual].Contains(variant.Description.AlleleIndex.ToString()); // should catch the case where it's -1 if the INFO isn't from SnpEff
@@ -145,12 +158,12 @@ public static List<TBioPolymerType> ApplyVariants<TBioPolymerType>(TBioPolymerTy
                     {
                         if (isDeepAlternateAllele && isDeepReferenceAllele)
                         {
-                            if (newVariantProteins.Count == 1 && maxAllowedVariantsForCombinitorics > 0)
+                            if (newVariantProteins.Count == 1 && maxAllowedVariantsForCombinatorics > 0)
                             {
                                 TBioPolymerType variantProtein = ApplySingleVariant(variant, newVariantProteins[0], individual);
                                 newVariantProteins.Add(variantProtein);
                             }
-                            else if (maxAllowedVariantsForCombinitorics > 0)
+                            else if (maxAllowedVariantsForCombinatorics > 0)
                             {
                                 newVariantProteins[1] = ApplySingleVariant(variant, newVariantProteins[1], individual);
                             }
@@ -159,7 +172,7 @@ public static List<TBioPolymerType> ApplyVariants<TBioPolymerType>(TBioPolymerTy
                                 // no heterozygous variants
                             }
                         }
-                        else if (isDeepAlternateAllele && maxAllowedVariantsForCombinitorics > 0)
+                        else if (isDeepAlternateAllele && maxAllowedVariantsForCombinatorics > 0)
                         {
                             newVariantProteins = newVariantProteins.Select(p => ApplySingleVariant(variant, p, individual)).ToList();
                         }
@@ -176,7 +189,7 @@ public static List<TBioPolymerType> ApplyVariants<TBioPolymerType>(TBioPolymerTy
 
                         foreach (var ppp in newVariantProteins)
                         {
-                            if (isDeepAlternateAllele && maxAllowedVariantsForCombinitorics > 0 && isDeepReferenceAllele)
+                            if (isDeepAlternateAllele && maxAllowedVariantsForCombinatorics > 0 && isDeepReferenceAllele)
                             {
                                 // keep reference allele
                                 if (variant.Description.Genotypes[individual].Contains("0"))
@@ -187,7 +200,7 @@ public static List<TBioPolymerType> ApplyVariants<TBioPolymerType>(TBioPolymerTy
                                 // alternate allele (replace all, since in heterozygous with two alternates, both alternates are included)
                                 combinitoricProteins.Add(ApplySingleVariant(variant, ppp, individual));
                             }
-                            else if (isDeepAlternateAllele && maxAllowedVariantsForCombinitorics > 0)
+                            else if (isDeepAlternateAllele && maxAllowedVariantsForCombinatorics > 0)
                             {
                                 combinitoricProteins.Add(ApplySingleVariant(variant, ppp, individual));
                             }
@@ -253,7 +266,14 @@ private static TBioPolymerType ApplySingleVariant<TBioPolymerType>(SequenceVaria
             Dictionary<int, List<Modification>> adjustedModifications = AdjustModificationIndices(variantGettingApplied, variantSequence, protein);
             List<SequenceVariation> adjustedAppliedVariations = AdjustSequenceVariationIndices(variantGettingApplied, variantSequence, appliedVariations);
 
-            return protein.CreateVariant(variantSequence, protein, adjustedAppliedVariations, adjustedProteolysisProducts, adjustedModifications, individual);
+            // Pass BOTH lists: current DB list (possibly already filtered on the interim proteoform) + newly applied
+            return protein.CreateVariant(
+                variantSequence,
+                protein,
+                adjustedAppliedVariations,
+                adjustedProteolysisProducts,
+                adjustedModifications,
+                individual);
         }
 
         /// <summary>
@@ -412,20 +432,34 @@ private static Dictionary<int, List<Modification>> AdjustModificationIndices(Seq
             return mods;
         }
 
+        /// <summary>
+        /// Deterministic ordering for variant-based naming/IDs to avoid order-of-application artifacts.
+        /// </summary>
+        private static IEnumerable<SequenceVariation> OrderForNaming(IEnumerable<SequenceVariation> variations) =>
+            variations
+                .OrderBy(v => v.OneBasedBeginPosition)
+                .ThenBy(v => v.OneBasedEndPosition)
+                .ThenBy(v => v.OriginalSequence ?? string.Empty)
+                .ThenBy(v => v.VariantSequence ?? string.Empty);
+
         /// <summary>
         /// Format string to append to accession
         /// </summary>
         private static string CombineSimpleStrings(IEnumerable<SequenceVariation>? variations)
         {
-            return variations.IsNullOrEmpty() ? "" : string.Join("_", variations.Select(v => v.SimpleString()));
+            return variations.IsNullOrEmpty()
+                ? ""
+                : string.Join("_", OrderForNaming(variations!).Select(v => v.SimpleString()));
         }
 
         /// <summary>
         /// Format string to append to protein names
         /// </summary>
         public static string CombineDescriptions(IEnumerable<SequenceVariation>? variations)
         {
-            return variations.IsNullOrEmpty() ? "" : string.Join(", variant:", variations.Select(d => d.Description));
+            return variations.IsNullOrEmpty()
+                ? ""
+                : string.Join(", variant:", OrderForNaming(variations!).Select(v => v.Description));
         }
         /// <summary>
         /// Applies all possible combinations of the provided SequenceVariation list to the base TBioPolymerType object,
@@ -451,7 +485,13 @@ public static IEnumerable<TBioPolymerType> ApplyAllVariantCombinations<TBioPolym
             count++;
             if (count >= maxCombinations)
                 yield break;
-
+            // Sort variations to ensure deterministic combination order
+            variations = variations
+                .OrderBy(v => v.OneBasedBeginPosition)
+                .ThenBy(v => v.OneBasedEndPosition)
+                .ThenBy(v => v.OriginalSequence ?? string.Empty)
+                .ThenBy(v => v.VariantSequence ?? string.Empty)
+                .ToList();
             int n = variations.Count;
             for (int size = 1; size <= n; size++)
             {

mzLib/Proteomics/Protein/Protein.cs

@@ -133,30 +133,35 @@ public Protein(Protein originalProtein, string newBaseSequence)
         /// <param name="applicableProteolysisProducts"></param>
         /// <param name="oneBasedModifications"></param>
         /// <param name="sampleNameForVariants"></param>
-        public Protein(string variantBaseSequence, Protein protein, IEnumerable<SequenceVariation> appliedSequenceVariations,
-            IEnumerable<TruncationProduct> applicableProteolysisProducts, IDictionary<int, List<Modification>> oneBasedModifications, string sampleNameForVariants)
+        public Protein(
+            string variantBaseSequence, 
+            Protein protein, 
+            IEnumerable<SequenceVariation> appliedSequenceVariations,
+            IEnumerable<TruncationProduct> applicableProteolysisProducts, 
+            IDictionary<int, List<Modification>> oneBasedModifications, 
+            string sampleNameForVariants)
             : this(
-                  variantBaseSequence,
-                  VariantApplication.GetAccession(protein, appliedSequenceVariations),
-                  organism: protein.Organism,
-                  geneNames: new List<Tuple<string, string>>(protein.GeneNames),
-                  oneBasedModifications: oneBasedModifications != null ? oneBasedModifications.ToDictionary(x => x.Key, x => x.Value) : new Dictionary<int, List<Modification>>(),
-                  proteolysisProducts: new List<TruncationProduct>(applicableProteolysisProducts ?? new List<TruncationProduct>()),
-                  name: VariantApplication.GetVariantName(protein.Name, appliedSequenceVariations),
-                  fullName: VariantApplication.GetVariantName(protein.FullName, appliedSequenceVariations), 
-                  isDecoy: protein.IsDecoy,
-                  isContaminant: protein.IsContaminant,
-                  databaseReferences: new List<DatabaseReference>(protein.DatabaseReferences),
-                  sequenceVariations: new List<SequenceVariation>(protein.SequenceVariations),
-                  disulfideBonds: new List<DisulfideBond>(protein.DisulfideBonds),
-                  spliceSites: new List<SpliceSite>(protein.SpliceSites),
-                  databaseFilePath: protein.DatabaseFilePath,
-                  dataset: protein.DatasetEntryTag,
-                  created: protein.CreatedEntryTag,
-                  modified: protein.ModifiedEntryTag,
-                  version: protein.VersionEntryTag,
-                  xmlns: protein.XmlnsEntryTag,
-                  uniProtSequenceAttributes: protein.UniProtSequenceAttributes)
+                variantBaseSequence,
+                VariantApplication.GetAccession(protein, appliedSequenceVariations),
+                organism: protein.Organism,
+                geneNames: new List<Tuple<string, string>>(protein.GeneNames),
+                oneBasedModifications: oneBasedModifications != null ? oneBasedModifications.ToDictionary(x => x.Key, x => x.Value) : new Dictionary<int, List<Modification>>(),
+                proteolysisProducts: new List<TruncationProduct>(applicableProteolysisProducts ?? new List<TruncationProduct>()),
+                name: VariantApplication.GetVariantName(protein.Name, appliedSequenceVariations),
+                fullName: VariantApplication.GetVariantName(protein.FullName, appliedSequenceVariations),
+                isDecoy: protein.IsDecoy,
+                isContaminant: protein.IsContaminant,
+                databaseReferences: new List<DatabaseReference>(protein.DatabaseReferences),
+                sequenceVariations: new List<SequenceVariation>(protein.SequenceVariations),
+                disulfideBonds: new List<DisulfideBond>(protein.DisulfideBonds),
+                spliceSites: new List<SpliceSite>(protein.SpliceSites),
+                databaseFilePath: protein.DatabaseFilePath,
+                dataset: protein.DatasetEntryTag,
+                created: protein.CreatedEntryTag,
+                modified: protein.ModifiedEntryTag,
+                version: protein.VersionEntryTag,
+                xmlns: protein.XmlnsEntryTag,
+                uniProtSequenceAttributes: protein.UniProtSequenceAttributes)
         {
             NonVariantProtein = protein.ConsensusVariant as Protein;
             OriginalNonVariantModifications = ConsensusVariant.OriginalNonVariantModifications;
@@ -603,8 +608,9 @@ public TBioPolymerType CreateVariant<TBioPolymerType>(string variantBaseSequence
             if (original is not Protein originalProtein)
                 throw new ArgumentException("The original BioPolymer must be Protein to create a protein variant");
 
-            var variantProtein =  new Protein(variantBaseSequence, originalProtein, appliedSequenceVariants, 
+            var variantProtein = new Protein(variantBaseSequence, originalProtein, appliedSequenceVariants,
                 applicableProteolysisProducts, oneBasedModifications, sampleNameForVariants);
+
             return (TBioPolymerType)(IHasSequenceVariants)variantProtein;
         }
         #endregion