# How to detect the frequency of gene pairings in lists of genetic profiles

We have a list of >1000 gene profiles each consisting typically of 4 genes. Many profiles share pairs of identical genes. Can we derive clusters of profiles which share pairs or even trios of genes? So.. replacing the gene symbols (e.g. ABL1) by letters for simplicity, let the list of gene profiles be {{a,b,c,d},{a,b,e,f},{e,f,g,h}}, we wish to derive the clusters {{{a,b,c,d},{a,b,e,f}},{{a,b,e,f},{e,f,g,h}}.

• Does the order matter? I.e. would {{a,b,c,d},{e,f,g,h},{a,b,e,f}} also yield {{{a,b,c,d},{a,b,e,f}},{{a,b,e,f},{e,f,g,h}} or not? – march Jun 23 '17 at 23:02
• No, the order is unimportant. – Boris Jun 23 '17 at 23:05

You can do it in several ways using Subsets combined with Pick, Select, Cases etc:

ClearAll[pairingsF1, pairingsF2, pairingsF3, pairingsF4, pairingsF5]

pairingsF1 = Pick[#, Length[Intersection @@ #] >= 2 & /@ #]& @ Subsets[#, {2, ∞}] &;
pairingsF2 = Select[Subsets[#, {2, ∞}], Length[Intersection @@ #] >= 2 &] &;
pairingsF3 = Cases[Subsets[#, {2, ∞}], _?(Length[Intersection @@ #] >= 2 &)] &;
pairingsF4 = DeleteCases[Subsets[#, {2, ∞}], _?(Length[Intersection @@ #] < 2 &)] &;
pairingsF5 = If[Length[Intersection @@ #] < 2, ## &[], #] & /@ Subsets[#, {2, ∞}] &;


Examples:

pairingsF1 @ {{a, b, c, d}, {a, b, e, f}, {e, f, g, h}}


{{{a, b, c, d}, {a, b, e, f}}, {{a, b, e, f}, {e, f, g, h}}}

Equal @@ (#[{{a, b, c, d}, {a, b, e, f}, {e, f, g, h}}] & /@
{pairingsF1, pairingsF2, pairingsF3, pairingsF4, pairingsF5})


True

pairingsF1@{{a, b, c, d}, {e, b, c, f}, {g, c, h, d}}


{{{a, b, c, d}, {e, b, c, f}}, {{a, b, c, d}, {g, c, h, d}}}

pairingsF1@{{a, b, c, d}, {a, b, e, f}, {e, f, g, h}, {a, b, c, h}}


{{{a, b, c, d}, {a, b, e, f}},
{{a, b, c, d}, {a, b, c, h}},
{{a, b, e, f}, {e, f, g, h}},
{{a, b, e, f}, {a, b, c, h}},
{{a, b, c, d}, {a, b, e, f}, {a, b, c, h}}}

• it seems that testing with {{a, b, c, d}, {a, b, e, f}, {e, f, g, h}, {a, b, c, h}} gives one redundant pair {{e, f, g, h}, {a, b, c, h}}. – garej Jun 24 '17 at 6:15
• Doesn't seem to work on something a bit more complicated e.g. {{a, b, c, d}, {e, b, c, f}, {g, c, h, d}} yields {{{a, b, c, d}, {e, b, c, f}}, {{a, b, c, d}, {g, c, h, d}}, {{e, b, c, f}, {g, c, h, d}}, {{a, b, c, d}, {e, b, c, f}, {g, c, h, d}}}, What am I missing? – Boris Jun 24 '17 at 6:37
• Boris, @garej, the original answer was checking if the Intersection is non-empty to group a set of profiles. The updated version fixes the issue using the correct condition (i.e. Length[Intersection@#]>=2) for a set of profiles to form a cluster. – kglr Jun 24 '17 at 9:45
• @kglr, do you think Subsets would work with >1000 sequences? For OP it seems better to avoid Infinity and work with subsets of {2} elements. – garej Jun 24 '17 at 10:44
• @kglr,@garej - before I go the full monty on my 1000+ set of profiles I tested your routines on uneven subsets and it was OK - can you confirm this is OK. Your efforts appreciated, I'll acknowledge you in the MS! – Boris Jun 24 '17 at 11:13

Here's a way to do it with Association.

OP's example. You want to sort the profiles so that the keys (= pairs of genes) are in canonical order. If the profiles in gprof are known to be ordered, Sort may be omitted, but it's not a big expense.

gprof = {{a, b, c, d}, {a, b, e, f}, {e, f, g, h}};
clusters = Merge[
Function[p, AssociationMap[p &, Subsets[p, {2}]]] /@ Sort /@ gprof,
Join];
Select[clusters, Length[#] >= 2 &]
(*
<|{a, b} -> {{a, b, c, d}, {a, b, e, f}},
{e, f} -> {{a, b, e, f}, {e, f, g, h}}|>
*)


An example similar to the size indicated in the OP. It seems quick enough.

SeedRandom[0];
genes = RandomSample[GenomeData[], 2000];
gprof = Table[RandomSample[genes, 4], {1001}];
(clusters = Merge[
Function[p, AssociationMap[p &, Subsets[p, {2}]]] /@ Sort /@ gprof,
Join];
Select[clusters, Length[#] >= 2 &]) // AbsoluteTiming
(*
{0.038202,
<|{"LOC100147773", "NPHP4"} ->
{{"LOC100147773", "NPHP4", "PROCR", "RAB11FIP1"},
{"LOC100147773", "LOC389049", "NPHP4", "OR2V2"}},
{"LOC729449", "TRBV102"} ->
{{"LOC401431", "LOC729449", "POLQ", "TRBV102"},
{"FAM129C", "LOC729449", "SLC22A9", "TRBV102"}},
{"LOC440606", "TRAF5"} ->
{{"CDH19", "LOC100132381", "LOC440606", "TRAF5"},
{"LOC440606", "LOC728351", "LOC730179", "TRAF5"}},
{"LOC100134062", "ZNF606"} ->
{{"LOC100134062", "LOC100134243", "SSR4", "ZNF606"},
{"ACOT4", "LOC100133564", "LOC100134062", "ZNF606"}},
{"LOC100133652", "RNU6ATAC4P"} ->
{{"LOC100133652", "MEPE", "RNU6ATAC4P", "SNORA11"},
{"LOC100127887", "LOC100130622", "LOC100133652", "RNU6ATAC4P"}}|>}
*)


A bigger example with clusters of three or more and profiles of 4-6 genes.

SeedRandom[2];
genes = RandomSample[GenomeData[], 2000];
gprof = Table[RandomSample[genes, RandomInteger[{4, 6}]], {2000}];
(clusters = Merge[
Function[p, AssociationMap[p &, Subsets[p, {2}]]] /@ Sort /@ gprof,
Join];
Select[clusters, Length[#] >= 3 &]) // AbsoluteTiming
(*
{0.125036,
<|{"IFT81", "LOC100134117"} ->
{{"IFT81", "LOC100134117", "Tcag71263", "TMEM16J"},
{"IFT81", "LOC100134117", "LOC401400", "LOC732377", "LOR", "RCAN2"},
{"IFT81", "LOC100134117", "LOC653758", "MRX35", "NSX", "UPLP"}}|>}
*)


Well, in fact, just one cluster of size exactly three.

• @Micheal E2: Here are the first 10 in my list of profiles. Substituting this for gprof in clusters misses the SRSF2,BCORL1 pairings. BTW @ kglr, your routines work but run out of memory for somewhere less than 50 profiles. – Boris Jun 26 '17 at 18:54
• Here is the list {{"NRAS, KIT, KIT"}, {"TP53, TP53"}, {"KRAS, TP53"}, {"ASXL1, CBLC, \ DNMT3A, HRAS, IDH2, PHF6, PHF6, PHF6, RUNX1, SETBP1"}, {"HRAS, \ TET2"}, {"BCORL1, DNMT3A, NPM1, SRSF2, WT1"}, {"CBLB, SRSF2, BCORL1, \ BCORL1"}, {"TET2, TET2, KMT2A, BCORL1"}, {"CUX1, NOTCH1, TP53, \ BCORL1"}, {"TET2, TET2, SMC3, SRSF2, RUNX1, RUNX1, BCORL1"}}. As I paste this list I notice that the profiles turn up as strings. Maybe this is an artefact. I'll check and split the strings if necessary. – Boris Jun 26 '17 at 18:58
• Thats the problem! – Boris Jun 26 '17 at 19:05
• @Boris I would figure that if they show up as single strings (such as "NRAS, KIT, KIT"`) then they were not interpreted/imported correctly by Mathematica. -- By "Thats the problem!" do you mean it now works correctly after the strings are split? – Michael E2 Jun 26 '17 at 22:44
• @ Michael E2 Yes, after the split it works perfectly on list length 840 - I'll go higher later. A few sets exceed 50 profiles. Food for thought. i imported "TSV" and got the single strings - is there another way to import thus avoiding the 'split' step? Maybe its the formulation of the source data. Your input appreciated. – Boris Jun 27 '17 at 7:42