from Bio import Entrez
import requests
import time
import pandas as pd
import networkx as nx
import matplotlib.pyplot as plt
from tqdm import tqdm

Entrez.email = "youremail@dot.com"  # Always tell NCBI who you are
url = "http://bern2.korea.ac.kr/pubmed"  # Remember 100 requests for 100 seconds


def get_pmid(query, num_list):
    """
    num_list should be below 10000
    """
    handle = Entrez.esearch(
        db="pubmed", term=query, retmax=num_list, sort="pub+date", retmode="xml"
    )
    records = Entrez.read(handle)
    return records["IdList"]


def query_pmid(pmids, url="http://bern2.korea.ac.kr/pubmed"):
    try:
        return requests.get(url + "/" + ",".join(pmids)).json()
    except:
        pass


def make_table(json_list):
    try:
        df = pd.json_normalize(json_list, record_path=["annotations"], meta=["pmid"])
        refine_df = df[["mention", "obj", "prob", "pmid"]]
        # use pd.astype() function for save memory
        refine_df = refine_df.astype(
            {
                "prob": "float16",
                "pmid": "int32",
            }
        )
        return refine_df
    except:
        return []


def get_bern2(pmids):
    temp_df = pd.DataFrame()
    for i in tqdm(pmids, unit="pmid"):
        # print(query_pmid(i))
        new_df = make_table(query_pmid([i]))  # only list
        time.sleep(1)  # delay for BERN2 API
        try:
            temp_df = pd.concat([temp_df, new_df], ignore_index=True)
        except:
            pass
    return temp_df

entrez_query = input("What do want to search in PubMed? ")
pmids = get_pmid(entrez_query, "10")  # string is not a mistake
num_pmids = len(pmids)
print(f"Number of search results is: {num_pmids}!")

df = get_bern2(pmids)
print(f"Shape of result table is: {df.shape}.")

df.tail()

df.info()

counts = df["mention"].value_counts()
num_nodes_to_inspect = 10
counts[:num_nodes_to_inspect].plot(kind="barh").invert_yaxis()

df["mention"] = df["mention"].replace(["CXCL8"], "IL-8")
df["mention"] = df["mention"].replace(["interleukin-8"], "IL-8")
df["mention"] = df["mention"].replace(["Interleukin-8"], "IL-8")
df["mention"] = df["mention"].replace(["interleukin 8"], "IL-8")
df["mention"] = df["mention"].replace(["IL-8 receptors"], "IL-8 receptor")
df["mention"] = df["mention"].replace(["calcium"], "Ca2+")
df["mention"] = df["mention"].replace(["PMNs"], "PMN")
df["mention"] = df["mention"].replace(["murine"], "mice")
df["mention"] = df["mention"].replace(["mouse"], "mice")
df["mention"] = df["mention"].replace(["patients"], "human")

df.drop(df[df["mention"] == "CXCR1"].index, inplace=True)  # Alternative name of CD128
df.drop(
    df[df["mention"] == "IL-8 receptor"].index, inplace=True
)  # Alternative name of CD128
df.drop(df[df["mention"] == "N"].index, inplace=True)
df.drop(df[df["mention"] == "C"].index, inplace=True)
df.drop(df[df["mention"] == "IL-8R"].index, inplace=True)

counts = df["mention"].value_counts()
num_nodes_to_inspect = 10
counts[:num_nodes_to_inspect].plot(kind="barh").invert_yaxis()

# df.to_csv(f"./output/{entrez_query}.csv") # CSV 파일 저장

G = nx.from_pandas_edgelist(df, source="pmid", target="mention", edge_attr="obj")
len(G.nodes())

degrees = dict(nx.degree(G))
nx.set_node_attributes(G, name="degree", values=degrees)
degree_df = pd.DataFrame(G.nodes(data="degree"), columns=["node", "degree"])
degree_df = degree_df.sort_values(by="degree", ascending=False)
degree_df.head(10)

degree_df.describe()

# remove low-degree nodes
low_degree = [n for n, d in G.degree() if d < 2]
G.remove_nodes_from(low_degree)

# Specify figure size
plt.figure(figsize=(20, 15))

# Compute node position using the default spring_layout
node_position = nx.spring_layout(G)
nx.draw(
    G,
    node_position,
    node_color="#F4ABAA",
    edge_color="gainsboro",
    with_labels=True,
    alpha=0.6,
)
plt.show()

G = nx.from_pandas_edgelist(df, source="pmid", target="mention", edge_attr="obj")

# largest connected component
components = nx.connected_components(G)
largest_component = max(components, key=len)
H = G.subgraph(largest_component)

# compute centrality
centrality = nx.betweenness_centrality(H, k=10, endpoints=True)

# compute community structure
lpc = nx.community.label_propagation_communities(H)
community_index = {n: i for i, com in enumerate(lpc) for n in com}

#### draw graph ####
fig, ax = plt.subplots(figsize=(20, 20))
pos = nx.spring_layout(H, k=0.15, seed=42)
node_color = [community_index[n] for n in H]
node_size = [v * 20000 for v in centrality.values()]
nx.draw_networkx(
    H,
    pos=pos,
    with_labels=True,
    node_color=node_color,
    node_size=node_size,
    edge_color="gainsboro",
    alpha=0.6,
)

# Title/legend
font = {"color": "k", "fontweight": "bold", "fontsize": 20}
ax.set_title(f"{entrez_query} network", font)


ax.text(
    0.80,
    0.10,
    "node color = community structure",
    horizontalalignment="center",
    transform=ax.transAxes,
    fontdict=font,
)
ax.text(
    0.80,
    0.08,
    "node size = betweeness centrality",
    horizontalalignment="center",
    transform=ax.transAxes,
    fontdict=font,
)

# Resize figure for label readibility
ax.margins(0.1, 0.05)
fig.tight_layout()
plt.axis("off")
plt.show()

%matplotlib inline

from pathlib import Path

import pandas as pd
import numpy as np
from scipy.stats import trim_mean
from statsmodels import robust
import wquantiles

import seaborn as sns
import matplotlib.pylab as plt

try:
    import common

    DATA = common.dataDirectory()
except ImportError:
    DATA = Path().resolve() / "data"

AIRLINE_STATS_CSV = DATA / "airline_stats.csv"
KC_TAX_CSV = DATA / "kc_tax.csv.gz"
LC_LOANS_CSV = DATA / "lc_loans.csv"
AIRPORT_DELAYS_CSV = DATA / "dfw_airline.csv"
SP500_DATA_CSV = DATA / "sp500_data.csv.gz"
SP500_SECTORS_CSV = DATA / "sp500_sectors.csv"
STATE_CSV = DATA / "state.csv"

state = pd.read_csv(STATE_CSV)
state.head(8)

state = pd.read_csv(STATE_CSV)
print(state["Population"].mean())

print(trim_mean(state["Population"], 0.1))

print(state["Population"].median())

print(np.average(state["Murder.Rate"], weights=state["Population"]))

print(wquantiles.median(state["Murder.Rate"], weights=state["Population"]))

state.head(8)

print(state["Population"].std())

print(robust.scale.mad(state["Population"]))
# print(abs(state['Population'] - state['Population'].median()).median() / 0.6744897501960817)

print(state["Murder.Rate"].quantile([0.05, 0.25, 0.5, 0.75, 0.95]))

# Table 1.4
percentages = [0.05, 0.25, 0.5, 0.75, 0.95]
df = pd.DataFrame(state["Murder.Rate"].quantile(percentages))
df.index = [f"{p * 100}%" for p in percentages]
print(df.transpose())

ax = (state["Population"] / 1_000_000).plot.box(figsize=(3, 4))
ax.set_ylabel("Population (millions)")

plt.tight_layout()
plt.show()

binnedPopulation = pd.cut(state["Population"], 10)
print(binnedPopulation.value_counts())

각각의 주가 어느 구간에 속하는지 좀 더 보기좋게 만들어봅니다.

# Table 1.5
binnedPopulation.name = "binnedPopulation"
df = pd.concat([state, binnedPopulation], axis=1)
df = df.sort_values(by="Population")

groups = []
for group, subset in df.groupby(by="binnedPopulation"):
    groups.append(
        {
            "BinRange": group,
            "Count": len(subset),
            "States": ",".join(subset.Abbreviation),
        }
    )
print(pd.DataFrame(groups))

ax = (state["Population"] / 1_000_000).plot.hist(figsize=(4, 4))
ax.set_xlabel("Population (millions)")

plt.tight_layout()
plt.show()

ax = state["Murder.Rate"].plot.hist(
    density=True, xlim=[0, 12], bins=range(1, 12), figsize=(4, 4)
)
state["Murder.Rate"].plot.density(ax=ax)
ax.set_xlabel("Murder Rate (per 100,000)")

plt.tight_layout()
plt.show()

# Table 1-6
dfw = pd.read_csv(AIRPORT_DELAYS_CSV)
print(100 * dfw / dfw.values.sum())

ax = dfw.transpose().plot.bar(figsize=(4, 4), legend=False)
ax.set_xlabel("Cause of delay")
ax.set_ylabel("Count")

plt.tight_layout()
plt.show()

sp500_sym = pd.read_csv(SP500_SECTORS_CSV)
sp500_px = pd.read_csv(SP500_DATA_CSV, index_col=0)

etfs = sp500_px.loc[
    sp500_px.index > "2012-07-01", sp500_sym[sp500_sym["sector"] == "etf"]["symbol"]
]
print(etfs.head())

fig, ax = plt.subplots(figsize=(5, 4))
ax = sns.heatmap(
    etfs.corr(),
    vmin=-1,
    vmax=1,
    cmap=sns.diverging_palette(20, 220, as_cmap=True),
    ax=ax,
)

plt.tight_layout()
plt.show()

ax = telecom.plot.scatter(x="T", y="VZ", figsize=(4, 4), marker="$\u25EF$", alpha=0.5)
# $\u25EF$는 속이 빈 원을 뜻한다.
ax.set_xlabel("ATT (T)")
ax.set_ylabel("Verizon (VZ)")
ax.axhline(0, color="grey", lw=1)
print(ax.axvline(0, color="grey", lw=1))

kc_tax = pd.read_csv(KC_TAX_CSV)
kc_tax0 = kc_tax.loc[
    (kc_tax.TaxAssessedValue < 750000)
    & (kc_tax.SqFtTotLiving > 100)
    & (kc_tax.SqFtTotLiving < 3500),
    :,
]
print(kc_tax0.shape)

ax = kc_tax0.plot.hexbin(
    x="SqFtTotLiving", y="TaxAssessedValue", gridsize=30, sharex=False, figsize=(5, 4)
)
ax.set_xlabel("Finished Square Feet")
ax.set_ylabel("Tax Assessed Value")

plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(4, 4))
sns.kdeplot(data=kc_tax0.sample(10000), x="SqFtTotLiving", y="TaxAssessedValue", ax=ax)
ax.set_xlabel("Finished Square Feet")
ax.set_ylabel("Tax Assessed Value")

plt.tight_layout()
plt.show()

lc_loans = pd.read_csv(LC_LOANS_CSV)

df = crosstab.copy().loc["A":"G", :]
df.loc[:, "Charged Off":"Late"] = df.loc[:, "Charged Off":"Late"].div(df["All"], axis=0)
df["All"] = df["All"] / sum(df["All"])
perc_crosstab = df
print(perc_crosstab)

airline_stats = pd.read_csv(AIRLINE_STATS_CSV)
airline_stats.head()
ax = airline_stats.boxplot(by="airline", column="pct_carrier_delay", figsize=(5, 5))
ax.set_xlabel("")
ax.set_ylabel("Daily % of Delayed Flights")
plt.suptitle("")

plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(5, 5))
sns.violinplot(
    data=airline_stats,
    x="airline",
    y="pct_carrier_delay",
    ax=ax,
    inner="quartile",
    color="white",
)
ax.set_xlabel("")
ax.set_ylabel("Daily % of Delayed Flights")

plt.tight_layout()
plt.show()

zip_codes = [98188, 98105, 98108, 98126]
kc_tax_zip = kc_tax0.loc[kc_tax0.ZipCode.isin(zip_codes), :]
kc_tax_zip


def hexbin(x, y, color, **kwargs):
    cmap = sns.light_palette(color, as_cmap=True)
    plt.hexbin(x, y, gridsize=25, cmap=cmap, **kwargs)


g = sns.FacetGrid(kc_tax_zip, col="ZipCode", col_wrap=2)
g.map(hexbin, "SqFtTotLiving", "TaxAssessedValue", extent=[0, 3500, 0, 700000])
g.set_axis_labels("Finished Square Feet", "Tax Assessed Value")
g.set_titles("Zip code {col_name:.0f}")

plt.tight_layout()
plt.show()