title: “Introduction to the XINA pagkage”
author: “Lang Ho Lee, Sasha A. Singh”
date: “February 6, 2019”
vignette: >
%
%
%
output:
knitr:::html_vignette:
df_print: kable
toc: true
number_sections: true

1. Introduction

Quantitative proteomics experiments, using for instance isobaric tandem mass tagging approaches, are conducive to measuring changes in protein abundance over multiple time points in response to one or more conditions or stimulations. The aim of XINA is to determine which proteins exhibit similar patterns within and across experimental conditions, since proteins with co-abundance patterns may have common molecular functions. XINA imports multiple datasets, tags dataset in silico, and combines the data for subsequent subgrouping into multiple clusters. The result is a single output depicting the variation across all conditions. XINA, not only extracts co-abundance profiles within and across experiments, but also incorporates protein-protein interaction databases and integrative resources such as KEGG to infer interactors and molecular functions, respectively, and produces intuitive graphical outputs.

1-1. Main contribution

An easy-to-use software for non-expert users of clustering and network analyses.

1-2. Data inputs

Any type of quantitative proteomics data, labeled or label-free

2. XINA websites

https://cics.bwh.harvard.edu/software http://bioconductor.org/packages/XINA/ https://github.com/langholee/XINA/

3. XINA installation

XINA requires R>=3.5.0.

# Install from Bioconductor
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("XINA")

# Install from Github
install.packages("devtools")
library(devtools)
install_github("langholee/XINA")

The first step is to call XINA

library(XINA)

To follow this vignette, you may need the following packages

install.packages("igraph")
install.packages("ggplot2")
BiocManager::install("STRINGdb")

4. Example theoretical dataset

We generated an example dataset to show how XINA can be used for your research. To demonstrate XINA functions and allow users to perform similar exercises, we included a module that can generate multiplexed time-series datasets using theoretical data. This data consists of three treatment conditions, ‘Control’, ‘Stimulus1’ and ‘Stimulus2’. Each condition has time series data from 0 hour to 72 hours. As an example, we chose the mTOR pathway to be differentially regulated across the three conditions.

# Generate random multiplexed time-series data
random_data_info <- make_random_xina_data()

# The number of proteins
random_data_info$size
## [1] 500
# Time points
random_data_info$time_points
## [1] "0hr"  "2hr"  "6hr"  "12hr" "24hr" "48hr" "72hr"
# Three conditions
random_data_info$conditions
## [1] "Control"   "Stimulus1" "Stimulus2"

Read and check the randomly generated data

Control <- read.csv("Control.csv", check.names=FALSE, stringsAsFactors = FALSE)
Stimulus1 <- read.csv("Stimulus1.csv", check.names=FALSE, stringsAsFactors = FALSE)
Stimulus2 <- read.csv("Stimulus2.csv", check.names=FALSE, stringsAsFactors = FALSE)

head(Control)
##   Accession                                                                  Description 0hr    2hr    6hr   12hr
## 1    FAM49A                                  family with sequence similarity 49 member A 0.5 0.4577 0.0243 0.3614
## 2     PXDC1                                                       PX domain containing 1 0.5 0.8565 0.0945 0.1238
## 3    VCPIP1                             valosin containing protein interacting protein 1 0.5 0.3380 0.7030 0.7362
## 4     KRT13                                                                   keratin 13 0.5 0.1865 0.6355 0.0558
## 5      PJA1                                         praja ring finger ubiquitin ligase 1 0.5 0.1222 0.0991 0.8449
## 6    MCIDAS multiciliate differentiation and DNA synthesis associated cell cycle protein 0.5 0.4273 0.0478 0.7328
##     24hr   48hr   72hr
## 1 0.2160 0.6698 0.4698
## 2 0.3542 0.4111 0.7061
## 3 0.8331 0.2453 0.7024
## 4 0.1011 0.7718 0.6007
## 5 0.7894 0.0749 0.2802
## 6 0.3324 0.3450 0.5003
head(Stimulus1)
##   Accession                                                                  Description 0hr    2hr    6hr   12hr
## 1    FAM49A                                  family with sequence similarity 49 member A 0.5 0.3034 0.0999 0.7214
## 2     PXDC1                                                       PX domain containing 1 0.5 0.0170 0.7694 0.7538
## 3    VCPIP1                             valosin containing protein interacting protein 1 0.5 0.5652 0.1085 0.3918
## 4     KRT13                                                                   keratin 13 0.5 0.4513 0.2071 0.5934
## 5      PJA1                                         praja ring finger ubiquitin ligase 1 0.5 0.2832 0.1420 0.3540
## 6    MCIDAS multiciliate differentiation and DNA synthesis associated cell cycle protein 0.5 0.2163 0.1614 0.9797
##     24hr   48hr   72hr
## 1 0.7297 0.7611 0.8759
## 2 0.0164 0.9957 0.1262
## 3 0.0741 0.8811 0.2479
## 4 0.6506 0.7603 0.4085
## 5 0.4437 0.2379 0.6623
## 6 0.4440 0.6513 0.8546
head(Stimulus2)
##   Accession                                                                  Description 0hr    2hr    6hr   12hr
## 1    FAM49A                                  family with sequence similarity 49 member A 0.5 0.0826 0.0478 0.8025
## 2     PXDC1                                                       PX domain containing 1 0.5 0.6798 0.7823 0.4658
## 3    VCPIP1                             valosin containing protein interacting protein 1 0.5 0.3959 0.9274 0.8680
## 4     KRT13                                                                   keratin 13 0.5 0.9661 0.8917 0.1466
## 5      PJA1                                         praja ring finger ubiquitin ligase 1 0.5 0.2513 0.1593 0.4043
## 6    MCIDAS multiciliate differentiation and DNA synthesis associated cell cycle protein 0.5 0.2287 0.7015 0.8684
##     24hr   48hr   72hr
## 1 0.7914 0.8905 0.8539
## 2 0.2469 0.7208 0.4229
## 3 0.1506 0.0450 0.4395
## 4 0.1027 0.2784 0.8220
## 5 0.5423 0.3286 0.4144
## 6 0.6455 0.5873 0.9634

Since XINA needs to know which columns have the kinetics data matrix, the user should give a vector containing column names of the kinetics data matrix. These column names have to be the same in all input datasets (Control input, Stimulus1 input and Stimulus2 input).

head(Control[random_data_info$time_points])
##   0hr    2hr    6hr   12hr   24hr   48hr   72hr
## 1 0.5 0.4577 0.0243 0.3614 0.2160 0.6698 0.4698
## 2 0.5 0.8565 0.0945 0.1238 0.3542 0.4111 0.7061
## 3 0.5 0.3380 0.7030 0.7362 0.8331 0.2453 0.7024
## 4 0.5 0.1865 0.6355 0.0558 0.1011 0.7718 0.6007
## 5 0.5 0.1222 0.0991 0.8449 0.7894 0.0749 0.2802
## 6 0.5 0.4273 0.0478 0.7328 0.3324 0.3450 0.5003

5. Package features

XINA is an R package and can examine, but is not limited to, time-series omics data from multiple experiment conditions. It has three modules: 1. Model-based clustering analysis, 2. coregulation analysis, and 3. Protein-protein interaction network analysis (we used STRING DB for this practice).

5.1 Clustering analysis using model-based clustering or k-means clustering algorithm

XINA implements model-based clustering to classify features (genes or proteins) depending on their expression profiles. The model-based clustering optimizes the number of clusters at minimum Bayesian information criteria (BIC). Model-based clustering is fulfilled by the ‘mclust’ R package [https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5096736/], which was used by our previously developed tool mIMT-visHTS [https://www.ncbi.nlm.nih.gov/pubmed/26232111]. By default, XINA performs sum-normalization for each gene/protein time-series profile [https://www.ncbi.nlm.nih.gov/pubmed/19861354]. This step is done to standardize all datasets. Most importantly, XINA assigns an electronic tag to each dataset’s proteins (similar to TMT) in order to combine the multiple datasets (Super dataset) for subsequent clustering.

XINA uses the ‘mclust’ package for the model-based clustering. ‘mclust’ requires the fixed random seed to get reproducible clustering results.

set.seed(0)

‘nClusters’ is the number of desired clusters. ‘mclust’ will choose the most optimized number of clusters by considering the Bayesian information criteria (BIC). BIC of ‘mclust’ is the negative of normal BIC, thus the higher BIC, the more optimized clustering scheme in ‘mclust’, while lowest BIC is better in statistics.

# Data files
data_files <- paste(random_data_info$conditions, ".csv", sep='')
data_files
## [1] "Control.csv"   "Stimulus1.csv" "Stimulus2.csv"
# time points of the data matrix
data_column <- random_data_info$time_points
data_column
## [1] "0hr"  "2hr"  "6hr"  "12hr" "24hr" "48hr" "72hr"

Run the model-based clustering

# Run the model-based clusteirng
clustering_result <- xina_clustering(data_files, data_column=data_column, nClusters=20)

XINA also supports k-means clustering as well as the model-based clustering

clustering_result_km <- xina_clustering(data_files, data_column=data_column, nClusters=20, chosen_model='kmeans')

For visualizing clustering results, XINA draws line graphs of the clustering results using ‘plot_clusters’.

library(ggplot2)
theme1 <- theme(title=element_text(size=8, face='bold'),
                axis.text.x = element_text(size=7),
                axis.text.y = element_blank(),
                axis.ticks.x = element_blank(),
                axis.ticks.y = element_blank(),
                axis.title.x = element_blank(),
                axis.title.y = element_blank())
plot_clusters(clustering_result, ggplot_theme=theme1)

XINA calculates sample condition composition, for example the sample composition in the cluster 28 is higher than 95% for Stimulus2. ‘plot_condition_composition’ plots these compositions as pie-charts. Sample composition information is insightful because we can find which specific patterns are closely related with each stimulus.

theme2 <- theme(legend.key.size = unit(0.3, "cm"),
                legend.text=element_text(size=5),
                title=element_text(size=7, face='bold'))
condition_composition <- plot_condition_compositions(clustering_result, ggplot_theme=theme2)

tail(condition_composition)
##    Cluster Condition  N Percent_ratio
## 41      14 Stimulus2  2          3.45
## 42      15   Control  1          1.67
## 43      15 Stimulus1 57         95.00
## 44      15 Stimulus2  2          3.33
## 45      16   Control  1          1.72
## 46      16 Stimulus2 57         98.28

5.2 coregulation analysis

XINA supposes that proteins that comigrate between clusters in response to a given condition are more likely to be coregulated at the biological level than other proteins within the same clusters. For this module, at least two datasets to be compared are needed. XINA supposes features assigned to the same cluster in an experiment condition as a coregulated group. XINA traces the comigrated proteins in different experiment conditions and finds signficant trends by 1) the number of member features (proteins) and 2) the enrichment test using the Fishers exact test. The comigrations are displayed via an alluvial plot. In XINA the comigration is defined as a condition of proteins that show the same expression pattern, classified and evaluated by XINA clustering, in at least two dataset conditions. If there are proteins that are assigned to the same cluster in more than two datasets, XINA considers those proteins to be comigrated. XINA’s ‘alluvial_enriched’ is designed to find these comigrations and draws alluvial plots for visualizing the found comigrations.

classes <- as.vector(clustering_result$condition)
classes
## [1] "Control"   "Stimulus1" "Stimulus2"
all_cor <- alluvial_enriched(clustering_result, classes)
## [1] "length(selected_conditions) > 2, so XINA can't apply the enrichment filter\n            Can't apply the enrichment filter, so pval_threshold is ignored"

head(all_cor)
##   Control Stimulus1 Stimulus2 Comigration_size RowNum PValue Pvalue.adjusted TP FP FN TN Alluvia_color
## 1       0         0         1                9      1     NA              NA NA NA NA NA       #BEBEBE
## 2       0         0         2                4      2     NA              NA NA NA NA NA       #BEBEBE
## 3       0         0         3                4      3     NA              NA NA NA NA NA       #BEBEBE
## 4       0         0         4                4      4     NA              NA NA NA NA NA       #BEBEBE
## 5       0         0         5                3      5     NA              NA NA NA NA NA       #BEBEBE
## 6       0         0         6                9      6     NA              NA NA NA NA NA       #BEBEBE

You can narrow down comigrations by using the size (the number of comigrated proteins) filter.

cor_bigger_than_5 <- alluvial_enriched(clustering_result, classes, comigration_size=5)
## [1] "length(selected_conditions) > 2, so XINA can't apply the enrichment filter\n            Can't apply the enrichment filter, so pval_threshold is ignored"

head(cor_bigger_than_5)
##    Control Stimulus1 Stimulus2 Comigration_size RowNum PValue Pvalue.adjusted TP FP FN TN Alluvia_color
## 1        0         0         1                9      1     NA              NA NA NA NA NA       #BEBEBE
## 6        0         0         6                9      6     NA              NA NA NA NA NA       #BEBEBE
## 7        0         0         7               12      7     NA              NA NA NA NA NA       #BEBEBE
## 8        0         0         8               12      8     NA              NA NA NA NA NA       #BEBEBE
## 11       0         0        11                7     11     NA              NA NA NA NA NA       #BEBEBE
## 12       0         0        12               13     12     NA              NA NA NA NA NA       #BEBEBE

5.3 Network analysis

XINA conducts protein-protein interaction (PPI) network analysis through implementing ‘igraph’ and ‘STRINGdb’ R packages. XINA constructs PPI networks for comigrated protein groups as well as individual clusters of a specific experiment (dataset) condition. In the constructed networks, XINA finds influential players by calculating various network centrality calculations including betweenness, closeness and eigenvector scores. For the selected comigrated groups, XINA can calculate an enrichment test based on gene ontology and KEGG pathway terms to help understanding comigrated groups.

XINA’s example dataset is from human gene names, so download human PPI database from STRING DB and run XINA PPI network analysis.

library(STRINGdb)
string_db <- STRINGdb$new( version="10", species=9606, score_threshold=0, input_directory="" )
string_db

xina_result <- xina_analysis(clustering_result, string_db)

You can draw PPI networks of all the XINA clusters using ‘xina_plots’ function easily. PPI network plots will be stored in the working directory

# XINA network plots labeled gene names
xina_plot_all(xina_result, clustering_result)

If you want to see more, please check “README.md” of our Github XINA repository.