<jupyter_start><jupyter_text>This notebook is for exploring completeness of VirusSeq Data Portal Metadata Last Update: March 16 2022  Import libraries<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
from copy import deepcopy
from datetime import datetime<jupyter_output><empty_output><jupyter_text>Read the sample data and see what it looks like<jupyter_code>seqs = pd.read_csv("data/files-archive-nov25.tsv", sep = "\t")
seqs.head()<jupyter_output><empty_output><jupyter_text>Print out all column names<jupyter_code>print(seqs.columns)<jupyter_output>Index(['study_id', 'specimen collector sample ID', 'sample collected by',
       'sequence submitted by', 'sample collection date',
       'sample collection date null reason', 'geo_loc_name (country)',
       'geo_loc_name (state/province/territory)', 'organism', 'isolate',
       'fasta header name', 'purpose of sampling',
       'purpose of sampling details', 'anatomical material', 'anatomical part',
       'body product', 'environmental material', 'environmental site',
       'collection device', 'collection method', 'host (scientific name)',
       'host disease', 'host age', 'host age null reason', 'host age unit',
       'host age bin', 'host gender', 'purpose of sequencing',
       'purpose of sequencing details', 'sequencing instrument',
       'sequencing protocol', 'raw sequence data processing method',
       'dehosting method', 'consensus sequence software name',
       'consensus sequence software version', 'breadth of coverage value',
       'depth of coverage value', 'r[...]<jupyter_text>Remove the columns I'm not going to examine here<jupyter_code>seqs.drop(columns=['study_id', 'specimen collector sample ID', 'sample collected by',
       'sequence submitted by',
       'sample collection date null reason', 'geo_loc_name (country)',
       'organism', 'isolate',
       'fasta header name',
       'purpose of sampling details', 'anatomical material', 'anatomical part',
       'body product', 'environmental material', 'environmental site',
       'collection device', 'collection method', 'host (scientific name)',
       'host disease','host age null reason', 'host age unit',
       'purpose of sequencing details',
       'sequencing protocol', 'raw sequence data processing method',
       'dehosting method', 'consensus sequence software name',
       'consensus sequence software version', 'reference genome accession',
       'bioinformatics protocol', 'gene name', 
       'diagnostic pcr Ct value null reason'], inplace = True)
seqs.head()<jupyter_output><empty_output><jupyter_text>Let's replace the many forms of 'Not Provided' or NaN with np.nan and convert collection date to a date time object<jupyter_code>seqs.replace({'Not Provided':np.nan, 'NaN':np.nan, 'Not Applicable':np.nan, 'Missing': np.nan}, inplace = True)
seqs['sample collection date'] = pd.to_datetime(seqs['sample collection date'])
seqs.head()<jupyter_output><empty_output><jupyter_text>What percent of the data is missing?<jupyter_code>seqs.isna().sum()/len(seqs)*100<jupyter_output><empty_output><jupyter_text>All samples have associated collection dates! Let's try grouping data by province to see if there's a pattern to what is missing.<jupyter_code>seqs.columns = ['collection_date', 'province', 'sampling_purpose', 'age', 'age_bin', 'gender', 'sequencing_purpose', 'sequencing_instrument',
'coverage_breadth', 'coverage_depth', 'Ct', 'gisaid_accession']
counts = seqs.groupby('province').count()
totals = counts['collection_date']
seqs.groupby('province').count().rsub(seqs.groupby('province').size(), axis = 0).divide(totals, axis = 0)<jupyter_output><empty_output><jupyter_text>Some of the provinces are missing very little data overall (NL, NS, MB), while those that are missing lots of data seem to have data missing from common fields (age, age bin, gender). Some missing fields are unique to provinces (only BC is missing sampling purpose).  How many samples are there from each province?<jupyter_code>seqs.groupby('province').size()<jupyter_output><empty_output><jupyter_text>How many complete values are there from each province in each field?<jupyter_code>counts<jupyter_output><empty_output><jupyter_text>What proportion of each field is complete for each province?<jupyter_code>seqs_comp = deepcopy(pd.DataFrame(seqs))
seqs_comp = seqs_comp.groupby('province').count()

seqs_comp.divide(totals, axis = 0)<jupyter_output><empty_output><jupyter_text>Save .csv file<jupyter_code>seq_proportions = seqs_comp.divide(totals, axis = 0)
seq_proportions.to_csv('data/data_portal_proportions.csv', index = True)<jupyter_output><empty_output><jupyter_text>Let's explore the column values in more detail<jupyter_code>print(seqs.columns)
cols = ['province', 'sampling_purpose', 'age', 'age_bin',
       'gender', 'sequencing_purpose', 'sequencing_instrument',
       'coverage_breadth', 'coverage_depth', 'Ct', 'gisaid_accession']
print(seqs.sampling_purpose.unique())
print(seqs.age_bin.unique())
print(seqs.gender.unique())
print(seqs.sequencing_purpose.unique())
print(seqs.sequencing_instrument.unique())
print(seqs.Ct.unique())<jupyter_output>[  nan 14.2  12.4  ... 28.83 14.41 29.28]<jupyter_text>Let's figure out how many samples are in the portal from each province by collection month<jupyter_code>seqs['year_month'] = pd.to_datetime(seqs['collection_date']).dt.to_period('M')
seqs.head()
months_counts = pd.DataFrame(seqs.groupby(['province', 'year_month']).size())
totals = months_counts[0]
totals<jupyter_output><empty_output><jupyter_text>Let's see if data completeness varies by collection date<jupyter_code>seqs_comp_month = deepcopy(pd.DataFrame(seqs))
seqs_comp_month = seqs_comp_month.groupby(['province', 'year_month']).count()
seqs_comp_month.divide(totals, axis = 0)<jupyter_output><empty_output><jupyter_text>Save to csv file<jupyter_code>seqs_comp_month.to_csv('data/data_portal_proportions_month.csv', index = True)<jupyter_output><empty_output><jupyter_text>As a preliminary sanity check, let's see whether genome coverage breadth correlates inversely with Ct value<jupyter_code>seqs_model = deepcopy(seqs)
seqs_model.drop(columns = ['sampling_purpose', 'sequencing_purpose', 'gisaid_accession'], inplace =True)
seqs_model<jupyter_output><empty_output><jupyter_text>Remove percentage signs and 'x' characters from coverage columns<jupyter_code>print(seqs_model.dtypes)
seqs_model['coverage_breadth'] = seqs_model['coverage_breadth'].astype('str')
seqs_model['coverage_depth'] = seqs_model['coverage_depth'].astype('str')
seqs_model['coverage_breadth'] = seqs_model['coverage_breadth'].map(lambda x: x.rstrip('%'))
seqs_model['coverage_depth'] = seqs_model['coverage_depth'].map(lambda x: x.rstrip('Xx'))
seqs_model<jupyter_output>collection_date          datetime64[ns]
province                         object
age                             float64
age_bin                          object
gender                           object
sequencing_instrument            object
coverage_breadth                 object
coverage_depth                   object
Ct                              float64
year_month                    period[M]
dtype: object<jupyter_text>Drop all rows where coverage info is absent, change data type to float<jupyter_code>seqs_model.dropna(subset=['coverage_breadth', 'coverage_depth', 'Ct'], inplace = True)
seqs_model['coverage_breadth'] = seqs_model['coverage_breadth'].astype('float64')
seqs_model['coverage_depth'] = seqs_model['coverage_depth'].astype('float64')
seqs_model<jupyter_output><empty_output><jupyter_text>Subset dataset to Ct values greater than 1 and less than 45<jupyter_code>seqs_model = seqs_model[seqs_model['Ct'] < 45]
seqs_model = seqs_model[seqs_model['Ct'] > 1]<jupyter_output><empty_output><jupyter_text>Down-sample dataset so that we can visualize well with a jointplot<jupyter_code>perc = 0.05

df_sub = seqs_model.sample(n=int(seqs_model.shape[0]*perc))
sns.set_theme(style = "darkgrid")
sns.set(rc={'figure.figsize':(20,20)})
sns.jointplot(x = 'Ct', y='coverage_breadth', kind = 'reg', data = df_sub)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12050116693280924, 'fraction_numerical': 0.007247266920525734, 'mean_word_length': 5.8076923076923075, 'pattern_counts': {'":': 0, '<': 79, '<?xml version=': 0, '>': 79, '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': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '622817', 'license': 'MIT', 'n_tokens_mistral': 2676, 'n_tokens_neox': 2557, 'n_words': 625, 'path': 'data_portal_contents.ipynb', 'repo_name': 'eringill/gisaid_uploads', 'size': '8141'}