Reproducible workflow for the analysis of CITE-seq data for the project

“Broad immune activation underlies shared set point signatures for vaccine responsiveness in healthy individuals and disease activity in patients with lupus”

Introduction

The following workflow was used to analyze the CITE-seq data and generate the CITE-seq related figures for our paper “Broad immune activation underlies shared set point signatures for vaccine responsiveness in healthy individuals and disease activity in lupus patients” published in Nature Medicine in 2020.

Prepare the working directory

Some steps of the workflow require data generated for the first part of the analysis of flow cytometry and microarrays data (see Workflow.Rmd).

Before running the workflow make sure you have the following subfloders in the working directory:

data
citeseq
--data
--figures
--R
--results
generated_data

• Download the archived CITE-seq data folder from and unpack it into the citeseq/data folder.
• Download the archived directory citeseq with R code and unpack it into the citeseq folder.
• Make sure you have the .Rprofile file in the working directory. Modify the PROJECT_DIR variable with the full path to the working directory.
• Create empty folders generated_data and figure_generation.
• If you are going to use the Singularity container download it from . It should be located outside of the working directory.

R packages used

Our pipeline requires the following R packages:

# CRAN:
install.packages(c("plyr", "tidyverse", "data.table", "pROC", "MASS", "limma", "mclust", "corrplot", "ggraph", "circlize", "pals", "ggsignif", "ggridges", "viridis", "clustree", "cowplot"))

# Bioconductor:
install.packages("BiocManager")
BiocManager::install()
BiocManager::install(c("SingleCellExperiment", "scater", "scran", "fgsea", "tmod", "ComplexHeatmap"))

# Seurat ver. 2.3.4
source("https://z.umn.edu/archived-seurat")

Notes

If while running some scripts you get X11 forwarding related error, check if you have X11 forwarding turned on (in Putty it is in Connection-SSH-X11) and turn it off. Sometime an error may appear due to a conflict between R packages. In this case restart R and rerun the script that previously returned the error.

---

0. Before running this workflow always change the working directory to citeseq.

```r
source(\R/remove_doublet_clusters.r\)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->

<br/>

# 1. Read Seurat object with demultiplexed data and filtered for singlet cells from H1N1 day0 samples only

Input data are stored in ```data``` folder. These files or other RDS files generated on further steps can be downloaded from the [Figshare repository](https://doi.org/10.35092/yhjc.c.4753772). 

  * RDS file with Seurat object
    * ```data/H1_day0_demultilexed_singlets.RDS```
  * RDS file with negative control cells
    * ```data/neg_control_object.rds```


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuc291cmNlKFxcUi9hZHRfY2x1c3RlcmluZ19jbHVzdGVyYXZlcmFnZXNfM2xldmVsc19maWd1cmUuclxcKVxuYGBgXG5gYGAifQ== -->

```r
```r
source(\R/adt_clustering_clusteraverages_3levels_figure.r\)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->

Intermediate RDS files containing lists of Seurat objects with normalized data for individual batches are stored in ```data/normalization_data```.

Output RDS file is stored as ```data/H1_day0_scranNorm_adtbatchNorm.rds```.
<br/>


# 2. Clustering all cells by ADT data at multiple resolution

The clustering was performed using as an input the matrix of distances between cells using ADT data only. 

## Input data

Input data are stored in ```data``` folder:

  * RDS file with Seurat object with RNA-seq data SCRAN-normalized and ADT data batch-normalized and batch-corrected      * ```data/H1_day0_scranNorm_adtbatchNorm.rds```

The calculations require OVER 64gb of RAM and was run on a high-performance cluster. The output file can be downloaded from the [Figshare repository](https://doi.org/10.35092/yhjc.c.4753772).


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuc291cmNlKFxcUi9jb21wdXRlX3BjYS5yXFwpXG5gYGBcbmBgYCJ9 -->

```r
```r
source(\R/compute_pca.r\)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Output: 

  * Table of clusters assignment for each cell: 
    * ```./data/H1N1_full_clustered_p3dist_multires_metadata.rds```.
  * Seurat object with clustering assignment:
    * ```./data/H1_day0_scranNorm_adtbatchNorm_dist_clustered.rds```.
  
<br/>


# 3. Clusters annotation

Generate a heatmap of average expression of selected protein markers in each of the cell clusters derived from the three different clustering resolutions


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuc291cmNlKFxcUi9jb21wdXRlX3RzbmUuclxcKVxuYGBgXG5gYGAifQ== -->

```r
```r
source(\R/compute_tsne.r\)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


The generated table in ```results/cluster_nodes.txt``` is used add cluster labelse and manual annotation of cluster cell types.

The generated figures are stored in ```figures/cluster_annotation``` and used for cluster labeling and annotation.

<br/>


# 4. Filter cells and generate final heatmap

After annotating the clusters we found that at resolution=3 clusters 34, 36, 38, 39 contain mostly doublet cells. We remove cells in these clusters from the Seurat object.


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuc291cmNlKFxcUi9jbHVzdHJlZV9sYWJlbHMuclxcKVxuYGBgXG5gYGAifQ== -->

```r
```r
source(\R/clustree_labels.r\)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Output file:
  ```./data/H1_day0_scranNorm_adtbatchNorm_dist_clustered_filtered.rds```

<br/>

The following script uses manual annottion table in ```data/clustree_node_labels_withCellTypeLabels.txt``` and clusters at first three resolutions to generate Figure 4c.


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuc291cmNlKFwiUi9hZHRfY2x1c3RlcmluZ19jbHVzdGVyYXZlcmFnZXNfM2xldmVsc19maWd1cmUuclwiKVxuYGBgIn0= -->

```r
source("R/adt_clustering_clusteraverages_3levels_figure.r")

The generated heatmap is stored in figures/cluster_annotation.

5. PCA and tSNE analysis of the ADT data

source("R/compute_pca.r")

Output file: ./data/H1_day0_scranNorm_adtbatchNorm_dist_clustered_PCA.rds

PCA elbow plot is stored in figures/TSNE. Using this plot We determined the number of principal components to use for tSNE as 7.

source("R/compute_tsne.r")

This script for tSNE computation at multiple perplexity was run on high-performance cluster. The output file can be downloaded from the Figshare repository.

Output file: ./data/H1_day0_scranNorm_adtbatchNorm_dist_clustered_TSNE.rds

6. Label and visualise the clusters

Re-label the clusters and create the final object.

source("R/clustree_labels.r")

Output file: ./data/H1_day0_scranNorm_adtbatchNorm_dist_clustered_TSNE_labels.rds. The output file can be downloaded from the Figshare repository.

This Seurat object was used for all further analysis.

Generate tSNE plots to check for batch effect and show the first-level clusters for Figure 4a

source("R/tsne_figures.r")

The generated figures are stored in figures/TSNE.

Generate the clustree for Figure 4b with nodes colored by level 1 cluters

source("R/clustree_vertical_clr_level1.r")

The generated figure is stored in figures/cluster_annotation.

Generate the supplemental figure with the distribution of denoised and background rescaled counts for the surface proteins

source("R/adt_clustering_protein_distribution.r")

The generated figures are stored in figures/cluster_annotation.

7. Test various gene signatures comparing low and high responders

Generate the list of signatures:

source("R/make_list_of_signatures.r")

Output file with the list of signatures is saved as sig/sig.list.RDS.

Check enrichemnt of BTM modules in CD40act gene signature:

source("R/CD40act_genes_HG.test.r")

Output plot is saved in figures/CD40act.

Test gene signatures comparing low vs. high responders using cells in pseudo-bulk and all cell clusters:

source("R/test_sigs_clusters.r")

Output file with the list of signatures is saved as results/test_sig.genes_in_clusters_high_vs_low_responders.txt.

Compute scores for selected gene signatures in pseudo-bulk and all cell clusters:

source("R/sig_scores_calc_by_cluster.r")

Output files are saved in results/sig_scores.

8. Visualize the results of gene signatures tests

Generate box-plot for pseudo-bulk for gene signatures with significant difference between low and high responders:

source("R/sig_scores_boxplot_pseudobulk.r")

Output plots (for Figure 4d,e) are saved in results/sig_test_boxplots.

Generate box-plot for level 1 clusters for gene signatures with significant difference between low and high responders:

source("R/sig_scores_boxplot_by_clusters.r")

Output plots (for Fig. 5a,b,e and Ext. Data Fig. 8f) are saved in results/sig_test_boxplots.

Generate clustree plots for selected gene signatures indicating significance of difference between low and high responders:

source("R/clustree_vertical_test_sigs.r")

Output plots (for Fig. 5a,b,e and Ext. Data Fig. 8f) are saved in results/sig_test_clustree.

Generate box-plots for cluster C9 (pDC) for selected gene signatures:

source("R/sig_scores_boxplot_C9.r")

Output plots (for Fig. 5c) are saved in results/sig_test_boxplots.

Drop-out test for selected gene signatures and selected cluster combinations:

source("R/test_sigs_clusters_OUT.r")

Result table is saved as results/test_sig.genes_in_clusters_high_vs_low_responders_OUT.txt.

Output plots (for Ext. Data Fig. 8g) are saved in figures/drop_out_test.

9. Manual gating of CD20+CD38++ B cells

source("R/hand_gating.r")

Output table of cell frequencies is saved as results/CITEseq_CD38hi_cell_data.txt. Output plots (for Ext. Data Fig. 8d,e) are saved in results/hand_gating.

10. Compile scores for multiple signatures and perform correlation analysis

Compute scores of SLE-Sig signature for microarray day 0 data:

source("R/SLE.sig_MA_sig_score.r")

Output file is saved as results/sig_scores/scores_SLE.sig_MA.txt.

Compile scores for selected gene signatures for CITE-seq and microarray data plus CD38++ cell frequency from flow cytometry data:

source("R/scores_for_correlation.r")

Output table of cell frequencies is saved as results/CITEseq_CD38hi_cell_data.txt.

Generate scatter plot of CD40act score in C3.1.0 cell cluster vs. CD20=CD38++ B cells frequency from flow cytometry (Figure 5f):

source("R/CD40act_vs_CD38hi.flow.r")

Output plot is saved in figures/sig_ranked_correlation.

Generate scatter plots of correlation between selected signature scores in selected clusters (for Figure 6a):

source("R/sig_correlation_scatterplots.r")

Output plots are saved in figures/sig_ranked_correlation.

Generate a correlation plot (for FIgure 6a):

source("R/scores_correlation_heatmap_custom_order.r")

Output plot is saved in figures/sig_ranked_correlation.

11. Differential Expression of selected surface proteins in pDC between low and high responders

source("R/CD86_HLA-DR_vs_response.r")

Output plot is saved in figures/adt_vs_response.

12. Comparing the TGSig and SLE-Sig scores in females versus males

source("R/TGSig_SLE.sig_vs_sex.r")

Output plot is saved in figures/male_vs_female.