This script run a PubMed-Comorbidities pipeline using the following characteristics:
In [1]:
addprocs(2);
In [2]:
using Revise #used during development to detect changes in module
using PubMedMiner
#Settings
const mh = "Epilepsy"
const concepts = ("Disease or Syndrome", "Mental or Behavioral Dysfunction", "Neoplastic Process");
In [3]:
overwrite = false
@time save_semantic_occurrences(mh, concepts...; overwrite = overwrite)
In [5]:
using FreqTables
occurrence_df = get_semantic_occurrences_df(mh, concepts...)
@time mesh_frequencies = freqtable(occurrence_df, :pmid, :descriptor);
info("Found ", size(occurrence_df, 1), " related descriptors")
In [6]:
using PlotlyJS
using NamedArrays
# Visualize frequency
topn = 50
mesh_counts = vec(sum(mesh_frequencies, 1))
count_perm = sortperm(mesh_counts, rev=true)
mesh_names = collect(keys(mesh_frequencies.dicts[2]))
#traces
#most frequent is epilepsy - remove from plot for better scaling
freq_trace = PlotlyJS.bar(; x = mesh_names[count_perm[2:topn]], y= mesh_counts[count_perm[2:topn]], marker_color="orange")
data = [freq_trace]
layout = Layout(;title="$(topn)-Most Frequent MeSH ",
showlegend=false,
margin= Dict(:t=> 70, :r=> 0, :l=> 50, :b=>200),
xaxis_tickangle = 90,)
plot(data, layout)
Out[6]:
In [10]:
using BCBIStats.COOccur
using StatsBase
#co-occurrance matrix - only for topp MeSH
# min_frequency = 5 -- alternatively compute topn based on min-frequency
top_occ = mesh_frequencies.array[:, count_perm[2:topn]]
top_mesh_labels = mesh_names[count_perm[2:topn]]
top_occ_sp = sparse(top_occ)
top_coo_sp = top_occ_sp' * top_occ_sp
#Point Mutual Information
pmi_sp = BCBIStats.COOccur.pmi_mat(top_coo_sp)
#chi2
top_chi2= BCBIStats.COOccur.chi2_mat(top_occ, min_freq=0);
In [12]:
using PlotlyJSFactory
p = create_chord_plot(top_coo_sp, labels = top_mesh_labels)
relayout!(p, title="Co-occurrances between top 50 MeSH terms")
JupyterPlot(p)
Out[12]:
In [13]:
using ARules
using DataTables
In [14]:
# Remove serch MeSH (set column to 0). Rules of lenght 2 with search term correspond to histogram,
# since it is in every transaction
mh_occ = convert(BitArray{2}, mesh_frequencies.array)
mh_col = mesh_frequencies.dicts[2][mh]
mh_occ[:, mh_col] = zeros(size(mh_occ,1))
mh_lkup = convert(DataStructures.OrderedDict{String,Int16}, mesh_frequencies.dicts[2])
@time mh_rules = apriori(mh_occ, supp = 0.001, conf = 0.1, maxlen = 9)
#Pretty print of rules
mh_lkup = Dict(zip(values(mesh_frequencies.dicts[2]), keys(mesh_frequencies.dicts[2])))
rules_dt= ARules.rules_to_datatable(mh_rules, mh_lkup, join_str = " | ");
In [15]:
println(head(rules_dt))
println("Found ", size(rules_dt, 1), " rules")
In [ ]:
supp_int = round(Int, 0.001 * size(mh_occ, 1))
@time root = frequent_item_tree(mh_occ, supp_int, 9);
supp_lkup = gen_support_dict(root, size(mh_occ, 1))
item_lkup = mesh_frequencies.dicts[2]
item_lkup_t = Dict(zip(values(item_lkup), keys(item_lkup)))
freq = ARules.suppdict_to_datatable(supp_lkup, item_lkup_t);
In [ ]:
function fill_sankey_data!(node, sources, targets, vals)
if length(node.item_ids) >1
push!(sources, node.item_ids[end-1]-1)
push!(targets, node.item_ids[end]-1)
push!(vals, node.supp)
end
if has_children(node)
for nd in node.children
fill_sankey_data!(nd, sources, targets, vals)
end
end
end
In [ ]:
sources = []
targets = []
vals = []
fill_sankey_data!(root, sources, targets, vals)
In [ ]:
trace=sankey(orientation="h",
node = attr(domain=attr(x=[0,1], y=[0,1]), pad=1/size(epilepsy_occ, 1), thickness=1/size(epilepsy_occ, 1), line = attr(color="black", width= 0.5),
label=mesh_names),
link = attr(source=sources, target=targets, value = vals))
plot([trace])
In [ ]: