<jupyter_start><jupyter_text>analyze_metadata_tableThis notebook builds on the previously produced data table, performs QC, and adds meta-features for visualization.<jupyter_code># import required modules
import numpy as np
import pandas as pd
from collections import Counter
import statistics

import datetime
from datetime import date
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib import rc<jupyter_output><empty_output><jupyter_text>Import data<jupyter_code># import table from create_metadata_table.ipynb
allFactors = pd.read_csv('../data_tables/metadataMatrix_raw.csv', low_memory=False)<jupyter_output><empty_output><jupyter_text>Conduct analysis on this table, producing new metrics<jupyter_code># add a column with the total use count of each converted_accession ID

# first, create a new dataframe with no double-counted accessions
# justAcc includes all unique pairings of PMC ID and Series/Study accession
justAcc = allFactors.loc[:, ['pmc_ID', 'converted_accession']]
justAcc = justAcc.drop_duplicates()

rc = Counter(justAcc['converted_accession'])
reuse_counts = pd.DataFrame.from_dict(rc, orient = 'index').reset_index()
reuse_counts.columns = ['converted_accession', 'total_use_count']
# reuse_counts
# merge reuse counts onto table
allFactors = pd.merge(allFactors, reuse_counts, how = 'left', on = 'converted_accession')
# allFactors
# normalize reuse count by amount of time public

# convert repository upload dates and today's date to datetime objects
allFactors.loc[:,'repository_date'] = allFactors['repository_date'].replace({' ':''})
fmt = '%Y-%m-%d'
allFactors.loc[:,'repository_date'] = pd.to_datetime(allFactors['repository_date'], format = fmt)
allFactors['today'] = date.today().strftime(fmt)
allFactors['today'] = pd.to_datetime(allFactors['today'], format = fmt)

# create timedelta objects representing amount of time public in YEARS
allFactors['time_public'] = allFactors['today'] - allFactors['repository_date']
allFactors['time_public'] = allFactors['time_public'].dt.days / 365

# normalize reuse count by years public
allFactors['reuse_count'] = (allFactors['total_use_count'] - 1)
allFactors['norm_reuse_count'] = allFactors['reuse_count'] / allFactors['time_public']
allFactors = allFactors.drop(labels = ['total_use_count'], axis = 1)
# allFactors
# Tag each paper as either a generator or a reuser

# take each converted accession along with the first date it appears in a paper
justPap = allFactors.loc[:, ['converted_accession', 'pmc_date']].drop_duplicates()
justPap = justPap.groupby('converted_accession', as_index = False).agg({'pmc_date': 'min'})
justPap = justPap.rename(columns = {'pmc_date':'min_pmc_date'})

# label all of the usages of these datasets on these dates with a 'G' (generators)
justPap.loc[:,'reuse_role'] = 'G'
# justPap
# merge these back onto the original table, matching dates back to papers
allFactors = pd.merge(allFactors, justPap, how = 'left', 
                      left_on = ['converted_accession', 'pmc_date'], 
                      right_on = ['converted_accession', 'min_pmc_date'])

# label every usage of a dataset that isn't 'G' (generator) as 'R' (reuser)
allFactors['allR'] = 'R'
allFactors['reuse_role'] = allFactors['reuse_role'].fillna(allFactors['allR'])

# drop redundant columns
allFactors = allFactors.drop(labels = ['min_pmc_date', 'allR'], axis = 1)

# allFactors
# perform QC on generator/reuser tags...

# cut out whitespace
allFactors.loc[:,'pmc_date'] = allFactors['pmc_date'].replace({' ':''})

# take slice of what we need from allFactors
# gens = allFactors[['pmc_ID', 'converted_accession', 'repository', 'pmc_date', 'repository_date', 'reuse_role']]
gens = allFactors.loc[allFactors['reuse_role'] == 'G', :]

# convert to date objects
fmt = '%Y-%m-%d'
gens.loc[:,'pmc_date'] = pd.to_datetime(gens['pmc_date'], format = fmt)
gens.loc[:,'repository_date'] = pd.to_datetime(gens['repository_date'], format = fmt)

# subtract columns of date objects, creating a column of timedelta objects
gens.loc[:,'pub_delay'] = gens['pmc_date'] - gens['repository_date']
gens = gens.dropna(subset = ['pub_delay'])

# convert timedeltas back to integers
gens.loc[:,'pub_delay'] = gens['pub_delay'].dt.days
gens = gens.sort_values(by = 'pub_delay', ascending = False)

# gens

# head = gens.head(n = 15)
# tail = gens.tail(n = 15)
# plot timedelta distribution
fig_dims = (18,12)
fig, ax = plt.subplots(figsize=fig_dims)
sns.distplot(gens['pub_delay'], ax = ax)
# select only papers in a reasonable timedelta window

DELAY_MIN = -270
DELAY_MAX = 270

gens_filtered = gens.loc[gens['pub_delay'] <= DELAY_MAX, :].loc[gens['pub_delay'] >= DELAY_MIN, :]
# gens_filtered
# select random papers inside/outside threshold for QC analysis
import random

# gens_filtered are inside threshold already
# get papers outside threshold, but within 365 days

DELAY_EDGE = 365

gens_aboveThresh = gens.loc[gens['pub_delay'] >= DELAY_MAX, :].loc[gens['pub_delay'] <= DELAY_EDGE]
gens_aboveThresh = gens_aboveThresh.loc[gens_aboveThresh['repository'] == "GEO", :]
gens_belowThresh = gens.loc[gens['pub_delay'] <= DELAY_MIN, :].loc[gens['pub_delay'] >= -DELAY_EDGE]
gens_belowThresh = gens_belowThresh.loc[gens_belowThresh['repository'] == "GEO", :]
gens_insideThresh = gens_filtered.loc[gens_filtered['repository'] == "GEO", :]

num_papers = 50
half_papers = int(num_papers / 2)

# take random samples inside, above, and below threshold
rand_in = pd.DataFrame(random.sample(gens_insideThresh['pmc_ID'].tolist(), num_papers))
rand_in.columns = ['pmc_ID']
rand_ab = pd.DataFrame(random.sample(gens_aboveThresh['pmc_ID'].tolist(), half_papers))
rand_ab.columns = ['pmc_ID']
rand_be = pd.DataFrame(random.sample(gens_belowThresh['pmc_ID'].tolist(), half_papers))
rand_be.columns = ['pmc_ID']

# match accessions to papers to make manual work easier
r_in = pd.merge(rand_in, gens_insideThresh[['pmc_ID', 'accession']], on = 'pmc_ID', how = 'left').drop_duplicates()
r_in['real_introducer'] = np.NaN
r_ab = pd.merge(rand_ab, gens_aboveThresh[['pmc_ID', 'accession']], on = 'pmc_ID', how = 'left').drop_duplicates()
r_be = pd.merge(rand_be, gens_belowThresh[['pmc_ID', 'accession']], on = 'pmc_ID', how = 'left').drop_duplicates()
r_out = pd.concat([r_ab, r_be])
r_out['real_introducer'] = np.NaN

# save incomplete files to csv
r_in.to_csv('../data_tables/randomIntroducers_INSIDE.csv', index = False)
r_out.to_csv('../data_tables/randomIntroducers_OUTSIDE.csv', index = False)
# not included in the script: have someone else (very nicely) manually 
# validate each of these randomly selected paper-dataset relationships

# import the manually produced data (now complete files)
rinf = pd.read_csv("../data_tables/randomIntroducers_INSIDE_270days_filledin.csv").dropna()
routf = pd.read_csv("../data_tables/randomIntroducers_OUTSIDE_270days_filledin.csv").dropna()

# create a true positive column
# true positive means PMC_ID == REAL_INTRODUCER
rinf["TP"] = np.where(rinf["pmc_ID"] == rinf["real_introducer"], True, False)
routf["TN"] = np.where(routf["pmc_ID"] != routf["real_introducer"], True, False)

# reduce the tables to just PMC_ID and TP/TN
rinf = rinf[["pmc_ID", "TP"]].drop_duplicates()
routf = routf[["pmc_ID", "TN"]].drop_duplicates()

tp = statistics.mean(rinf["TP"])
tn = statistics.mean(routf["TN"])

print("True Positive rate (inside threshold): " + str(tp))
print("True Negative rate (outside threshold, inside year): " + str(tn))
not_gens = allFactors.loc[allFactors['reuse_role'] != 'G', :]
filtered_matrix = pd.concat(objs = [gens_filtered, not_gens], ignore_index = True, sort = False)
# filtered_matrix
lsc = Counter(filtered_matrix.loc[filtered_matrix['repository'] == "SRA", ]["library_strategy"])
lsc_df = pd.DataFrame.from_dict(lsc, orient = 'index').reset_index()
sraLS = list(lsc_df["index"])
sraLS.remove("RNA-Seq")
# sraLS
# add a column for "class"... separating GEO/SRA distinction into:
# (1) GEO gene expression
# (2) RNA-seq (GEO and SRA)
# (3) SRA other omics technologies

key = {"repository": ["GEO", "GEO", "SRA"] + ["SRA"]*(len(sraLS)), 
       "library_strategy": ["Expression_Array", "RNA-Seq", "RNA-Seq"] + sraLS,
      "class": ["GEO_expr", "RNA-Seq", "RNA-Seq"] + ["SRA_other"]*(len(sraLS))}
key = pd.DataFrame(key)
key
filtered_matrix = pd.merge(filtered_matrix, key, how = 'left', on = ['repository', 'library_strategy'])
filtered_matrix
# what are the SRA library strategies?
lc = Counter(filtered_matrix.loc[filtered_matrix["repository"] == "SRA", :]["library_strategy"])
libCntr = pd.DataFrame.from_dict(lc, orient = "index")
libCntr.to_csv("../data_tables/SRAlibraryStrategies.csv", header = False)
# how much data coverage do we get when using the class column?
Counter(filtered_matrix['class'])
filtered_matrix.to_csv('../data_tables/metadataMatrix_filtered.csv', index = False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14624773960216997, 'fraction_numerical': 0.004294755877034358, 'mean_word_length': 6.442388561816653, 'pattern_counts': {'":': 3, '<': 16, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '373370', 'license': 'MIT', 'n_tokens_mistral': 3159, 'n_tokens_neox': 3012, 'n_words': 718, 'path': 'jupyter_notebooks/analyze_metadata_table.ipynb', 'repo_name': 'smangul1/data_reusability', 'size': '8848'}