pkgdown/extra.css

Skip to contents

This vignette will guide you through the initial steps of a proteomics pipeline, which is the exploratory data analysis. Except from HDAnalyzeR we will import viridis for the color palettes.

📓 This vignette is just a basic example of how to explore the data and perform dimensionality reduction to identify potential issues and trends. However, every dataset is different and may require different exploration techniques other than the ones contained in the package.

⚠️ Please make sure to preprocess the data and metadata before proceeding with the analysis. For the most common proteomics data formats, we provide some recommendations in the package documentation under articles “Data Preparation for HDAnalyzeR: What You Need Before Using the Package”.

Loading the Data

Let’s start with loading the example data and metadata that come with the package and initialize the HDAnalyzeR object.

hd_obj <- hd_initialize(dat = example_data, 
                        metadata = example_metadata, 
                        is_wide = FALSE, 
                        sample_id = "DAid",
                        var_name = "Assay",
                        value_name = "NPX")

Exploratory Data Analysis

Now, we will perform an automated exploratory data analysis (QC check) on the data. hd_qc_summary() will return the percentage of missing values for each column and row and histograms of their distributions both for data and metadata, the variable-variable (in this case protein-protein) correlation and the metadata variable distributions.

⚠️ This should not be confused with the quality control of the data, which is a different process that is usually performed right after the proteomics experiments. The exploratory data analysis is a first step to understand the data and identify potential issues.

As variable we should pass the column name of the metadata that contains the different classes, for example the case and control groups, the different diseases, or the different stages of a disease. The palettes are optional and can be used to color the plots of the metadata variable distributions. You should pass a list with the variable name and the palette name, for example list(Sex = c(“F” = “red”, “M” = “blue”), Disease = “cancers12”). As you can see, the palettes can be either a named vector or a character vector with the name of a palette from the package. You can display all available HDAnalyzeR palettes by running hd_show_palettes().

qc_data <- hd_qc_summary(hd_obj, 
                         variable = "Disease", 
                         palette = list(Disease = "cancers12", Sex = "sex"), 
                         cor_threshold = 0.7)
#> [1] "Summary:"
#> [1] "Note: In case of long output, only the first 10 rows are shown. To see the rest display the object with view()"
#> [1] "Number of samples: 586"
#> [1] "Number of variables: 101"
#> [1] "--------------------------------------"
#> [1] "categorical : 1"
#> [1] "continuous : 100"
#> [1] "--------------------------------------"
#> [1] "NA percentage in each column:"
#> # A tibble: 91 × 2
#>    column   na_percentage
#>    <chr>            <dbl>
#>  1 ACE2               6.1
#>  2 ACTA2              6.1
#>  3 ACTN4              6.1
#>  4 ADAM15             6.1
#>  5 ADAMTS16           6.1
#>  6 ADH4               6.1
#>  7 AKR1C4             6.1
#>  8 AMBN               6.1
#>  9 AMN                6.1
#> 10 AOC1               6.1
#> # ℹ 81 more rows
#> [1] "--------------------------------------"
#> [1] "NA percentage in each row:"
#> # A tibble: 144 × 2
#>    DAid    na_percentage
#>    <chr>           <dbl>
#>  1 DA00450          57.4
#>  2 DA00482          53.5
#>  3 DA00542          53.5
#>  4 DA00003          50.5
#>  5 DA00463          46.5
#>  6 DA00116          43.6
#>  7 DA00475          42.6
#>  8 DA00578          42.6
#>  9 DA00443          41.6
#> 10 DA00476          35.6
#> # ℹ 134 more rows
#> [1] "--------------------------------------"
#> [1] "Protein-protein correlations above 0.7:"
#>   Protein1 Protein2 Correlation
#> 1  ATP5IF1    AIFM1        0.76
#> 2    AXIN1 ARHGEF12        0.76
#> 3    AIFM1  ATP5IF1        0.76
#> 4 ARHGEF12    AXIN1        0.76
#> 5 ARHGEF12    AIFM1        0.71
#> 6    AIFM1 ARHGEF12        0.71
#> [1] "--------------------------------------"
#> [1] "Summary:"
#> [1] "Note: In case of long output, only the first 10 rows are shown. To see the rest display the object with view()"
#> [1] "Number of samples: 586"
#> [1] "Number of variables: 9"
#> [1] "--------------------------------------"
#> [1] "categorical : 7"
#> [1] "continuous : 2"
#> [1] "--------------------------------------"
#> [1] "NA percentage in each column:"
#> # A tibble: 1 × 2
#>   column na_percentage
#>   <chr>          <dbl>
#> 1 Grade           91.5
#> [1] "--------------------------------------"
#> [1] "NA percentage in each row:"
#> # A tibble: 536 × 2
#>    DAid    na_percentage
#>    <chr>           <dbl>
#>  1 DA00001          11.1
#>  2 DA00002          11.1
#>  3 DA00003          11.1
#>  4 DA00004          11.1
#>  5 DA00005          11.1
#>  6 DA00006          11.1
#>  7 DA00007          11.1
#>  8 DA00008          11.1
#>  9 DA00009          11.1
#> 10 DA00010          11.1
#> # ℹ 526 more rows
#> [1] "--------------------------------------"

qc_data$data_summary$na_col_hist

qc_data$data_summary$na_row_hist

qc_data$data_summary$cor_heatmap

qc_data$metadata_summary$na_col_hist

qc_data$metadata_summary$na_row_hist

qc_data$metadata_summary$Sex

qc_data$metadata_summary$Stage

qc_data$metadata_summary$Age

qc_data$metadata_summary$BMI

From the EDA results we can see that all assays have less than 10% missing values, while only a few samples have more than 20% missing values. The protein-protein correlation plot shows that there are no extremely highly correlated proteins. Also, the metadata variable Stage has a high percentage of missing values that are not correctly stated as NA values but as “” or “Unknown”. These are only a few examples of the information that we can extract from the EDA results.

Based on these we can take decisions on filtering samples with high missing values, excluding whole assays and metadata variables from statistical analysis.

Dimensionality Reduction

Finally we will run a PCA and UMAP analysis on the data to check if the data contain any outliers or cluster in an unexpected way. From the PCA analysis we can also check how the variance in the data is explained by the different principal components. In this example, we will color the PCA plot based on Disease, while the UMAP plot based on Age.

The hd_pca() and hd_umap() functions will run the respective analysis while the hd_plot_dim() will plot the data on a 2D plane. For the PCA analysis we can also use hd_plot_pca_loadings() and hd_plot_pca_variance() to plot the loadings of the different proteins on the principal components and the variance explained by each principal component respectively.

In the hd_plot_dim() function we can pass the color parameter to color the data points based on a metadata variable. We should not forget to pass again the hd_obj object to the metadata argument so the package can search the metadata and identify the color variable.

pca_res <- hd_pca(hd_obj, components = 15) |> 
  hd_plot_dim(hd_obj, "PC1", "PC2", color = "Disease", palette = "cancers12", axis_variance = TRUE) |> 
  hd_plot_pca_loadings(displayed_pcs = 6, displayed_features = 10) |> 
  hd_plot_pca_variance()

head(pca_res$pca_res)
#> # A tibble: 6 × 16
#>   DAid     PC1     PC2    PC3    PC4    PC5    PC6     PC7    PC8    PC9    PC10
#>   <fct>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl>  <dbl>  <dbl>   <dbl>
#> 1 DA00… -3.67  -4.28   -2.34  -3.10  -2.65  -2.73  -2.78    0.433 -3.23  -0.399 
#> 2 DA00…  4.11  -2.64    2.04  -0.441 -4.43  -1.91  -0.897   1.54  -0.283  0.166 
#> 3 DA00… -3.34   4.72    1.41   0.881 -0.561  0.308 -0.0612 -0.267  1.42   0.0840
#> 4 DA00… -4.78   0.443   1.41   0.107 -1.10  -0.262  0.350   3.33  -0.336 -0.748 
#> 5 DA00… -4.98  -3.67    0.711 -5.70  -0.807 -3.77  -0.969   1.83  -1.29  -0.991 
#> 6 DA00…  0.395  0.0572 -1.90  -7.75   0.707 -2.70  -0.681   0.301 -0.486  1.32  
#> # ℹ 5 more variables: PC11 <dbl>, PC12 <dbl>, PC13 <dbl>, PC14 <dbl>,
#> #   PC15 <dbl>
pca_res$pca_plot

pca_res$pca_loadings_plot

pca_res$pca_variance_plot

This time we want to color based on a continuous variable, the Age. In this case, we can either bin the data into categories using hd_bin_columns() and use a categorical palette as before or use a continuous color palette directly after the plot is created. We will do the second by setting the palette argument to NULL and using the viridis palette afterwards.

umap_res <- hd_umap(hd_obj, components = 2) |> 
  hd_plot_dim(hd_obj, "UMAP1", "UMAP2", color = "Age", palette = NULL)

umap_res$umap_plot + scale_color_viridis()

We can also run a UMAP where the data points are Proteins (or any feature) instead of samples. The only thing we need to do is to set the by_sample parameter to FALSE and pass the Assay column as the plot_color parameter. We can also pass a custom palette to color the different assays. This time we will use the wrapper function hd_auto_umap(). The wrapper functions are a bit less flexible but they are easier to use and require less code.

umap_res <- hd_auto_umap(hd_obj, 
                         by_sample = FALSE, 
                         plot_color = "Assay",
                         plot_palette = c("ADA" = "darkblue", 
                                          "ABL1" = "red3", 
                                          "ACAN" = "green3"))

head(umap_res$umap_res)
#> # A tibble: 6 × 3
#>   Assay   UMAP1  UMAP2
#>   <fct>   <dbl>  <dbl>
#> 1 AARSD1  3.31  -2.38 
#> 2 ABL1    1.52  -1.85 
#> 3 ACAA1   0.238 -1.50 
#> 4 ACAN   -0.289  1.69 
#> 5 ACE2   -0.307 -0.731
#> 6 ACOX1  -2.58   1.23
umap_res$umap_plot

📓 Remember that these data are a dummy-dataset with artificial data and the results in this guide should not be interpreted as real results. The purpose of this vignette is to show you how to use the package and its functions.