# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Preprocessing data from the COVID-19 multiomics project
#
# This notebook is provided as a similar case study for our pipeline. [It is conceptually similar to our previous integrative analysis of proteome and translatome data, but contains a larger quantity of omics data and samples](https://gitlab.com/tyagilab/sars-cov-2/-/blob/master/README.md).
#
# Full details of the original study can be found in the original publication: [https://dx.doi.org/10.1101%2F2020.07.17.20156513](https://dx.doi.org/10.1101%2F2020.07.17.20156513)
#
# [Authors provided a sql database with accession number `MSV000085703`](ftp://massive.ucsd.edu/MSV000085703/other/Covid-19%20Study%20DB.sqlite) from which multiomics data measurements were extracted.
#
# ## Summary of the original study
#
# *This information is taken from the original study by [Overmyer et al 2020](https://dx.doi.org/10.1101%2F2020.07.17.20156513)*
#
# The authors intended to investigate (a) biological insight into the host's response to SARS-CoV-2, and (b) pathways influencing its severity. Authors integrated molecule measurements with clinical attributes and identified features associated with (1) COVID-19 status and (2) `HFD-45` (hospital free days, this is a composite metric).
#
# The authors performed multiple experiments to obtain their results (summarised below). We focus specifically on the experiment of `covid severity`. As it has a relatively straightforward experimental design, we illustrate the usage of our pipeline on this case study.
#
# ### Covid state
#
# To identify differentially abundant molecules, authors performed **ANOVA and log-likelihood ratio tests** to discover: 
# - 2,537 leukocyte transcripts, 
# - 146 plasma proteins, 
# - 168 plasma lipids,
# - 13 plasma metabolites
# associated with COVID-19 status (Table S1). 
#
# To discover enriched biological processes associated with the differing biomolecules, the authors used **GO and molecular class enrichment analysis** (Tables S2A and S2B). The authors showed that these included: 
# - mitotic cell cycle, 
# - phagocytosis recognition, 
# - positive recognition of B cell activation, 
# - complement activation (classical pathway), 
# - innate immune response.
#
# ### Covid severity
#
# To discover biomolecules associated with severity, authors used **univariate regression** of HFD-45 against abundance of each biomolecule. Authors accounted for sex and age and found disease associated molcules (Table S1):
# - 6,202 transcripts, 
# - 189 plasma proteins, 
# - 218 plasma lipids,  
# - 35 plasma small molecules
#
# To further refine these authors performed **multivariate linear regression** on HFD-45 using the elastic net penalty (Zou and Hastie, 2005) as predictive features for HFD-45 (Table S1):
# - 497 transcripts, 
# - 382 proteins, 
# - 140 lipids, 
# - 60 metabolites 
#
# ### Combined information
#
# To generate a list of 219 features that were most significantly associated with COVID-19 status and severity (Figure 2C; Table S1), the authors combined: 
# 1. significance with COVID-19 status, 
# 2. significance with HFD-45, 
# 3. elastic net feature selection
#
# ## Summary of our analysis
#
# We integrated lipidomics, metabolomics, proteomics and transcriptomics data for 100 samples. Two classes were included: `less severe` vs `more severe` covid states, represented by the variable `HFD-45`. Classes are balanced but age and sex of patients are not consistent.
#
# This jupyter notebook describes the steps taken to download and parse the input data as well as metadata, resulting in matrices of continuous values suitable for input into our pipeline. The appendix contains suggestions for handling cases where experimental design is less straightforward.

import os
import matplotlib.pyplot as plt
import pandas as pd
from functools import reduce
from scipy import stats
from scipy.stats import ttest_ind
from statsmodels.stats.multitest import multipletests
# %matplotlib inline

# Tables in sqlite database:
#
# ```
# biomolecules                   
# deidentified_patient_metadata  
# lipidomics_measurements        proteomics_measurements
# lipidomics_runs                proteomics_runs
# metabolomics_measurements      pvalues
# metabolomics_runs              rawfiles
# metadata                       transcriptomics_measurements
# omes                           transcriptomics_runs
# ```
#
# `.csv` files were extracted from the publicly accessible `sqlite` database. The following set of commands were repeated for each of the tables listed above.
#
# ```
# .file 'Covid-19 Study DB.sqlite'
# .header on
# .mode csv
#
# .output ${TABLE}.csv
# select * from ${TABLE};
# ```
#
# `md5` sums are shown below for reproducibility:

# !for i in $(find ../../data/MSV000085703 -name *csv.gz | sort); do md5 $i; done
# !find ../../data/MSV000085703 -name *csv.gz -exec gzip -d {} \;

files = [x for x in os.listdir("../../data/MSV000085703") if x.endswith("csv")]
data = [pd.read_csv("/".join(["../../data/MSV000085703", x]), sep=",") for x in files]
data = dict(zip(files, data))
data.keys()


# Data is spread across individual `sql` tables. We remap the numerical codes to biologically meaningful identifiers. Database schema is available in Supp Fig S3 of the original manuscript.

# +
def biomolecules_to_omes(data):
    gene = pd.DataFrame([[5, "Gene"]], columns=["omics_id", "omics_name"])
    omes = pd.concat([data["omes.csv"], gene])
    biomolecules = data["biomolecules.csv"]
    biomolecules = biomolecules.merge(omes, left_on="omics_id", right_on="omics_id", how="outer")
    return biomolecules

biomolecules = biomolecules_to_omes(data)
keep = biomolecules[biomolecules["keep"] == "1"]["biomolecule_id"].tolist()
biomolecules
# -

# ## Selection of Outcome Measure
#
# The authors constructed a composite variable `HFD-45` (hospital-free days at day 45). This assigns zero value to patients requiring admission longer than 45 days or who die during the admission, and progressively more free days depending on the hospitalization length. A lower value indicates higher severity.
#
# The variable is intended to:
#
# 1. be able to combine severity of disease with mortality in one single metric; 
# 2. be amenable to both ICU and medical floor populations; 
# 3. use a timeframe that accounts for the fact that COVID-19 patients with respiratory failure require longer hospitalizations compared with non-COVID-19 individuals (Wang et al., 2020a, 2020b); 
# 4. consider that COVID-19 causes a linear disease’s deterioration pattern that transition from mild respiratory compromise to respiratory failure, followed by respiratory distress requiring mechanical ventilatory support and eventually death.
#
# ## Examining the raw data

patient_metadata = data["deidentified_patient_metadata.csv"]
raw_files = data["rawfiles.csv"]
patient_files = raw_files.merge(patient_metadata, left_on="sample_id", right_on="sample_id", how="left")
patient_files.rename(columns={"keep": "keep_patient"}, inplace=True)
patient_files


# The raw data is assigned to individual categories based on their omics category (lipidome, metabolome, proteome, transcriptome).

# +
def biomolecules_to_data(biomolecules, data):
    # biomolecules = biomolecules[biomolecules["keep"] == "1"]
    return data.merge(biomolecules, left_on="biomolecule_id", right_on="biomolecule_id", how="left")

def measurements_to_runs(data, runs):
    return data.merge(runs, left_on="replicate_id", right_on="replicate_id", how="left")

def runs_to_patient(runs, patient):
    return runs.merge(patient, left_on="rawfile_id", right_on="rawfile_id", how="left")

omics = {"L": data['lipidomics_measurements.csv'],
         "M": data['metabolomics_measurements.csv'],
         "P": data['proteomics_measurements.csv'],
         "T": data['transcriptomics_measurements.csv'],}

runs = {"L": data['lipidomics_runs.csv'],
        "M": data['metabolomics_runs.csv'],
        "P": data['proteomics_runs.csv'],
        "T": data['transcriptomics_runs.csv'],}

omics = [measurements_to_runs(omics[i], runs[i]) for i in omics.keys()]
omics = dict(zip(runs.keys(), omics))
omics = [runs_to_patient(omics[i], patient_files) for i in omics.keys()]
omics = dict(zip(runs.keys(), omics))
data_omics = [biomolecules_to_data(biomolecules, omics[i]) for i in omics.keys()]
patient_status = "Albany_sampleID" # unique_identifier
# -

# ## Lipidomics
#
# ### Experiment notes
#
# "The LC–MS data were processed using Compound Discoverer 2.1 (Thermo Scientific) and LipiDex (Hutchins et al., 2018) (v. 1.1.0). All peaks between 1 min and 45 min retention time and 100 Da to 5000 Da MS1 precursor mass were grouped into distinct chromatographic profiles (i.e., compound groups) and aligned using a 10-ppm mass and 0.3 min retention time tolerance. Profiles not reaching a minimum peak intensity of 5x10ˆ5, a maximum peak-width of 0.75, a signal-to-noise (S/N) ratio of 3, and a 3-fold intensity increase over blanks were excluded from further processing. MS/MS spectra were searched against an in-silico generated lipid spectral library containing 35,000 unique molecular compositions representing 48 distinct lipid classes (LipiDex library “LipiDex_HCD_Formic”, with a full range of acyl-chains included). Spectral matches with a dot product score greater than 500 and a reverse dot product score greater than 700 were retained for further analysis, with a minimum 75% spectral purity for designating fatty acid composition. Removed from the data set were adducts, class IDs greater than 3.5 median absolute retention time deviation (M.A.D. RT) of each other, and features found in less than 3 files. Data were additionally searched with Compound Discoverer 3.1 with the discovery metabolomics nodes for additional spectral matching to mzCloud and mzVault libraries but retaining the feature group and peak picking settings as detailed for the Compound Discoverer 2.1 analysis." - [Overmyer et al 2020](https://dx.doi.org/10.1101%2F2020.07.17.20156513)
#
# ### Data notes
#
# Data is a matrix of continuous values. Column names correspond to lipid names and row names correspond to unique sample identifiers. Lipids with no known annotation are also recorded and assigned a unique identifier based on their mass and charge.

# +
l = data_omics[0][["unique_identifier", "standardized_name", "normalized_abundance", "keep"]]
l = l[l["keep"] == "1"]
l.drop("keep", axis=1, inplace=True)
l = l.pivot(index="unique_identifier", columns="standardized_name")

# this maps the sample to patient
l_id = data_omics[0][["unique_identifier", patient_status]]
l_id.drop_duplicates(inplace=True)
l_id.set_index("unique_identifier", inplace=True)

l = l.merge(l_id, left_index=True, right_index=True, how="left")
l.reset_index(inplace=True)
l.set_index(["unique_identifier", patient_status], inplace=True)
l.sort_index(level=patient_status, inplace=True)

# remove null and assign sample groups
l_class = l.reset_index()
l_class.dropna(inplace=True)

l_class["covid_state"] = l_class[patient_status]
l_class["covid_state"].replace(regex=r'^C.*', value='Covid', inplace=True)
l_class["covid_state"].replace(regex=r'^NC.*', value='NonCovid', inplace=True)
l_map = l_class[["unique_identifier", patient_status, "covid_state"]]

l_data = l_class.drop([patient_status, "covid_state"], axis=1)
l_data.set_index("unique_identifier", inplace=True)
l_data.columns = [y for x, y in l_data.columns]
l_data
# -

# ## Metabolomics
#
# ### Experiment notes
#
# "GC-MS raw files were processed using a software suite developed in-house that is available at https://github.com/coongroup. Following data acquisition, raw EI-GC/MS spectral data was deconvolved into chromatographic features and then grouped into features based on co-elution. Only features with at least 10 fragment ions and present in 33% of samples were kept. Feature groups from samples and background were compared, and only feature groups greater than 3-fold higher than background were retained. Compound identifications for the metabolites analyzed were assigned by comparing deconvolved high-resolution spectra against unit-resolution reference spectra present in the NIST 12 MS/EI library as well as to authentic standards run in-house. To calculate spectral similarity between experimental and reference spectra, a weighted dot product calculation was used. Metabolites lacking a confident identification were classified as “Unknown metabolites” and appended a unique identifier based on retention time. Peak heights of specified quant m/z were used to represent feature (metabolite) abundance. The data set was also processed through, where we applied a robust linear regression approach, rlm() function (Marazzi et al., 1993), non- log2 transformed intensity values versus run order, to normalize for run order effects on signal. AEX-LC-MS/MS: raw files were processed using Xcalibur Qual Browser (v4.0.27.10, Thermo Scientific) with results exported and further processed using Microsoft Excel 2010. The prepared standard solution was used to locate appropriate peaks for peak area analysis." - [Overmyer et al 2020](https://dx.doi.org/10.1101%2F2020.07.17.20156513)
#
# ### Data notes
#
# Data is a matrix of continuous values. Column names correspond to metabolite names and row names correspond to unique sample identifiers. Metabolites with no known annotation are also recorded and assigned a unique identifier based on their mass and charge. There are two subsets of metabolites - discovery and targeted. These were originally split into two sets for preliminary investigation but are later recombined before entering the pipeline.

# +
m = data_omics[1][["unique_identifier", "standardized_name", "normalized_abundance", "keep", "omics_id"]]
m = m[m["keep"] == "1"]
m.drop("keep", axis=1, inplace=True)

m_d = m[m["omics_id"] == 3]
m_d.drop("omics_id", axis=1, inplace=True)

m_t = m[m["omics_id"] == 4]
m_t.drop("omics_id", axis=1, inplace=True)

m_d = m_d.pivot(index="unique_identifier", columns="standardized_name")
m_d.columns = [x[1] for x in m_d.columns]

# this maps the sample to patient
m_id = data_omics[1][["unique_identifier", patient_status]]
m_id.drop_duplicates(inplace=True)
m_id.set_index("unique_identifier", inplace=True)

m_d = m_d.merge(m_id, left_index=True, right_index=True, how="left")
m_d.reset_index(inplace=True)
m_d.set_index(["unique_identifier", patient_status], inplace=True)
m_d.sort_index(level=patient_status, inplace=True)

# remove null and assign sample groups
m_d_class = m_d.reset_index()
m_d_class.dropna(inplace=True)

m_d_class["covid_state"] = m_d_class[patient_status]
m_d_class["covid_state"].replace(regex=r'^C.*', value='Covid', inplace=True)
m_d_class["covid_state"].replace(regex=r'^NC.*', value='NonCovid', inplace=True)

m_d_data = m_d_class.drop([patient_status, "covid_state"], axis=1)
m_d_data.set_index("unique_identifier", inplace=True)
keep = ~m_d_data.eq(m_d_data.iloc[:, 0], axis=0).all(1)
m_d_data = m_d_data[keep]

m_d_map = m_d_class[["unique_identifier", patient_status, "covid_state"]]
m_d_map.set_index("unique_identifier", inplace=True)
m_d_map = m_d_map[keep].reset_index()
m_d_data

# +
m_t = m_t.pivot(index="unique_identifier", columns="standardized_name")
m_t.columns = [x[1] for x in m_t.columns]

m_t = m_t.merge(m_id, left_index=True, right_index=True, how="left")
m_t.reset_index(inplace=True)
m_t.set_index(["unique_identifier", patient_status], inplace=True)
m_t.sort_index(level=patient_status, inplace=True)

# remove null and assign sample groups
m_t_class = m_t.reset_index()
m_t_class.dropna(inplace=True)

m_t_class["covid_state"] = m_t_class[patient_status]
m_t_class["covid_state"].replace(regex=r'^C.*', value='Covid', inplace=True)
m_t_class["covid_state"].replace(regex=r'^NC.*', value='NonCovid', inplace=True)
m_t_map = m_t_class[["unique_identifier", patient_status, "covid_state"]]

m_t_data = m_t_class.drop([patient_status, "covid_state"], axis=1)
m_t_data.set_index("unique_identifier", inplace=True)
m_t_data
# -

# ## Proteomics
#
# ### Experiment notes
#
# "Shotgun proteomics raw files were searched using MaxQuant quantitative software package (Cox et al., 2014) (version 1.6.10.43) against UniProt Homo Sapiens database (downloaded on 6.18.2019), containing protein isoforms and computationally predicted proteins. If not specified, default MaxQuant settings were used. LFQ quantification was performed using LFQ minimum ratio count of 1 and no MS/MS requirement for LFQ comparisons. iBAQ quantitation and “match between runs” were enabled with default settings. ITMS MS/MS tolerance was set to 0.35 Da. Lists of quantified protein groups were filtered to remove reverse identifications, potential contaminants, and proteins identified only by a modification site. LFQ abundance values were log2 transformed. Missing quantitative values were imputed by randomly drawing values from the left tail of the normal distribution of all measured protein abundance values (Tyanova et al., 2016). Protein groups that contained more than 50% missing values were removed from final analyses. Relative standard deviations (RSDs) for each protein group quantified across all seven technical replicates of healthy plasma controls were calculated, and proteins with RSD greater than 30% were removed from final analyses. PRM: identification and quantification of targeted peptides for PRM analysis were performed using Skyline open access software package (version 20.1). 4-5 most intense and specific transitions were used to quantify peptide abundances, and area-under-the-curve measurements for each peptide were exported for further analysis." - [Overmyer et al 2020](https://dx.doi.org/10.1101%2F2020.07.17.20156513)
#
# ### Data notes
#
# Data is a matrix of continuous values. Column names correspond to protein names and row names correspond to unique sample identifiers.

# +
p = data_omics[2][["unique_identifier", "standardized_name", "normalized_abundance", "keep"]]
p = p[p["keep"] == "1"]
p.drop("keep", axis=1, inplace=True)
p = p.pivot(index="unique_identifier", columns="standardized_name")
p.columns = [x[1] for x in p.columns]

# this maps the sample to patient
p_id = data_omics[2][["unique_identifier", patient_status]]
p_id.drop_duplicates(inplace=True)
p_id.set_index("unique_identifier", inplace=True)

p = p.merge(p_id, left_index=True, right_index=True, how="left")
p.reset_index(inplace=True)
p.set_index(["unique_identifier", patient_status], inplace=True)
p.sort_index(level=patient_status, inplace=True)

# remove null and assign sample groups
p_class = p.reset_index()
p_class.dropna(inplace=True)

p_class["covid_state"] = p_class[patient_status]
p_class["covid_state"].replace(regex=r'^C.*', value='Covid', inplace=True)
p_class["covid_state"].replace(regex=r'^NC.*', value='NonCovid', inplace=True)
p_map = p_class[["unique_identifier", patient_status, "covid_state"]]

p_data = p_class.drop([patient_status, "covid_state"], axis=1)
p_data.set_index("unique_identifier", inplace=True)
p_data
# -

# ## RNA-Seq Data Processing
#
# "All RNA transcripts were downloaded from the NCBI refseq ftp site (wget ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/mRNA_Prot/∗.rna.fna.gz ). Only mRNA (accessions NM_xxxx and XM_xxxx) and rRNA (excluding 5.8S) was then extracted, and immunoglobulin transcripts were downloaded from ENSEMBL (IG_C, IG_D, IG_J and IG_V ). We created a file mapping accession numbers to gene symbols, and then used rsem-prepare-reference to build a bowtie-2 reference database. Fastq files were trimmed and filtered using a custom algorithm tailored to improve quality scores and maximize retained reads in paired-end data. RNA-Seq expression estimation was performed by RSEM v 1.3.0 (parameters: seed-length=20, no-qualities, bowtie2-k=200, bowtie2-sensitivity-level=sensitive) (Li and Dewey, 2011), with bowtie-2 (v 2.3.4.1) for the alignment step (Langmead and Salzberg, 2012), using the custom hg38 reference described above. After the collation of expression estimates, hemoglobin transcripts were removed from further analysis, and TPM values were rescaled to total 1,000,000 in each sample. Differential Expression analysis was performed using the EBSeq package (v 1.26.0) (Leng et al., 2013) in R (v 3.6.2)." - [Overmyer et al 2020](https://dx.doi.org/10.1101%2F2020.07.17.20156513)
#
# ### Data notes
#
# Data is a matrix of continuous values. Zero values exist. Column names correspond to unique gene names and row names correspond to sample identifiers.

# +
t = data_omics[3][["unique_identifier", "standardized_name", "normalized_abundance", "keep"]]
t = t[t["keep"] == "1"]
t.drop("keep", axis=1, inplace=True)
t = t.pivot(index="unique_identifier", columns="standardized_name")
t.columns = [x[1] for x in t.columns]

# this maps the sample to patient
t_id = data_omics[3][["unique_identifier", patient_status]]
t_id.drop_duplicates(inplace=True)
t_id.set_index("unique_identifier", inplace=True)

t = t.merge(t_id, left_index=True, right_index=True, how="left")
t.reset_index(inplace=True)
t.set_index(["unique_identifier", patient_status], inplace=True)

t.sort_index(level=patient_status, inplace=True)

# remove null and assign sample groups
t_class = t.reset_index()
t_class.dropna(inplace=True)

t_class["covid_state"] = t_class[patient_status]
t_class["covid_state"].replace(regex=r'^C.*', value='Covid', inplace=True)
t_class["covid_state"].replace(regex=r'^NC.*', value='NonCovid', inplace=True)
t_map = t_class[["unique_identifier", patient_status, "covid_state"]]

t_data = t_class.drop([patient_status, "covid_state"], axis=1)
t_data.set_index("unique_identifier", inplace=True)
t_data
# -

# ## Sample filtering
#
# Our pipeline and underlying methods require all samples to be identical across the blocks of omics data. These steps filter out all samples which are not represented in each omics data block.

# our framework only ingests data with matched samples
data_clean = [l_data, m_d_data, m_t_data, p_data, t_data]
maps = [l_map, m_d_map, m_t_map, p_map, t_map]
maps = [x.set_index("Albany_sampleID") for x in maps]
common = [set(x.index.tolist()) for x in maps]
common = set.intersection(*common)
maps = [x.loc[list(common)] for x in maps]
matched = [i.loc[j["unique_identifier"]] for i, j in tuple(zip(data_clean, maps))]
[x.shape for x in matched]


# ## Sample mapping (Covid state)
#
# Unique sample identifiers are matched to their biologically relevant sample categories. Long feature names are also shortened, and mapping files containing the original names are generated for later recovery if necessary.

# +
# remap identifiers to sample types
def map_common_id(metadata, data):
    metadata = pd.DataFrame(metadata.set_index("unique_identifier")["Albany_sampleID"])
    data = pd.merge(metadata, data, left_index=True, right_index=True)
    data = data.reset_index().set_index("Albany_sampleID")
    data.index.name = None
    return data.drop("unique_identifier", axis=1)

# long feature names break the pipeline
def shorten_colnames(data):
    name_long = pd.Series(data.columns, name="long")
    data.columns = [x[:24] for x in data.columns]
    name_short = pd.Series(data.columns, name="short")
    name_map = pd.concat([name_short, name_long], axis=1)
    return data, name_map

classes_clean = [l_class, m_d_class, m_t_class, p_class, t_class]
data_final = [map_common_id(x, y) for x, y in tuple(zip(classes_clean, matched))]
data_short = [shorten_colnames(x) for x in data_final]
data_map = [x[1] for x in data_short]
data_final = [x[0] for x in data_short]
# -

# ## Data for pipeline input
#
# Output the preprocessed data into the relevant files.

# +
outcovid = [
    "../../data/MSV000085703/covid/data_lipidomics.tsv" ,
    "../../data/MSV000085703/covid/data_metabolomicsdiscovery.tsv" ,
    "../../data/MSV000085703/covid/data_metabolomicstargeted.tsv" ,
    "../../data/MSV000085703/covid/data_proteomics.tsv" ,
    "../../data/MSV000085703/covid/data_transcriptomics.tsv" ,
]

outhfd45 = [
    "../../data/MSV000085703/hfd45/data_lipidomics.tsv" ,
    "../../data/MSV000085703/hfd45/data_metabolomicsdiscovery.tsv" ,
    "../../data/MSV000085703/hfd45/data_metabolomicstargeted.tsv" ,
    "../../data/MSV000085703/hfd45/data_proteomics.tsv" ,
    "../../data/MSV000085703/hfd45/data_transcriptomics.tsv" ,
]

# !mkdir -p ../../data/MSV000085703/covid ../../data/MSV000085703/hfd45

[x.to_csv(y, sep="\t") for x, y in tuple(zip(data_final, outcovid))]
[x[:100].to_csv(y, sep="\t") for x, y in tuple(zip(data_final, outhfd45))]

# +
outmeta = [
    "../../data/MSV000085703/lipidomics_class.tsv" ,
    "../../data/MSV000085703/metabolomicsdiscovery_class.tsv" ,
    "../../data/MSV000085703/metabolomicstargeted_class.tsv" ,
    "../../data/MSV000085703/proteomics_class.tsv" ,
    "../../data/MSV000085703/transcriptomics_class.tsv" ,
]

def map_sample_type(metadata, data):
    data = pd.merge(metadata, data.set_index("Albany_sampleID"), left_index=True, right_index=True)
    data = data[~data.index.duplicated(keep='first')]
    return pd.DataFrame(data["covid_state"])

classes_final = [map_sample_type(x, y) for x, y in tuple(zip(data_final, classes_clean))]
[x.to_csv(y, sep="\t") for x, y in tuple(zip(classes_final, outmeta))]

# all these files are identical, choose one to copy to avoid confusion
classes_final[0].to_csv("../../data/MSV000085703/covid/classes_diablo.tsv", sep="\t")
classes_final[0][:100].to_csv("../../data/MSV000085703/hfd45/classes_diablo.tsv", sep="\t")

# +
outmap = [
    "../../data/MSV000085703/lipidomics_featuremap.tsv" ,
    "../../data/MSV000085703/metabolomicsdiscovery_featuremap.tsv" ,
    "../../data/MSV000085703/metabolomicstargeted_featuremap.tsv" ,
    "../../data/MSV000085703/proteomics_featuremap.tsv" ,
    "../../data/MSV000085703/transcriptomics_featuremap.tsv" ,
]

[x.to_csv(y, sep="\t") for x, y in tuple(zip(data_map, outmap))]
# -

# ## Combine metabolomics data
#
# The authors combined both sets of `discovery` and `targeted` metabolomics data. These were originally separated to investigate the possibility that the individual blocks of omics data would be informative during integration.

# +
m_covid = "../../data/MSV000085703/covid/data_metabolomics.tsv"
m_hfd45 = "../../data/MSV000085703/hfd45/data_metabolomics.tsv"

metabolomics = data_final[1].merge(data_final[2], left_index=True, right_index=True)
metabolomics.to_csv(m_covid, sep="\t")
metabolomics[:100].to_csv(m_hfd45, sep="\t")
# -

# ## Sample mapping (Covid severity)
#
# The `HFD-45` metric created by the authors measures the days a patient is in hospital. The authors use this to classify patients by disease severity.

# +
# transcriptomics covid vs noncovid
dge_path_hfd45 = "../../data/MSV000085703/hfd45/data_transcriptomics_dge_hfd45.tsv"

def extract_hfd45(outpath, data):
    data = data[:100]
    data.to_csv(outpath, sep="\t")
    return data

hfd45_paths = [
    "../../data/MSV000085703/hfd45/data_lipidomics.tsv",
    "../../data/MSV000085703/hfd45/data_metabolomicsdiscovery.tsv",
    "../../data/MSV000085703/hfd45/data_metabolomicstargeted.tsv",
    "../../data/MSV000085703/hfd45/data_proteomics.tsv",
    "../../data/MSV000085703/hfd45/data_transcriptomics.tsv",
]

# transcriptomics hfd45 spectrum
patient_info = patient_metadata.set_index("Albany_sampleID")
patient_info = classes_final[4].merge(patient_info, left_index=True, right_index=True)
patient_info = patient_info[["covid_state", "Hospital_free_days_45"]]
patient_info = patient_info[patient_info["covid_state"] == "Covid"].drop("covid_state", axis=1)

# assign severity score based on median value (26) reported by authors
patient_info[patient_info["Hospital_free_days_45"] <= 26] = 1
patient_info[patient_info["Hospital_free_days_45"] > 26] = 0
patient_info[patient_info["Hospital_free_days_45"] == 1] = "More severe"
patient_info[patient_info["Hospital_free_days_45"] == 0] = "Less severe"
patient_info.to_csv("../../data/MSV000085703/hfd45/classes_diablo.tsv", sep="\t")

# subset data accordingly (take only the 100 covid patients)
covid_only = patient_info.index
hfd45_data = [extract_hfd45(x, y) for x, y in list(zip(hfd45_paths, data_final))]


# -

# ## Data imputation
#
# Transcriptomics data was sparse, and was imputed as part of the pipeline. For reproducibility, we included the imputed data file used, which can be provided as an input file directly.
#
# Note that the `--icomp` flag can be enabled in the main pipeline script, which will perform the imputation step internally given `N` components.

# # Appendix
#
# This appendix is included for completeness only and describes a few ways to work around common problems in data.
#
# ## Handling class imbalance in data
#
# There are class imbalances between `COVID` and `NON-COVID` samples. It is possible to subsampled a set of samples from the `COVID` sample type. In this example, we first investigated if the distribution of the samples are equal. The subset contains `25` samples while the full set contains `100` samples.

# +
# check that the distribution of data in all covid cases is similar to the subset
def check_dist(data, meta):
    f = plt.figure(figsize=(10,3))
    ax1 = f.add_subplot(121)
    ax1.set_title(meta + " subset")
    ax2 = f.add_subplot(122)
    ax2.set_title(meta + " full set")
    pd.DataFrame(data[75:100].values.flatten()).boxplot(ax=ax1)
    pd.DataFrame(data[:100].values.flatten()).boxplot(ax=ax2)
    plt.show()

meta = [
    "lipidomics",
    "metabolomics discovery",
    "metabolomics targeted",
    "proteomics",
    "transcriptomics",
]
   
[check_dist(x, y) for x, y in list(zip(data_final, meta))]
# -

# ## Feature selection
#
# To improve run time and performance, feature selection can be performed. Ideally, features can be iteratively removed and the classification accuracy observed. However, in biological datasets where features exceed samples (for example a sample may have 10000 associated genes), this is not always possible.

# Features were selected by performing a t test on the data after checking assumptions of independence, distributions and variance. Because of the absence of raw patient data, the t test was performed on normalised abundance measures and did not go through a conventional differential expression pipeline with `edgeR` and `limma`. We note that the purpose of this step is to lower the feature count and not to use these results in downstream analysis.

# +
dge_path_covid = "../../data/MSV000085703/data_transcriptomics_dge_covid.tsv"
dge_path_covid_balanced = "../../data/MSV000085703/balanced/data_transcriptomics_dge_covid.tsv"

# we want to reduce features in transcriptomics data only
df = classes_final[4].merge(data_final[4], left_index=True, right_index=True)
df = df.set_index("covid_state")

# perform a welch t test for covid vs noncovid
ttest = stats.ttest_ind(df.loc["Covid"], df.loc["NonCovid"], equal_var=False)
mean_covid = df.loc["Covid"].mean()
mean_noncovid = df.loc["NonCovid"].mean()

# get bh adjusted p values
ttest = pd.DataFrame(multipletests(ttest[1], method="fdr_bh")[1])
ttest.index = df.columns
dgelist = ttest[ttest < 0.05].dropna().index

# independent (no patients are represented in more than 1 sample group)

# variances (similar but not identical)
stats.levene(mean_covid, mean_noncovid, center="mean")

# normal distribution 
df.loc["Covid"].mean().hist(alpha=0.5)
df.loc["NonCovid"].mean().hist(alpha=0.5, color="orange")
plt.show()

sampling_difference = mean_covid - mean_noncovid
stats.shapiro(sampling_difference)
stats.probplot(sampling_difference, plot=plt, rvalue= True)
plt.show()
# -



StarcoderdataJupyterScriptsDedupFiltered