Tutorial: Single Cell RNA-Seq Analysis of NK, CD4+ and CD8+ T cells isolated from LCMV infected WT and Ifng - CTCF binding site mutant mice, using SCassist - an AI-Based Workflow Assistant
1. Quality Control
This document takes the cellranger count data for two samples (a WT and an Ifng CBS mutant mice) and performs standard single cell rnaseq analysis, using the SCassist package.
# Install hdf5 library if its not already installed
#sudo apt-get install libhdf5-dev or brew install hdf5
#install.packages("hdf5r")
# Load below libraries
library(Seurat)
library(ggplot2)
library(plotly)
library(tidyverse)
library(cowplot)
library(scProportionTest)
library(rollama)
library(SCassist)
library(httr)
library(jsonlite)1.1. Data setup
# Set the input, output paths
setwd("/singlecellassistant/exampleOutput_G/LCMV/")
#datapath="exampleData/LCMV"
outputpath="/singlecellassistant/exampleOutput_G/LCMV/"
# Load RData for example dataset
load("/singlecellassistant/exampleOutput_G/LCMV/LCMV_G.RData")
# Specify samples for analysis.
# Define LCMV Day4, CD4, CD8, NK samples
samples=c("WT","KO")
# Read the data file
WT <- Read10X_h5(file.path("/singlecellassistant/exampleData/LCMV/GSM6625298_scRNA_LCMV_Day4_CD4_CD8_NK_WT_filtered_feature_bc_matrix.h5"), use.names = T)
KO <- Read10X_h5(file.path("/singlecellassistant/exampleData/LCMV/GSM6625299_scRNA_LCMV_Day4_CD4_CD8_NK_KO_filtered_feature_bc_matrix.h5"), use.names = T)
# Add sample names to columnnames
colnames(WT$`Gene Expression`) <- paste(sapply(strsplit(colnames(WT$`Gene Expression`),split="-"),'[[',1L),"WT",sep="-")
colnames(KO$`Gene Expression`) <- paste(sapply(strsplit(colnames(KO$`Gene Expression`),split="-"),'[[',1L),"KO",sep="-")
# Create a Seurat data object from the gex matrix
WT <- CreateSeuratObject(counts = WT[["Gene Expression"]], names.field = 2,names.delim = "\\-")
KO <- CreateSeuratObject(counts = KO[["Gene Expression"]], names.field = 2,names.delim = "\\-")
# Merge seurat objects
allsamples <- merge(WT, KO, project = "LCMV")
allsamples## An object of class Seurat
## 32285 features across 10551 samples within 1 assay
## Active assay: RNA (32285 features, 0 variable features)
## 2 layers present: counts.1, counts.2head(allsamples)## orig.ident nCount_RNA nFeature_RNA
## AAACCCAAGAGTTGTA-WT WT 17880 4120
## AAACCCAAGATTTGCC-WT WT 7916 2395
## AAACCCAAGTTGTAGA-WT WT 35993 5094
## AAACCCACAAGTATCC-WT WT 9059 2865
## AAACCCAGTCCAGGTC-WT WT 10948 2939
## AAACCCATCCTGGGAC-WT WT 16768 3484
## AAACGAACAGAGAGGG-WT WT 10467 3159
## AAACGAACAGCGACAA-WT WT 12996 2853
## AAACGAACAGTGTGGA-WT WT 26674 4654
## AAACGAATCTAGCCTC-WT WT 11434 3100table(Idents(allsamples))##
## WT KO
## 5666 4885# Merge counts data
allsamples = JoinLayers(allsamples)1.2. Quality Analysis
This part of the workflow uses SCassist_analyze_quality(), to identify potential filtering cutoff values. SCassist provides the recommendation based on the summary statistics and the quantile statistic of the data.
# Generate Percent mito for each of the cells
grep("^mt-",rownames(allsamples@assays$RNA$counts),value = TRUE)## [1] "mt-Nd1" "mt-Nd2" "mt-Co1" "mt-Co2" "mt-Atp8" "mt-Atp6" "mt-Co3"
## [8] "mt-Nd3" "mt-Nd4l" "mt-Nd4" "mt-Nd5" "mt-Nd6" "mt-Cytb"allsamples$percent.mito <- PercentageFeatureSet(allsamples, pattern = "^mt-")
summary(allsamples$percent.mito)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.819 2.425 5.053 3.310 83.600head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGAGTTGTA-WT WT 17880 4120 3.333333
## AAACCCAAGATTTGCC-WT WT 7916 2395 2.627590
## AAACCCAAGTTGTAGA-WT WT 35993 5094 2.008724
## AAACCCACAAGTATCC-WT WT 9059 2865 1.953858
## AAACCCAGTCCAGGTC-WT WT 10948 2939 2.831567
## AAACCCATCCTGGGAC-WT WT 16768 3484 0.942271
## AAACGAACAGAGAGGG-WT WT 10467 3159 3.668673
## AAACGAACAGCGACAA-WT WT 12996 2853 1.569714
## AAACGAACAGTGTGGA-WT WT 26674 4654 1.570818
## AAACGAATCTAGCCTC-WT WT 11434 3100 1.687948# Generate Percent ribo for each of the cells
allsamples$percent.ribo <- PercentageFeatureSet(allsamples, pattern = "^Rp[sl]")
summary(allsamples$percent.ribo)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.002 18.156 24.731 24.894 32.868 56.275head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGAGTTGTA-WT WT 17880 4120 3.333333
## AAACCCAAGATTTGCC-WT WT 7916 2395 2.627590
## AAACCCAAGTTGTAGA-WT WT 35993 5094 2.008724
## AAACCCACAAGTATCC-WT WT 9059 2865 1.953858
## AAACCCAGTCCAGGTC-WT WT 10948 2939 2.831567
## AAACCCATCCTGGGAC-WT WT 16768 3484 0.942271
## AAACGAACAGAGAGGG-WT WT 10467 3159 3.668673
## AAACGAACAGCGACAA-WT WT 12996 2853 1.569714
## AAACGAACAGTGTGGA-WT WT 26674 4654 1.570818
## AAACGAATCTAGCCTC-WT WT 11434 3100 1.687948
## percent.ribo
## AAACCCAAGAGTTGTA-WT 22.80761
## AAACCCAAGATTTGCC-WT 16.28348
## AAACCCAAGTTGTAGA-WT 21.85147
## AAACCCACAAGTATCC-WT 16.81201
## AAACCCAGTCCAGGTC-WT 25.30142
## AAACCCATCCTGGGAC-WT 35.07276
## AAACGAACAGAGAGGG-WT 22.92921
## AAACGAACAGCGACAA-WT 42.22068
## AAACGAACAGTGTGGA-WT 26.22029
## AAACGAATCTAGCCTC-WT 28.44149# Generate Percent hemoglobin for each of the cells
allsamples$percent.hb <- PercentageFeatureSet(allsamples, pattern = "^Hb[^(p)]")
summary(allsamples$percent.hb)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.004744 0.005529 2.162162head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGAGTTGTA-WT WT 17880 4120 3.333333
## AAACCCAAGATTTGCC-WT WT 7916 2395 2.627590
## AAACCCAAGTTGTAGA-WT WT 35993 5094 2.008724
## AAACCCACAAGTATCC-WT WT 9059 2865 1.953858
## AAACCCAGTCCAGGTC-WT WT 10948 2939 2.831567
## AAACCCATCCTGGGAC-WT WT 16768 3484 0.942271
## AAACGAACAGAGAGGG-WT WT 10467 3159 3.668673
## AAACGAACAGCGACAA-WT WT 12996 2853 1.569714
## AAACGAACAGTGTGGA-WT WT 26674 4654 1.570818
## AAACGAATCTAGCCTC-WT WT 11434 3100 1.687948
## percent.ribo percent.hb
## AAACCCAAGAGTTGTA-WT 22.80761 0.011185682
## AAACCCAAGATTTGCC-WT 16.28348 0.000000000
## AAACCCAAGTTGTAGA-WT 21.85147 0.002778318
## AAACCCACAAGTATCC-WT 16.81201 0.000000000
## AAACCCAGTCCAGGTC-WT 25.30142 0.009134088
## AAACCCATCCTGGGAC-WT 35.07276 0.005963740
## AAACGAACAGAGAGGG-WT 22.92921 0.009553836
## AAACGAACAGCGACAA-WT 42.22068 0.000000000
## AAACGAACAGTGTGGA-WT 26.22029 0.000000000
## AAACGAATCTAGCCTC-WT 28.44149 0.000000000# Check the total number of cells in this experiment
table(allsamples$orig.ident)##
## KO WT
## 4885 56661.3. Quality Filtering
Analyze quality of the cells using SCassist and use the SCassist recommended QC cutoff values obtained in the previous step to filter and visualize the before and after QC data.
# Generate before filtering quality plots
qcbefore=VlnPlot(allsamples,features = c("nFeature_RNA", "nCount_RNA","percent.mito","percent.ribo"),ncol = 2, pt.size = 0) +
NoLegend()
qcbefore
# Save the plot in current working directory
ggsave(paste(outputpath,"LCMV-before-qc-violineplot.pdf",sep=""), plot = qcbefore, width = 15, height = 20, units = "cm")
api_key_file = "/singlecellassistant/api_keys.txt"
# Analyze quality of the cells using SCassist
allsamplesquality <- SCassist_analyze_quality("allsamples", percent_mt = "percent.mito", percent_ribo = "percent.ribo", percent_hb = "percent.hb", api_key_file=api_key_file, llm_server="google")## Based on the data summary, below are my recommendations for the quality filtering of the data:
##
## **nCount_RNA:**
##
## * **Lower Cutoff:** 1500. This value is chosen to be slightly above the 5th percentile (947) to remove cells with very low counts, but still capture a significant portion of the data.
## * **Upper Cutoff:** 35000. This value is chosen to be slightly below the 95th percentile (26528) to remove cells with extremely high counts, which could indicate potential doublets or other artifacts.
##
## **nFeature_RNA:**
##
## * **Lower Cutoff:** 700. This value is chosen to be slightly above the 5th percentile (528) to remove cells with very few detected genes, but still capture a significant portion of the data.
## * **Upper Cutoff:** 5500. This value is chosen to be slightly below the 95th percentile (4793) to remove cells with an unusually high number of detected genes, which could indicate potential doublets or other artifacts.
##
## **percent.mito:**
##
## * **Upper Cutoff:** 15. This value is chosen to be significantly lower than the 95th percentile (23.01) to remove cells with a high percentage of mitochondrial reads, which could indicate cell stress or damage.
##
## **percent.ribo:**
##
## * **Upper Cutoff:** 45. This value is chosen to be slightly lower than the 95th percentile (41.97) to remove cells with a high percentage of ribosomal reads, which could indicate cell stress or damage.
##
## **percent.hb:**
##
## * **Upper Cutoff:** 0.025. This value is chosen to be slightly higher than the 95th percentile (0.019) to remove cells with a high percentage of hemoglobin reads, which could indicate contamination from blood cells.
##
## **Important Note:** These are just recommendations, and the researcher should test a range of values around these recommendations to determine the optimal cutoffs for their specific dataset. The optimal cutoffs will depend on the specific characteristics of the data and the research question being addressed.# Filter data based on the above plot.
allsamplesgood <- subset(allsamples, subset = nCount_RNA > 1500 & nCount_RNA < 35000 & nFeature_RNA > 700 & nFeature_RNA < 5500 & percent.mito < 15 & percent.ribo < 45 & percent.hb < 0.025)
# Counts before filtering
table(allsamples$orig.ident)##
## KO WT
## 4885 5666# Counts after filtering
table(allsamplesgood$orig.ident)##
## KO WT
## 3979 45761.4. Generate after QC plots
# Generate after filtering quality plots
qcafter=VlnPlot(allsamplesgood,features = c("nFeature_RNA", "nCount_RNA","percent.mito","percent.ribo"),ncol = 2, pt.size = 0) +
NoLegend()
qcafter
# Save the plot in current working directory
ggsave(paste(outputpath,"LCMV-after-qc-violineplot.pdf",sep=""), plot = qcafter, width = 15, height = 20, units = "cm")2. Normalization
Here we ask SCassist_recommend_normalization to use the number of cells, mean gene expression, standard deviation of gene expression, library size and the coefficient of library size variation, to analyze the data properties and recommend an appropriate normalization method among “LogNormalize”, “CLR”, “RC”, “SCTransform”, along with its reasoning for the recommendation.
# Identify suitable normalization method
normalization_recommendation<-SCassist_recommend_normalization("allsamplesgood", api_key_file = api_key_file, llm_server="google")## ## Recommended Normalization Method: SCTransform
##
## Based on the provided characteristics of the single-cell RNA-seq dataset, **SCTransform** is the most appropriate normalization method for this dataset. Here's why:
##
## **1. Large Number of Cells:** SCTransform is designed to handle large datasets efficiently. With 8555 cells, computational efficiency becomes a significant factor, and SCTransform excels in this regard.
##
## **2. High Variability in Gene Expression:** The high standard deviation of 14.8757353917764 in gene expression values indicates significant variability across cells. SCTransform accounts for this variability by modeling the technical noise associated with each gene and cell, leading to more accurate normalization.
##
## **3. Moderate Library Size Variation:** The coefficient of variation of 0.538648217306346 for library sizes suggests moderate variation. While not extremely high, SCTransform effectively handles library size variations by regressing out the effect of library size during normalization.
##
## **Why SCTransform is Preferred:**
##
## * **Comprehensive Normalization:** SCTransform combines normalization, variance stabilization, and feature selection into a single step, making it a powerful and efficient approach.
## * **Cell-Specific Modeling:** It models the technical noise associated with each cell and gene, leading to more accurate normalization and downstream analysis.
## * **Improved Downstream Analysis:** SCTransform has been shown to improve the performance of downstream analyses like dimensionality reduction, clustering, and differential gene expression analysis.
##
## **Potential Alternatives:**
##
## * **LogNormalize:** This method is a simple and commonly used normalization technique. However, it does not account for cell-specific variability or library size variations as effectively as SCTransform.
## * **CLR (Centered Log Ratio):** CLR is a robust normalization method that can handle zero values. However, it is not as widely used as SCTransform and may not be as effective in handling high variability in gene expression.
## * **RC (Relative Counts):** RC is a simple normalization method that divides each gene count by the total number of counts in the cell. It is not as sophisticated as SCTransform and may not be suitable for datasets with high variability.
##
## **Conclusion:**
##
## While other normalization methods are available, SCTransform is the most recommended option for this dataset due to its ability to handle large datasets, high variability in gene expression, and moderate library size variations. It provides a comprehensive and efficient approach to normalization, leading to improved downstream analysis.# Normalize Data
# Below we increase memory to 1GB for SCTransform run
options(future.globals.maxSize = 2000 * 1024^2)
# SCTransform normalization based on SCassist recommendation
allsamplesgood <- SCTransform(allsamplesgood, vars.to.regress = c("percent.mito", "nFeature_RNA", "percent.ribo", "nCount_RNA"))3. Variable features
Explore the variable features using the LLMs abilities to glean insights in to a list of genes, to explore what drives the cell-to-cell variation. Here we ask SCassist_analyze_variable_features to explore the top 30 variable features, along with the experimental design statement and provide potential insights on the underlying system being studied.
# Experimental design statement
experimental_design = "NK, CD4+ and CD8+ T cells isolated from lymphocytic choriomeningitis virus infected WT and Ifng - CTCF binding site mutant mice"
# List top 10 variable genes
top10 <- head(VariableFeatures(allsamplesgood), 10)
top10## [1] "Gzma" "Ccl5" "Hist1h2ap" "Hist1h1b" "Ccl4" "Gzmb"
## [7] "Hist1h2ae" "Tyrobp" "Klra4" "Xcl1"# Plot variable features, label the top 10, save the plot
vfp1 <- VariableFeaturePlot(allsamplesgood)
vfp1 <- LabelPoints(plot = vfp1, points = top10, repel = TRUE)
vfp1
ggsave(paste(outputpath,"LCMV-variablefeatureplot.pdf",sep=""), plot = vfp1, width = 20, height = 15, units = "cm")
# Analyze variable features (top 30 default genes)
variable_feature_analysis <- SCassist_analyze_variable_features("allsamplesgood", experimental_design = experimental_design, api_key_file = api_key_file, llm_server = "google")## ## Enriched Gene Ontologies and Pathways:
##
## The provided list of genes suggests enrichment in several key pathways and functions related to **immune response, cell cycle, and inflammation**.
##
## **1. Immune Response:**
##
## * **NK cell activation and cytotoxicity:** Genes like *Gzma, Gzmb, Nkg7, Klra1, Klra4, Klra7, Klra8, Klra9* are all involved in NK cell function. *Gzma* and *Gzmb* are granzyme proteins that induce apoptosis in target cells, while *Nkg7* is a receptor involved in NK cell activation. The *Klra* genes encode killer cell lectin-like receptors, which play a role in NK cell recognition and activation.
## * **Chemokine signaling:** Genes like *Ccl3, Ccl4, Ccl5, Cxcl10, Xcl1* are chemokines involved in attracting and activating immune cells, particularly T cells and NK cells.
## * **Antigen presentation:** *Tyrobp* encodes a protein involved in the signaling pathway of the Fc receptor, which plays a role in antigen presentation and immune cell activation.
## * **T cell activation and differentiation:** *Ikzf2* is a transcription factor involved in T cell development and differentiation.
##
## **2. Cell Cycle and Proliferation:**
##
## * **Cell cycle regulation:** *Mki67* is a marker of proliferating cells, while *Top2a* is involved in DNA replication and repair.
## * **Cell growth and differentiation:** *Stmn1* is a protein involved in microtubule dynamics and cell growth.
##
## **3. Inflammation and Tissue Repair:**
##
## * **Inflammation:** *S100a4, S100a6, Hmgb2* are involved in inflammation and immune response.
## * **Extracellular matrix remodeling:** *Spp1* encodes osteopontin, a protein involved in extracellular matrix remodeling and inflammation.
##
## **Relevance to Experimental Design:**
##
## The identified gene ontologies and pathways are highly relevant to the experimental design involving NK, CD4+, and CD8+ T cells isolated from lymphocytic choriomeningitis virus (LCMV) infected WT and Ifng - CTCF binding site mutant mice.
##
## * **LCMV infection:** LCMV is a virus that induces a robust immune response, involving both NK cells and T cells. The genes identified are likely involved in the immune response to LCMV infection.
## * **NK cell function:** The enrichment of genes involved in NK cell activation and cytotoxicity suggests that NK cells play a significant role in the immune response to LCMV infection.
## * **T cell function:** The presence of genes involved in T cell activation, differentiation, and chemokine signaling indicates that T cells are also actively involved in the immune response.
## * **Ifng - CTCF binding site mutant mice:** The mutation in the Ifng - CTCF binding site likely affects the expression of genes involved in the immune response, potentially altering the activation and function of NK and T cells.
##
## By analyzing the expression of these genes in WT and mutant mice, researchers can gain insights into the role of NK and T cells in the immune response to LCMV infection and how the mutation in the Ifng - CTCF binding site affects these processes.4. PCA
Here we perform PCA and ask SCassist_recommend_pcs to analzye the PCs to provide recommendation for the number of PCs to use for downstream analysis.
# Run PCA to generate 50 pcs
allsamplesgood <- RunPCA(object = allsamplesgood, npcs=50)
# Print genes in the top 5 pcs
print(allsamplesgood[["pca"]], dims = 1:5, nfeatures = 5)## PC_ 1
## Positive: Gzma, Ccl5, Gzmb, Tyrobp, Fcer1g
## Negative: Ly6c1, Ltb, Ikzf2, Rps27, Cd3d
## PC_ 2
## Positive: Hist1h2ap, Hist1h1b, Hist1h2ae, Top2a, Mki67
## Negative: Ccl5, Gzma, Gzmb, Tyrobp, Ccl4
## PC_ 3
## Positive: Ccl5, Gzma, Cd8b1, Rps27, Rps24
## Negative: C1qa, Slc40a1, C1qc, Hmox1, Vcam1
## PC_ 4
## Positive: Ly6c2, Cd8b1, Ly6c1, Dapl1, Plac8
## Negative: Ikzf2, Tnfrsf4, Foxp3, Tnfrsf9, Izumo1r
## PC_ 5
## Positive: Ccl4, Xcl1, Ccl3, Gzmb, Spp1
## Negative: Ccl5, Gzma, Cma1, S100a6, Cenpf4.1. Choosing appropriate number of PCs
Here we ask SCassist_recommend_pcs to analzye the PCs to provide recommendation for the number of PCs to use for downstream analysis. We also ask SCassist_analyze_pcs to analzye the SCassist recommended PCs and provide insights on the gene sets identified.
# Identify appropriate PCs to use
pc_recommendation=SCassist_recommend_pcs("allsamplesgood", experimental_design = experimental_design, api_key_file = api_key_file, llm_server = "google")## Based on the provided variance explained by each PC, the optimal number of PCs to use for downstream analysis is **10**.
##
## Here's the reasoning:
##
## * **Elbow point:** While there isn't a super clear "elbow" in the scree plot, there's a noticeable change in slope around PC10. The variance explained drops significantly after PC10, indicating that the subsequent PCs capture less substantial information.
## * **Variance explained:** The first 10 PCs capture approximately 70% of the total variance (34.35% + 15.63% + 6.67% + 3.91% + 2.88% + 2.5% + 2.35% + 2.11% + 2.01% + 1.9% = 70.31%). This is a good balance between capturing most of the variation and avoiding excessive complexity.
## * **Balance between complexity and interpretability:** Using 10 PCs provides a reasonable level of detail while maintaining interpretability. Including more PCs might capture subtle variations but could lead to overfitting and make the analysis harder to understand.
##
## Therefore, using 10 PCs strikes a balance between capturing most of the variance and maintaining a manageable level of complexity for downstream analysis.# Identify number of PCs that explains majority of variations
ElbowPlot(allsamplesgood, ndims=50, reduction = "pca")
# Visualize the genes in the first PC
VizDimLoadings(allsamplesgood, dims = 1, ncol = 1) + theme_minimal(base_size = 8)
# Visualize the genes in the second PC
VizDimLoadings(allsamplesgood, dims = 2, ncol = 1) + theme_minimal(base_size = 8)
# Set the identity column to WT vs KO
Idents(allsamplesgood) <- "orig.ident"
# Visualize the cells after pca
DimPlot(object = allsamplesgood, reduction = "pca", group.by = "orig.ident")
# Plot heatmap with cells=500 plotting cells with extreme cells on both ends of spectrum
DimHeatmap(object = allsamplesgood, dims = 1:5, cells = 500, balanced = TRUE)
# Here we are also analyzing the gene sets that make up the top 5 PCs, top_n_pc_contributing_genes = 50
pc_analyzed=SCassist_analyze_pcs("allsamplesgood", num_pcs = 5, experimental_design = experimental_design, api_key_file = api_key_file, llm_server = "google")##
## **PC1 Summary:**
## The top contributing genes for PC1 are primarily associated with cytotoxic T cell function and inflammation. Genes like *Gzma*, *Gzmb*, *Prf1*, *Nkg7*, and *Klra* family members are all involved in granule exocytosis and killing of target cells. *Ccl3*, *Ccl4*, and *Ccl5* are chemokines that attract immune cells to sites of inflammation. *Tyrobp*, *Fcer1g*, and *Irf8* are involved in immune signaling and activation. The presence of multiple histone genes (e.g., *Hist1h2ap*, *Hist1h1b*, *Hist1h2ae*, *Hist1h3c*) suggests a potential role for changes in gene expression regulation.
##
## Taken together, these genes suggest that PC1 might be capturing variations in the cytotoxic activity and inflammatory response of CD4+ and CD8+ T cells in the context of lymphocytic choriomeningitis virus infection. The differences in gene expression between WT and Ifng-CTCF binding site mutant mice could reflect variations in the ability of these cells to control viral infection, potentially due to altered cytokine production, cytotoxic activity, or immune signaling.
##
##
## **PC2 Summary:**
## The top contributing genes for PC2 are primarily involved in cell cycle regulation, proliferation, and immune response. Genes like *Mki67*, *Ccna2*, *Cdk1*, *Cenpf*, *Kif11*, and *Tuba1b* are known to be crucial for cell division and proliferation. *Gzma*, *Gzmb*, *Klra8*, *Ccl5*, and *Ccl4* are associated with cytotoxic T cell activity and cytokine production, indicating an immune response. The presence of histone genes like *Hist1h2ap*, *Hist1h1b*, *Hist1h2ae*, and *Hist1h3c* suggests active transcription and chromatin remodeling.
##
## Therefore, PC2 likely captures variations in the cellular state related to T cell activation, proliferation, and cytotoxic activity. This could be driven by differences in the immune response to lymphocytic choriomeningitis virus infection between WT and Ifng - CTCF binding site mutant mice.
##
##
## **PC3 Summary:**
## The top contributing genes for PC3 are primarily involved in immune response, inflammation, and macrophage activation. Genes like C1qa, C1qc, C1qb, and Cfb are components of the complement system, which plays a crucial role in innate immunity and inflammation. Other genes like Hmox1, Vcam1, and Cd5l are involved in macrophage activation and recruitment. Additionally, genes like Ifitm3, Ifi207, and Ccr3 are associated with antiviral responses. These genes suggest that PC3 might be capturing variations in the immune response to lymphocytic choriomeningitis virus (LCMV) infection, particularly the activation and recruitment of macrophages and the induction of antiviral responses. The presence of genes like Apoe and Lpl, involved in lipid metabolism, might also indicate a potential role of lipid metabolism in the immune response to LCMV infection.
##
##
## **PC4 Summary:**
## The top contributing genes for PC4 include several key players in T cell activation, differentiation, and function. Notably, we see genes like **Ikzf2 (Helios)**, a transcription factor associated with regulatory T cell (Treg) development, **Foxp3**, the master regulator of Treg function, and **Ctla4**, a negative regulator of T cell activation. Additionally, genes like **Cd8a**, **Cd8b1**, and **Gzmb** point towards cytotoxic T cell (CTL) activity. The presence of **Tnfrsf4 (OX40)** and **Tnfrsf9 (4-1BB)**, both costimulatory molecules involved in T cell activation and survival, further supports this notion.
##
## Taken together, these genes suggest that PC4 might be capturing variations in the balance between Treg and CTL populations, potentially reflecting differences in immune responses to lymphocytic choriomeningitis virus (LCMV) infection between WT and Ifng-CTCF binding site mutant mice. This could be driven by altered T cell differentiation pathways, changes in T cell activation and suppression, or variations in the expression of effector molecules like granzyme B.
##
##
## **PC5 Summary:**
## The top contributing genes for PC5 are primarily associated with immune response and inflammation. Genes like *Ccl4, Xcl1, Ccl3, Ccl5, Cxcl10, Ifng, Gzmb, Gzma, Tnfrsf9* are chemokines and cytokines involved in attracting and activating immune cells, particularly T cells and NK cells. *Egr1, Egr3, Nr4a1, Nr4a3, Nr4a2* are transcription factors that regulate the expression of genes involved in immune responses. *Irf8, Isg15, Ifi204, Ifitm3* are interferon-stimulated genes involved in antiviral responses. *Cd69, Cd7, Cd160* are cell surface receptors involved in T cell activation and differentiation. *Prf1, Vim, Lgals1, S100a4, S100a6, S100a10* are involved in cell adhesion, migration, and cytotoxicity. These genes suggest that PC5 might be capturing variations in the activation, differentiation, and effector functions of T cells and NK cells during lymphocytic choriomeningitis virus infection. The presence of genes involved in both innate and adaptive immune responses suggests that PC5 might be reflecting the interplay between these two arms of the immune system in controlling viral infection.
##
##
## ## Overall Summary of Principal Components (PCs)
##
## The five principal components (PCs) identified in this study capture distinct aspects of the immune response to lymphocytic choriomeningitis virus (LCMV) infection, particularly in the context of differences between wild-type (WT) and Ifng-CTCF binding site mutant mice.
##
## **PC1:** Primarily reflects variations in **cytotoxic T cell function and inflammation**. This component highlights genes involved in granule exocytosis, target cell killing, chemokine production, and immune signaling. Differences in PC1 scores between WT and mutant mice could indicate variations in their ability to control viral infection through altered cytokine production, cytotoxic activity, or immune signaling.
##
## **PC2:** Captures variations in **T cell activation, proliferation, and cytotoxic activity**. This component includes genes crucial for cell division, proliferation, and immune response, suggesting differences in the cellular state of T cells between WT and mutant mice.
##
## **PC3:** Focuses on **immune response, inflammation, and macrophage activation**. This component highlights genes involved in the complement system, macrophage activation and recruitment, and antiviral responses. Differences in PC3 scores could reflect variations in the activation and recruitment of macrophages and the induction of antiviral responses during LCMV infection.
##
## **PC4:** Reflects variations in the **balance between regulatory T cells (Tregs) and cytotoxic T cells (CTLs)**. This component includes genes associated with Treg development and function, as well as CTL activity. Differences in PC4 scores could indicate altered T cell differentiation pathways, changes in T cell activation and suppression, or variations in the expression of effector molecules like granzyme B.
##
## **PC5:** Captures variations in the **activation, differentiation, and effector functions of T cells and NK cells**. This component includes genes involved in chemokine and cytokine production, T cell activation and differentiation, and antiviral responses. The presence of genes involved in both innate and adaptive immune responses suggests that PC5 might be reflecting the interplay between these two arms of the immune system in controlling viral infection.
##
## **Overall, these PCs provide a comprehensive view of the complex immune response to LCMV infection, highlighting the distinct roles of different immune cell populations and pathways in controlling viral infection. The differences observed between WT and Ifng-CTCF binding site mutant mice suggest that the CTCF binding site in the Ifng gene plays a crucial role in regulating these immune responses.**5. FindNeighbors
Here we ask SCassist_recommend_k to identify optimum number for k, to use with FindNeighbors function, given the number of pcs we are going to use for the downstream analysis.
# Run SCassist_recommend_k to identify data based K
recommended_k <- SCassist_recommend_k("allsamplesgood", num_pcs = 10, api_key_file = api_key_file, llm_server = "google", experimental_design = experimental_design)## ## Recommended K: 15-30
##
## ## Reasoning:
##
## The `k.param` value in `FindNeighbors()` determines the number of nearest neighbors considered for each cell when constructing the k-nearest neighbor graph. A higher `k.param` value leads to a more connected graph, potentially capturing broader relationships between cells. However, too high a value can blur distinct cell populations.
##
## Given your dataset of approximately 8555 cells and the use of 10 principal components, a range of `k.param` values between 15 and 30 is recommended. This range allows for a balance between capturing local relationships within clusters while still considering broader connections between different cell types. Starting with a lower value like 15 and gradually increasing to 30 can help identify the optimal `k.param` for your specific data and clustering goals.# Run findneighbors
allsamplesgood <- FindNeighbors(allsamplesgood, dims = 1:10, k.param = 15, return.neighbor = TRUE)6. Clustering
Here we ask SCassist_recommend_res to analyze the number of cells, number of genes, mean expression variability and median neighbor distance, and recommend appropriate cluster resolutions.
# Run SCassist_recommend_res
recommended_res <- SCassist_recommend_res("allsamplesgood", api_key_file = api_key_file, llm_server = "google")## Based on the data characteristics, I recommend the following resolution range:
##
## **Recommended Resolution:** seq(0.4, 1.2, 0.1)
##
## **Reasoning:**
##
## The mean expression variability of 8.66843315573589 suggests a moderate level of heterogeneity in your dataset. This indicates that there might be subtle differences in gene expression between cell populations. The median neighbor distance of 4.83575773239136 in the k-nearest neighbor graph further supports this notion, as it implies a relatively dense and interconnected cell landscape.
##
## Therefore, a resolution range of seq(0.4, 1.2, 0.1) is recommended. This range allows for the identification of both distinct and subtle cell populations. Lower resolutions (0.4-0.7) will capture broader, more general cell types, while higher resolutions (0.8-1.2) will enable the identification of finer, more specific subpopulations. By exploring this range, you can effectively delineate the cellular landscape of your dataset and uncover potentially interesting biological insights.allsamplesgood <- FindNeighbors(allsamplesgood, dims = 1:10, k.param = 15)
# Perform clustering using the above identified resolution.
allsamplesgood <- FindClusters(allsamplesgood, resolution = seq(0.4,1.2,0.1))## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8555
## Number of edges: 241353
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8999
## Number of communities: 10
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8555
## Number of edges: 241353
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8877
## Number of communities: 12
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8555
## Number of edges: 241353
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8781
## Number of communities: 13
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8555
## Number of edges: 241353
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8688
## Number of communities: 14
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8555
## Number of edges: 241353
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8598
## Number of communities: 14
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8555
## Number of edges: 241353
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8519
## Number of communities: 17
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8555
## Number of edges: 241353
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8441
## Number of communities: 21
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8555
## Number of edges: 241353
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8370
## Number of communities: 21
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 8555
## Number of edges: 241353
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8317
## Number of communities: 23
## Elapsed time: 0 secondshead(allsamplesgood)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCAAGAGTTGTA-WT WT 17880 4120 3.333333
## AAACCCAAGATTTGCC-WT WT 7916 2395 2.627590
## AAACCCACAAGTATCC-WT WT 9059 2865 1.953858
## AAACCCAGTCCAGGTC-WT WT 10948 2939 2.831567
## AAACCCATCCTGGGAC-WT WT 16768 3484 0.942271
## AAACGAACAGAGAGGG-WT WT 10467 3159 3.668673
## AAACGAACAGCGACAA-WT WT 12996 2853 1.569714
## AAACGAACAGTGTGGA-WT WT 26674 4654 1.570818
## AAACGAATCTAGCCTC-WT WT 11434 3100 1.687948
## AAACGCTAGCGTTCCG-WT WT 28638 4519 3.579161
## percent.ribo percent.hb nCount_SCT nFeature_SCT
## AAACCCAAGAGTTGTA-WT 22.80761 0.011185682 10914 3831
## AAACCCAAGATTTGCC-WT 16.28348 0.000000000 9352 2394
## AAACCCACAAGTATCC-WT 16.81201 0.000000000 9513 2863
## AAACCCAGTCCAGGTC-WT 25.30142 0.009134088 10388 2936
## AAACCCATCCTGGGAC-WT 35.07276 0.005963740 10691 3307
## AAACGAACAGAGAGGG-WT 22.92921 0.009553836 10239 3158
## AAACGAACAGCGACAA-WT 42.22068 0.000000000 10582 2853
## AAACGAACAGTGTGGA-WT 26.22029 0.000000000 10379 2877
## AAACGAATCTAGCCTC-WT 28.44149 0.000000000 10519 3100
## AAACGCTAGCGTTCCG-WT 43.16642 0.000000000 9760 2387
## SCT_snn_res.0.4 SCT_snn_res.0.5 SCT_snn_res.0.6
## AAACCCAAGAGTTGTA-WT 5 4 2
## AAACCCAAGATTTGCC-WT 4 6 5
## AAACCCACAAGTATCC-WT 5 4 2
## AAACCCAGTCCAGGTC-WT 3 5 1
## AAACCCATCCTGGGAC-WT 0 2 7
## AAACGAACAGAGAGGG-WT 7 7 8
## AAACGAACAGCGACAA-WT 2 1 0
## AAACGAACAGTGTGGA-WT 8 8 9
## AAACGAATCTAGCCTC-WT 3 5 1
## AAACGCTAGCGTTCCG-WT 0 2 6
## SCT_snn_res.0.7 SCT_snn_res.0.8 SCT_snn_res.0.9
## AAACCCAAGAGTTGTA-WT 5 3 8
## AAACCCAAGATTTGCC-WT 4 5 4
## AAACCCACAAGTATCC-WT 13 13 14
## AAACCCAGTCCAGGTC-WT 1 1 3
## AAACCCATCCTGGGAC-WT 6 7 6
## AAACGAACAGAGAGGG-WT 7 8 7
## AAACGAACAGCGACAA-WT 2 0 2
## AAACGAACAGTGTGGA-WT 9 9 10
## AAACGAATCTAGCCTC-WT 1 1 3
## AAACGCTAGCGTTCCG-WT 6 6 5
## SCT_snn_res.1 SCT_snn_res.1.1 SCT_snn_res.1.2
## AAACCCAAGAGTTGTA-WT 6 6 9
## AAACCCAAGATTTGCC-WT 5 9 6
## AAACCCACAAGTATCC-WT 17 18 18
## AAACCCAGTCCAGGTC-WT 12 13 13
## AAACCCATCCTGGGAC-WT 4 5 3
## AAACGAACAGAGAGGG-WT 8 7 8
## AAACGAACAGCGACAA-WT 1 1 19
## AAACGAACAGTGTGGA-WT 9 8 11
## AAACGAATCTAGCCTC-WT 3 3 4
## AAACGCTAGCGTTCCG-WT 7 4 5
## seurat_clusters
## AAACCCAAGAGTTGTA-WT 9
## AAACCCAAGATTTGCC-WT 6
## AAACCCACAAGTATCC-WT 18
## AAACCCAGTCCAGGTC-WT 13
## AAACCCATCCTGGGAC-WT 3
## AAACGAACAGAGAGGG-WT 8
## AAACGAACAGCGACAA-WT 19
## AAACGAACAGTGTGGA-WT 11
## AAACGAATCTAGCCTC-WT 4
## AAACGCTAGCGTTCCG-WT 5# Count number of clusters at each resolution
sapply(grep("res",colnames(allsamplesgood@meta.data),value = TRUE),
function(x) length(unique(allsamplesgood@meta.data[,x])))## SCT_snn_res.0.4 SCT_snn_res.0.5 SCT_snn_res.0.6 SCT_snn_res.0.7 SCT_snn_res.0.8
## 10 12 13 14 14
## SCT_snn_res.0.9 SCT_snn_res.1 SCT_snn_res.1.1 SCT_snn_res.1.2
## 17 21 21 23# change default identity
Idents(allsamplesgood) <- "SCT_snn_res.1.1"
# list cell number in each cluster for HC vs UV
table(Idents(allsamplesgood),allsamplesgood$orig.ident)##
## KO WT
## 0 505 400
## 1 438 417
## 2 531 280
## 3 128 618
## 4 478 209
## 5 277 399
## 6 232 375
## 7 303 154
## 8 135 317
## 9 197 162
## 10 111 231
## 11 149 183
## 12 95 187
## 13 17 209
## 14 69 95
## 15 110 51
## 16 58 101
## 17 70 56
## 18 36 62
## 19 30 53
## 20 10 17set.seed(20)
# Run UMAP with identified dimensions
allsamplesgood <- RunUMAP(allsamplesgood, dims = 1:10)
# Plot umap
DimPlot(object = allsamplesgood, pt.size=0.5, reduction = "umap", label = T)
# Color by WT vs KO
DimPlot(object = allsamplesgood, pt.size=0.5, reduction = "umap", group.by = "orig.ident")