Functional dysconnectivity of visual and somatomotor networks yields a simple and robust biomarker for psychosis

2.1. Data

The data was collected from the 10x genomics website (https://www.10xgenomics.com/datasets/frozen-pbm-cs-donor-c-1-standard-1-1-0) and sequenced using 10x genomics chromium technology (Home Page - 10x Genomics, n.d.. In this technique, the Gel bead in Emulsion (GEM) method is used to put single cells in the shape of droplets (which is why this is one of the droplet-based sequencing techniques), labelling each gel bead with oligonucleotides unique to each bead and having a length of 10 bp, also known as barcodes or UMIs (unique molecular identifiers). Additionally, they have sequencing primers and adapters and an immobile attached oligo-dT of 30 bp.

The datasets were presented as a raw counts sparse grid array. (pre-analyzed, primary analysis done using Cell Ranger 1.1.0), with 9519 cells in the column and 32738 genes as rows.

2.2. Analysis

Figure 1.

The pipeline for the analysis of the scRNA seq derived PBMC data.

2.3. Pre-processing and HVG selection

The data was first “pre-processed,” which means it was TPM-normalized and log-transformed. Then, the highly variable gene selection method was used to find the most informative genes. For the clustering step, only those informative genes which showed maximum variation across cells in the data were chosen. So, the gene selection step not only picked the genes that changed a lot, but it also eliminated genes that were present in very few cells (because they were caught by chance or noise). Elimination of the housekeeping genes and the genes which were expressed uniformly throughout was done at this stage as they won’t be useful in the following step of cluster separation.

“Semi-soft clustering of single-cell data (SOUP)” (Zhu et al., 2018) was used for gene selection and clustering steps. This tool provides two methods for selecting highly variable and informative genes. First is the sparse PCA algorithm (SPCA), which is a lasso (elastic net)-based PCA method specialized for sparse loadings such as single-cell datasets (Witten et al., 2009). PCA is used for dimensionality reduction of the data by selecting a smaller number of variables that retain most of the sample’s information. It often happens that in the case of classic PCA, the input variables are typically combined linearly to form the principal components. This makes the result unpredictable and makes it hard to determine what the principal components mean since they are all not zero. This issue is resolved with sparse PCA by looking for linear combinations that use only a few input variables. Thus, in the case of gene expression data, each variable corresponds to a specific gene. With only a small number of genes, sparse PCA can generate a principal component so that the analysis can be narrowed down to just these genes. This method is better suited for single-cell datasets because it finds linear combinations while efficiently utilizing a small number of input variables, unlike the classic PCA method, which is typically a linear combination of all input variables. Also, the sparse PCA can retain consistency, which the classic PCA cannot.

The second method of gene selection to reduce noise and filter out less informative genes is the DESCEND algorithm (Wang et al., 2018). DESCEND (“deconvolution of single-cell expression distribution”, Wang et al., 2018) utilizes Gini coefficients to provide accurate estimations of measures based on distributions, such as the actual gene expression variance and the likelihood that true gene expression is positive. While preserving the genes that are more plausible to serve as cell subtype identifiers, low-variance genes are filtered away to reduce noise. The coefficient can take on values between 0 and 1, with 0 denoting the least possible variation and 1 denoting highly variable expression throughout. Thus, the greater the Gini index, the more distributed the data.

Together, these two methods choose a group of informative genes that shows high variance in their expressions in different cells and are also informative. These genes are then further subjected to clustering.

2.4. Clustering

The next step is clustering the cells using the selected genes. SOUP is a method for semi-soft clustering scRNA datasets that fills a void in characterizing cells with ambiguous labels, such as those that transition from one type of cell to another (Zhu et al., 2018). The SOUP procedure begins with identifying a collection of pure cells, followed by estimating the non-negative membership matrix. Considering that the data has (i) Pure cells belonging to a common group that requires a hard cluster allocation and (ii) Developmental cells (transitional cells) that are shifting among different cell types and should be assigned lenient allocations, SOUP first finds K different clusters (or cell subsets) from data and then assigns membership scores to mixed cells for each cluster K, forming a “membership matrix”. The membership matrix’s every row has positive numbers that add up to one. These numbers show how many cells are in each of the K clusters. Specifically, a type k non-transitional cell has a membership score of 1 and zeros elsewhere. Here, hard clustering is utilized for this data. An elbow plot (optimal K) was used to figure out K, and Scater (McCarthy et al., 2017) was used to make PCA, TSNE (t-distributed stochastic neighbour embeddings), and UMAP (Uniform manifold approximation and projection) plots. These plots are typically used to visualize single-cell data in low-dimensional space. They are also known as “non-linear dimensional reduction methods.”

2.5. Differential Expression Analysis of the Clusters

After the clustering, edgeR was used for in-depth analysis of the clusters In edgeR analysis, each cluster (including all of the cell types in the cluster) was treated as a sample, similar to other conventional differential expression analyses. The Differential analysis was carried out to determine potential key markers for each cluster. The contrasts were made to compare one cluster and the mean of all the remaining clusters. This procedure yielded a list of the most significant genes for each cluster. Then, various databases like the Human Protein Atlas, the CellMarker Database, and the EBI were used to find the standard markers of immune cells. The top genes from each cluster were searched for in the above-mentioned databases. According to the search, the cell subsets were assigned to the cluster based on which immune cell their top genes were markers for.

2.6. Construction of the Signature Matrix

After obtaining the matrix of cell types and their marker genes, the tool CIBERSORTx was implemented to construct a custom signature matrix. A .txt or .tsv file containing a tab-delimited matrix of scRNA-seq expression data, with rows representing genes and columns representing single-cell transcriptomic data, makes up the input file. It was checked to make sure the sum of all the genes expressed in a given cell wasn’t zero. The first column of the file was filled with gene names. Redundant genes were removed. Also, each cell in row 1 was given a cell phenotype, such as “CD8 T cell” or “B cell,” and it was made sure that there were at least three cells of each phenotype. The highly expressed G marker genes from every cell subset were integrated into a signature matrix BG, and CIBERSORT sorted significant genes for each cell type subset in decreasing order of fold change in comparison to other cell subsets. The algorithm preserved the signature matrix, identifying the one with the minimum condition number after iterating G from 50 to 200 across all subgroups. To avoid genes expressed on non-hematopoietic cell types confusing the findings of deconvolution, the option to filter out non-hematopoietic genes were activated. After we obtained the signature matrix for the immune cells, visualizations were drawn by means of a heatmap and hierarchical clustering to ensure the correct groupings and sub-groupings.

3.1. Pre-processing

The downloaded PBMC data had 9519 cells in the column and 32738 genes as rows. Since the data for these genes were raw counts, it was TPM normalized and log transformed. To find out the existing grouping in this data, single-cell linear and non-linear dimensionality reduction approaches were used. We utilized PCA, TSNE, and UMAP plots to view the data in lower dimensions (figure 2.).

Figure 2.

TSNE, PCA, and UMAP of PBMC data. These plots show that the data is divided into at least eight distinct groups.

3.2. HVG selection

Single-cell data usually requires clustering of the samples to find out the structure present in them, but before clustering, it is important to filter out the non-informative genes. In the first stage, we removed the sparsely expressed genes or genes whose expression was detected in less than ten cells. Then DESCEND, and SPCA was applied to find out highly variable genes from the data. The DESCEND selected 609 genes, and the SPCA selected 158 genes using default parameters.

Fig. 3. Shows the distribution and variance of the genes. The curve was made using quantile smooth regression (the default quantile is 0.5), demonstrating the correlation between the mean of deconvolved distributions and the Gini coefficient of the genes. The Gini distance from the curve is used to figure out each gene’s dispersion score. This score is normalized by the gene’s SE (squared error). Genes with a Gini coefficient and a normalized score higher than a threshold (the default is 10) are marked as having a lot of variation as shown in red. Total of 678 genes out of 32738 were identified as highly variable and the analysis was further continued with this set of highly variable genes only (along with all the cells). SOUP selected 678 genes when the gene selection was performed using default parameters. However, the number of selected genes may be increased or decreased when the parameters such as “threshold” and “sumabs” are changed in the gene selection step. The “threshold” is the cutoff for the Gini index of DESCEND. Fewer genes are picked when the threshold is higher. Similarly, “sumabs” is a measurement of the sparsity of genes in SPCA, between 1/sqrt(number of genes) and 1. Smaller values result in sparser results, hence fewer selected genes.

Figure 3.

HVG selected by DESCEND. The distribution and variance of the genes. The curve was made using quantile smooth regression (the default quantile is 0.5), showing the relationship between the mean of deconvolved distributions and Gini across genes. The Gi distance from the curve is used to figure out each gene’s dispersion score. This score is normalized by the gene’s squared error). Genes with a Gini coefficient and a normalized score higher than a threshold (the default is 10) are marked as having a lot of variation. The data points in red show this. Out of 32738 genes, 678 were designated as highly variable.

3.3. Clustering

Clustering techniques were performed as a next step after getting highly variable genes. SOUP allows users to specify the value of K for clustering. For the identification of optimal K in our data, we considered the elbow plot. To choose the ideal K, the number of clusters, SOUP offers a cross-validation process. Specifically, using 10-fold cross-validation, users can search across K > 2, > 10, (or more), and the mistakes are averaged over nCV repeats (Zhu et al., 2018). Cross-validation errors were visualized in the log scale when comparing several K values and discovered that K = 8 had the lowest cross-validation error. This suggests that there are at least eight groups in our PBMC data. The selected highly variable genes are able to segregate the data into 8 distinct groups. If the number of genes selected as highly variable is increased or decreased, then, maybe the number of groupings in our data also changes because when the data was visualized in PCA, t-SNE, and UMAP more than 8 groupings of cells were found. These other groupings may be subtypes of cells for example activated cells or resting cells and so on.

Figure 4.

Elbow plot showing the optimal K in the PBMC data. The elbow falls the lowest cross-validation error to 8, which shows that our data has at least 8 clusters.

Therefore, we continued to cluster the data by increasing K until we obtained clusters with fewer than five cells. Thus, we tested K = 10 and increased K until this condition was met, stopping at K = 15.

We further proceeded with K = 15 results, and the clusters were annotated on the above visualizations.

Figure 5.

TSNE, PCA, and UMAP plot showing clusters identified in the PBMC data plots showing clusters identified in PBMC data.

3.4. Differential expression analysis of the clusters

After the clustering methods r, differential expression analysis was done on the clustered groups. The goal was to find the top significant genes that could separate a cluster from the rest of the clusters. For this purpose, we considered using the classic differential expression analysis tool, edgeR (Robinson et al., 2010). EdgeR favors the raw integer read counts and finds a scaling (normalization) factor for every sample to normalize the read counts for different library sizes (sample-specific effect) (Robinson et al., 2010). The TMM (Trimmed Mean of M-values) (Robinson et al., 2010) between-sample normalization approach was used for the normalization. TMM presupposes that most of the genes do not express themselves differently. The scaling factors were then used to calculate the actual library sizes. CR or Cox-Reid profile-adjusted likelihood (CR) is used in the EdgeR approach to estimate dispersion for trials with many components. After getting estimates of dispersion and fitting negative binomial generalized linear models, we moved on to the testing steps to find differential expressions using the likelihood ratio test.

A design matrix is necessary for the GLM technique in order to specify the treatment conditions. To create the design matrix, we used the “model.matrix” function. From the groupings of clusters, we made a model matrix in which each cell had a 1 under the cluster where it was found and a 0 everywhere else, as shown below: Using the QL F-test, we looked for genes that showed significant differential expression. The “makeContrasts” function was used to specify the contrast of interest. The comparisons were set up so that one cluster was put up against the average of every other cluster, as shown below: This testing process provided the significant genes for each cluster (positively and negatively expressed) (Table 2. top 10 significant genes from each cluster).

Figure 6.

Heatmap of top 10 significant genes from each cluster and their differential expression across the clusters.

3.5. Cell type classification

We obtained the top significant genes for each cluster. These were positively or negatively expressed in the clusters, but only positive and significantly expressed genes were considered for classifying the cell types. For reference, we used databases such as the CellMarker database (Hu et al., 2022), The human protein atlas (The Human Protein Atlas, n.d., Karlsson et al., 2021, Uhlen et al., 2019), and EBI-Single cell expression atlas (Home < Single Cell Expression Atlas < EMBL-EBI, n.d.). The top genes were searched for cell-type-specific expression in these databases. The next table (Table 3.)shows the canonical markers chosen from each cluster to show that they are a certain type of immune cell.

Figure 7.

Dotplot of canonical markers and their differential expression across the clusters.

Cluster 15: Activated/Matured B cells

This cluster had high expression of B cell-specific markers like CD79A, CD79B, MS4A1, BANK1, maturation marker CD27, activation markers like CD22, FCER2, FCRLA, RALGPS2, and antigen-presenting markers like HLA-DQB1, HLA-DQA1, HLA-DOB, and HLA-DQA2. Therefore, this group was classified as Activated/Matured B cells.

The figure 8. shows final annotation of the clusters identified, into the cell types classified using the canonical markers. The groups and their subgroups fall closer in the plots and overlap. For example, both NK cell clusters are shown close together in the plot. The B cell subsets which are the memory B cells and activated B cells form another close group. We also observed the myeloids tightly bound together and similarly, all the T cells fall together closely in the above plots. The plasmacytoid DC although a myeloid cell is still similar to Plasma cells which are matured B cells, therefore this cell type cluster falls closer to the B cell group.

Figure 8.

t-SNE, PCA and UMAP plots showing the cell types identified in the PBMC data. Note that the T cells are grouped together; similarly, the B cells are also grouped together; the CD16 monocyte is a distinct cluster, as is the one with NK cells. All the myeloid cells, such as DCs, monocytes, macrophages, and neutrophils, are grouped closely. Another thing that stands out is that plasmacytoid DCs are close to B cells, which makes sense because plasmacytoid DCs are similar to plasma cells.

3.6. Construction of the Signature Matrix

The signature matrix was built after classifying the clusters as the types of cells with genes that were highly variable across the cells. We had 14 cell types and 678 genes and constructed the signature matrix by the method mentioned above. The reference file was prepared to have each cell labelled as the cell type it belonged to. Non-hematopoietic genes were set to be filtered from the resultant matrix. We obtained a signature matrix with 14 columns representing 14 cell types and 275 genes. The matrix had values for canonical markers the highest for their respective cell types (up to thousands) and 1 everywhere else. For example, The B cell marker, MS4A1, was the highest in the two B cell subtypes (memory B cells and activated/matured B cells) but the least (1) for other cell types. Similarly, the memory markers for B cells were high compared to those in activated B cells, but they were the least, i.e., 1 for every other cell type. The same was true for other cell types and subtypes as well.

Figure 9. shows a heatmap for the differential expression of the canonical markers in the custom signature matrix. It also shows the hierarchical clustering of the cell types. Heatmap shows the myeloid group separated from the lymphoid group. Under the Myeloid group fall the monocytes, DCs and macrophages. It also shows the macrophages and Monocyte DCs grouped together which is correct as they both are antigen presenting cell and originate from monocytes. Similarly, the neutrophils and monocytes grouped together since they both are phagocytes that act against infections. They also share some of the markers due to similar type of functions, thus, this grouping holds. The plasmacytoid DCs in the myeloid group but separate from other myeloid cell types. This is because although myeloid, the plasmacytoid DCs have similarities with the plasma cells (or matured B cells), they fall in the myeloid group because they have antigen-presenting markers, a characteristic of dendritic cells and other myeloid and dendritic cell markers.

Figure 9.

Heatmap showing the differential expression of the canonical markers in the custom signature matrix.

In the lymphoid group, we see three different sections. One is for B cells described earlier, the otehr is T cells and one more which is the NK cells group. The NK cell group has both the NK cell subtypes which holds true for obvious reasons. The T cell group has the cytotoxic and effector T cells grouped together and both the naive groups together. This grouping is true because the effector T cells are also capable of killer functions just like the cytotoxic T cells, thus they have similar expression profiles and many common markers. The naive group has both the naive cell forms of T cells.

SOUP selected 678 genes when the gene selection was performed using default parameters. However, the number of selected genes may be increased or decreased when the parameters such as “threshold” and “sumabs” are changed in the gene selection step. The “threshold” is the cutoff for the Gini index of DESCEND. Fewer genes are picked when the threshold is higher. Similarly, “sumabs” is a measurement of the sparsity of genes in SPCA, between 1/sqrt(number of genes) and 1. Smaller values result in sparser results, hence fewer selected genes.

The 678 genes are able to segregate the data into 8 distinct clusters. If, by changing the parameters in the gene selection step, the number of highly variable genes is increased or decreased, the number of groupings may change when identified by an elbow plot. However, it won’t really change the overall grouping structure in the 3 visualizations (PCA,t-SNE, UMAP). Thus, to sum up, the number of genes selected for further clustering the cells, play a crucial role in cluster segregation.

Since there were more than 8 groupings in the visualization plots (PCA, t-SNE, UMAP). These groupings could be cell subtypes, which is why the clustering was further continued till K15, after that, there were clusters with less than five cells. Therefore, it was stopped at k15.