Tutorial: Single Cell RNA-Seq Analysis of PBMC's from a healthy human Vs an Uveitis patient, using SCassist - an AI-Based Workflow Assistant
1. Quality Control
This document takes the cellranger count data for two samples (a healthy control and a Uveitis patient) 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)
library(dplyr)1.1. Data setup
# Set the input, output paths
setwd("/singlecellassistant")
datapath="exampleData"
outputpath="/singlecellassistant/exampleOutput_G/BCRUV/"
sampleinfopath="/singlecellassistant/exampleData/BCRUV/BCRUV_SampleInformation.tab"
# Load RData for example dataset
load("/singlecellassistant/exampleOutput_G/BCRUV/bcruv_G_new.RData")
# Read sample information file
sampleinformation=read.csv(sampleinfopath,sep="\t")
head(sampleinformation)## RunID file_name
## 1 HVHL7DRXX_19022128 HVHL7DRXX_19022128_S3_L001_R1_001.fastq.gz
## 2 HVHL7DRXX_19022128 HVHL7DRXX_19022128_S3_L001_R2_001.fastq.gz
## 3 HVHL7DRXX_19022126 HVHL7DRXX_19022126_S2_L001_R1_001.fastq.gz
## 4 HVHL7DRXX_19022126 HVHL7DRXX_19022126_S2_L001_R2_001.fastq.gz
## 5 HVHL7DRXX_19022124 HVHL7DRXX_19022124_S4_L001_R1_001.fastq.gz
## 6 HVHL7DRXX_19022124 HVHL7DRXX_19022124_S4_L001_R2_001.fastq.gz
## folder_location_in_biowulf submitted_name
## 1 /data/../TotalSeq/RawData/201230_INGENS_0255_B_HVHL7DRXX/ NS3R189BTS
## 2 /data/../TotalSeq/RawData/201230_INGENS_0255_B_HVHL7DRXX/ NS3R189BTS
## 3 /data/../TotalSeq/RawData/201230_INGENS_0255_B_HVHL7DRXX/ NS3R189BTS
## 4 /data/../TotalSeq/RawData/201230_INGENS_0255_B_HVHL7DRXX/ NS3R189BTS
## 5 /data/../TotalSeq/RawData/201230_INGENS_0255_B_HVHL7DRXX/ NS3R189BTS
## 6 /data/../TotalSeq/RawData/201230_INGENS_0255_B_HVHL7DRXX/ NS3R189BTS
## project machine_type pcnt_optical_duplicates clusters read_length
## 1 Sen Uveitis NovaSeq6000_SP NA 87268107 28,91
## 2 Sen Uveitis NovaSeq6000_SP NA 87268107 28,91
## 3 Sen Uveitis NovaSeq6000_SP NA 130870111 28,91
## 4 Sen Uveitis NovaSeq6000_SP NA 130870111 28,91
## 5 Sen Uveitis NovaSeq6000_SP NA 113258165 28,91
## 6 Sen Uveitis NovaSeq6000_SP NA 113258165 28,91
## experiment_type run_date RunType SampleType ShortSubmittedName
## 1 SCRNA-Seq 12/30/20 0:00 Gene Expression HC 3R189
## 2 SCRNA-Seq 12/30/20 0:00 Gene Expression HC 3R189
## 3 SCRNA-Seq 12/30/20 0:00 Gene Expression HC 3R189
## 4 SCRNA-Seq 12/30/20 0:00 Gene Expression HC 3R189
## 5 SCRNA-Seq 12/30/20 0:00 Gene Expression HC 3R189
## 6 SCRNA-Seq 12/30/20 0:00 Gene Expression HC 3R189
## ShortNameType SampleDescription Age Sex Race SampleBatch Sample.prepared
## 1 3R189_HC Healthy Control 37 FALSE C 1 12/2/20
## 2 3R189_HC Healthy Control 37 FALSE C 1 12/2/20
## 3 3R189_HC Healthy Control 37 FALSE C 1 12/2/20
## 4 3R189_HC Healthy Control 37 FALSE C 1 12/2/20
## 5 3R189_HC Healthy Control 37 FALSE C 1 12/2/20
## 6 3R189_HC Healthy Control 37 FALSE C 1 12/2/20
## DateReceived TreatmentCategory CMV_IgG CMV_IgM EBVCA_IgG
## 1 1/4/21 Healthy negative negative positive
## 2 1/4/21 Healthy negative negative positive
## 3 1/4/21 Healthy negative negative positive
## 4 1/4/21 Healthy negative negative positive
## 5 1/4/21 Healthy negative negative positive
## 6 1/4/21 Healthy negative negative positive# Specify samples for analysis.
samples=c("NS3R189BTS","NS7R65BBTS")
# Read the data file
NS3R189BTS <- Read10X_h5(file.path("/BCRUV/BCRUV_NS3R189BTS_filtered_feature_bc_matrix.h5"), use.names = T)
NS7R65BBTS <- Read10X_h5(file.path("/BCRUV/BCRUV_NS7R65BBTS_filtered_feature_bc_matrix.h5"), use.names = T)
# Add sample names to columnnames
colnames(NS3R189BTS$`Gene Expression`) <- paste(sapply(strsplit(colnames(NS3R189BTS$`Gene Expression`),split="-"),'[[',1L),"NS3R189BTS",sep="-")
colnames(NS7R65BBTS$`Gene Expression`) <- paste(sapply(strsplit(colnames(NS7R65BBTS$`Gene Expression`),split="-"),'[[',1L),"NS7R65BBTS",sep="-")
# Create a Seurat data object from the gex matrix
NS3R189BTS <- CreateSeuratObject(counts = NS3R189BTS[["Gene Expression"]], names.field = 2,names.delim = "\\-")
NS7R65BBTS <- CreateSeuratObject(counts = NS7R65BBTS[["Gene Expression"]], names.field = 2,names.delim = "\\-")
# Merge seurat objects
allsamples <- merge(NS3R189BTS, NS7R65BBTS, project = "UV")
allsamples## An object of class Seurat
## 36601 features across 39090 samples within 1 assay
## Active assay: RNA (36601 features, 0 variable features)
## 2 layers present: counts.1, counts.2head(allsamples)## orig.ident nCount_RNA nFeature_RNA
## AAACCCACAAGTCATC-NS3R189BTS NS3R189BTS 2431 1431
## AAACCCACACTGGAAG-NS3R189BTS NS3R189BTS 2676 1078
## AAACCCACAGCACACC-NS3R189BTS NS3R189BTS 3394 1760
## AAACCCACAGGAGGAG-NS3R189BTS NS3R189BTS 2160 1279
## AAACCCACATGACGTT-NS3R189BTS NS3R189BTS 2329 1271
## AAACCCAGTAGTTACC-NS3R189BTS NS3R189BTS 5362 2261
## AAACCCAGTCCAAATC-NS3R189BTS NS3R189BTS 7926 2660
## AAACCCAGTCCATACA-NS3R189BTS NS3R189BTS 2239 1283
## AAACCCAGTGAGTTTC-NS3R189BTS NS3R189BTS 2339 1291
## AAACCCAGTTCTTAGG-NS3R189BTS NS3R189BTS 3537 1608table(Idents(allsamples))##
## NS3R189BTS NS7R65BBTS
## 24790 14300# 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
allsamples$percent.mito <- PercentageFeatureSet(allsamples, pattern = "^MT-")
summary(allsamples$percent.mito)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.300 3.517 4.084 5.197 86.180head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCACAAGTCATC-NS3R189BTS NS3R189BTS 2431 1431 4.1135335
## AAACCCACACTGGAAG-NS3R189BTS NS3R189BTS 2676 1078 0.8221226
## AAACCCACAGCACACC-NS3R189BTS NS3R189BTS 3394 1760 5.1856217
## AAACCCACAGGAGGAG-NS3R189BTS NS3R189BTS 2160 1279 6.5277778
## AAACCCACATGACGTT-NS3R189BTS NS3R189BTS 2329 1271 4.3795620
## AAACCCAGTAGTTACC-NS3R189BTS NS3R189BTS 5362 2261 3.8418501
## AAACCCAGTCCAAATC-NS3R189BTS NS3R189BTS 7926 2660 4.6177139
## AAACCCAGTCCATACA-NS3R189BTS NS3R189BTS 2239 1283 5.0468959
## AAACCCAGTGAGTTTC-NS3R189BTS NS3R189BTS 2339 1291 5.3014109
## AAACCCAGTTCTTAGG-NS3R189BTS NS3R189BTS 3537 1608 2.2900763# 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.
## 0.5956 20.4335 25.8115 26.1630 32.7232 68.7378head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCACAAGTCATC-NS3R189BTS NS3R189BTS 2431 1431 4.1135335
## AAACCCACACTGGAAG-NS3R189BTS NS3R189BTS 2676 1078 0.8221226
## AAACCCACAGCACACC-NS3R189BTS NS3R189BTS 3394 1760 5.1856217
## AAACCCACAGGAGGAG-NS3R189BTS NS3R189BTS 2160 1279 6.5277778
## AAACCCACATGACGTT-NS3R189BTS NS3R189BTS 2329 1271 4.3795620
## AAACCCAGTAGTTACC-NS3R189BTS NS3R189BTS 5362 2261 3.8418501
## AAACCCAGTCCAAATC-NS3R189BTS NS3R189BTS 7926 2660 4.6177139
## AAACCCAGTCCATACA-NS3R189BTS NS3R189BTS 2239 1283 5.0468959
## AAACCCAGTGAGTTTC-NS3R189BTS NS3R189BTS 2339 1291 5.3014109
## AAACCCAGTTCTTAGG-NS3R189BTS NS3R189BTS 3537 1608 2.2900763
## percent.ribo
## AAACCCACAAGTCATC-NS3R189BTS 18.510901
## AAACCCACACTGGAAG-NS3R189BTS 42.526158
## AAACCCACAGCACACC-NS3R189BTS 19.151444
## AAACCCACAGGAGGAG-NS3R189BTS 8.842593
## AAACCCACATGACGTT-NS3R189BTS 18.978102
## AAACCCAGTAGTTACC-NS3R189BTS 18.985453
## AAACCCAGTCCAAATC-NS3R189BTS 30.052990
## AAACCCAGTCCATACA-NS3R189BTS 19.830281
## AAACCCAGTGAGTTTC-NS3R189BTS 21.675930
## AAACCCAGTTCTTAGG-NS3R189BTS 23.353124# 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.00000 0.00000 0.00000 0.07317 0.00000 96.37250head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCACAAGTCATC-NS3R189BTS NS3R189BTS 2431 1431 4.1135335
## AAACCCACACTGGAAG-NS3R189BTS NS3R189BTS 2676 1078 0.8221226
## AAACCCACAGCACACC-NS3R189BTS NS3R189BTS 3394 1760 5.1856217
## AAACCCACAGGAGGAG-NS3R189BTS NS3R189BTS 2160 1279 6.5277778
## AAACCCACATGACGTT-NS3R189BTS NS3R189BTS 2329 1271 4.3795620
## AAACCCAGTAGTTACC-NS3R189BTS NS3R189BTS 5362 2261 3.8418501
## AAACCCAGTCCAAATC-NS3R189BTS NS3R189BTS 7926 2660 4.6177139
## AAACCCAGTCCATACA-NS3R189BTS NS3R189BTS 2239 1283 5.0468959
## AAACCCAGTGAGTTTC-NS3R189BTS NS3R189BTS 2339 1291 5.3014109
## AAACCCAGTTCTTAGG-NS3R189BTS NS3R189BTS 3537 1608 2.2900763
## percent.ribo percent.hb
## AAACCCACAAGTCATC-NS3R189BTS 18.510901 0.00000000
## AAACCCACACTGGAAG-NS3R189BTS 42.526158 0.00000000
## AAACCCACAGCACACC-NS3R189BTS 19.151444 0.00000000
## AAACCCACAGGAGGAG-NS3R189BTS 8.842593 0.00000000
## AAACCCACATGACGTT-NS3R189BTS 18.978102 0.04293688
## AAACCCAGTAGTTACC-NS3R189BTS 18.985453 0.01864976
## AAACCCAGTCCAAATC-NS3R189BTS 30.052990 0.02523341
## AAACCCAGTCCATACA-NS3R189BTS 19.830281 0.04466280
## AAACCCAGTGAGTTTC-NS3R189BTS 21.675930 0.00000000
## AAACCCAGTTCTTAGG-NS3R189BTS 23.353124 0.00000000# Check the total number of cells in this experiment
table(allsamples$orig.ident)##
## NS3R189BTS NS7R65BBTS
## 24790 14300api_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:** 1000
## * **Upper Cutoff:** 15000
##
## **Rationale:** The 5th percentile suggests a lower limit of 1283, but considering the minimum value of 500, a lower cutoff of 1000 might be more appropriate to avoid removing potentially valid cells with low counts. The 95th percentile suggests an upper limit of 9600.55, but considering the maximum value of 80702, an upper cutoff of 15000 seems more reasonable to remove extreme outliers while retaining cells with high expression.
##
## **nFeature_RNA:**
##
## * **Lower Cutoff:** 800
## * **Upper Cutoff:** 4000
##
## **Rationale:** The 5th percentile suggests a lower limit of 780, but considering the minimum value of 46, a lower cutoff of 800 might be more appropriate to avoid removing potentially valid cells with low feature counts. The 95th percentile suggests an upper limit of 3069, but considering the maximum value of 8873, an upper cutoff of 4000 seems more reasonable to remove extreme outliers while retaining cells with high feature counts.
##
## **percent.mito:**
##
## * **Upper Cutoff:** 10
##
## **Rationale:** The 95th percentile suggests an upper limit of 8.68, but considering the maximum value of 86.18, an upper cutoff of 10 seems more reasonable to remove cells with extremely high mitochondrial content while retaining cells with moderate mitochondrial expression.
##
## **percent.ribo:**
##
## * **Upper Cutoff:** 50
##
## **Rationale:** The 95th percentile suggests an upper limit of 41.89, but considering the maximum value of 68.74, an upper cutoff of 50 seems more reasonable to remove cells with extremely high ribosomal content while retaining cells with moderate ribosomal expression.
##
## **percent.hb:**
##
## * **Upper Cutoff:** 0.1
##
## **Rationale:** The 95th percentile suggests an upper limit of 0.0499, but considering the maximum value of 96.37, an upper cutoff of 0.1 seems more reasonable to remove cells with extremely high hemoglobin content while retaining cells with moderate hemoglobin expression.
##
## **Important Note:** These are just recommendations based on the provided data. It is crucial to test a range of values around these recommendations to determine the optimal cutoffs for your specific dataset and analysis goals. Visualizing the distributions of each metric can help in making informed decisions about the cutoffs.# The SCassist recommendation is also stored in "allsamplesquality"
#cat(allsamplesquality)1.3. Quality Filtering
Here we use the SCassist recommended QC cutoff values obtained in the previous step to filter and visualize the before and after QC data.
# Filter data based on SCassists recommendation
allsamplesgood <- subset(allsamples, subset = nCount_RNA > 1000 & nCount_RNA < 15000 & nFeature_RNA > 800 & nFeature_RNA < 4000 & percent.mito < 10 & percent.ribo < 50 & percent.hb < 0.1)
# Counts before filtering
table(allsamples$orig.ident)##
## NS3R189BTS NS7R65BBTS
## 24790 14300# Counts after filtering
table(allsamplesgood$orig.ident)##
## NS3R189BTS NS7R65BBTS
## 22195 127721.4. Add metadata
# Assign sample types
allsamples$type <- plyr::mapvalues(
x = allsamples$orig.ident,
from = sampleinformation$submitted_name,
to = sampleinformation$SampleType
)
head(allsamples)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCACAAGTCATC-NS3R189BTS NS3R189BTS 2431 1431 4.1135335
## AAACCCACACTGGAAG-NS3R189BTS NS3R189BTS 2676 1078 0.8221226
## AAACCCACAGCACACC-NS3R189BTS NS3R189BTS 3394 1760 5.1856217
## AAACCCACAGGAGGAG-NS3R189BTS NS3R189BTS 2160 1279 6.5277778
## AAACCCACATGACGTT-NS3R189BTS NS3R189BTS 2329 1271 4.3795620
## AAACCCAGTAGTTACC-NS3R189BTS NS3R189BTS 5362 2261 3.8418501
## AAACCCAGTCCAAATC-NS3R189BTS NS3R189BTS 7926 2660 4.6177139
## AAACCCAGTCCATACA-NS3R189BTS NS3R189BTS 2239 1283 5.0468959
## AAACCCAGTGAGTTTC-NS3R189BTS NS3R189BTS 2339 1291 5.3014109
## AAACCCAGTTCTTAGG-NS3R189BTS NS3R189BTS 3537 1608 2.2900763
## percent.ribo percent.hb type
## AAACCCACAAGTCATC-NS3R189BTS 18.510901 0.00000000 HC
## AAACCCACACTGGAAG-NS3R189BTS 42.526158 0.00000000 HC
## AAACCCACAGCACACC-NS3R189BTS 19.151444 0.00000000 HC
## AAACCCACAGGAGGAG-NS3R189BTS 8.842593 0.00000000 HC
## AAACCCACATGACGTT-NS3R189BTS 18.978102 0.04293688 HC
## AAACCCAGTAGTTACC-NS3R189BTS 18.985453 0.01864976 HC
## AAACCCAGTCCAAATC-NS3R189BTS 30.052990 0.02523341 HC
## AAACCCAGTCCATACA-NS3R189BTS 19.830281 0.04466280 HC
## AAACCCAGTGAGTTTC-NS3R189BTS 21.675930 0.00000000 HC
## AAACCCAGTTCTTAGG-NS3R189BTS 23.353124 0.00000000 HCallsamplesgood$type <- plyr::mapvalues(
x = allsamplesgood$orig.ident,
from = sampleinformation$submitted_name,
to = sampleinformation$SampleType
)
head(allsamplesgood)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCACAAGTCATC-NS3R189BTS NS3R189BTS 2431 1431 4.1135335
## AAACCCACACTGGAAG-NS3R189BTS NS3R189BTS 2676 1078 0.8221226
## AAACCCACAGCACACC-NS3R189BTS NS3R189BTS 3394 1760 5.1856217
## AAACCCACAGGAGGAG-NS3R189BTS NS3R189BTS 2160 1279 6.5277778
## AAACCCACATGACGTT-NS3R189BTS NS3R189BTS 2329 1271 4.3795620
## AAACCCAGTAGTTACC-NS3R189BTS NS3R189BTS 5362 2261 3.8418501
## AAACCCAGTCCAAATC-NS3R189BTS NS3R189BTS 7926 2660 4.6177139
## AAACCCAGTCCATACA-NS3R189BTS NS3R189BTS 2239 1283 5.0468959
## AAACCCAGTGAGTTTC-NS3R189BTS NS3R189BTS 2339 1291 5.3014109
## AAACCCAGTTCTTAGG-NS3R189BTS NS3R189BTS 3537 1608 2.2900763
## percent.ribo percent.hb type
## AAACCCACAAGTCATC-NS3R189BTS 18.510901 0.00000000 HC
## AAACCCACACTGGAAG-NS3R189BTS 42.526158 0.00000000 HC
## AAACCCACAGCACACC-NS3R189BTS 19.151444 0.00000000 HC
## AAACCCACAGGAGGAG-NS3R189BTS 8.842593 0.00000000 HC
## AAACCCACATGACGTT-NS3R189BTS 18.978102 0.04293688 HC
## AAACCCAGTAGTTACC-NS3R189BTS 18.985453 0.01864976 HC
## AAACCCAGTCCAAATC-NS3R189BTS 30.052990 0.02523341 HC
## AAACCCAGTCCATACA-NS3R189BTS 19.830281 0.04466280 HC
## AAACCCAGTGAGTTTC-NS3R189BTS 21.675930 0.00000000 HC
## AAACCCAGTTCTTAGG-NS3R189BTS 23.353124 0.00000000 HC# View type counts
table(allsamples$type)##
## BCR-UV HC
## 14300 24790table(allsamplesgood$type)##
## BCR-UV HC
## 12772 22195# View data
head(allsamplesgood)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCACAAGTCATC-NS3R189BTS NS3R189BTS 2431 1431 4.1135335
## AAACCCACACTGGAAG-NS3R189BTS NS3R189BTS 2676 1078 0.8221226
## AAACCCACAGCACACC-NS3R189BTS NS3R189BTS 3394 1760 5.1856217
## AAACCCACAGGAGGAG-NS3R189BTS NS3R189BTS 2160 1279 6.5277778
## AAACCCACATGACGTT-NS3R189BTS NS3R189BTS 2329 1271 4.3795620
## AAACCCAGTAGTTACC-NS3R189BTS NS3R189BTS 5362 2261 3.8418501
## AAACCCAGTCCAAATC-NS3R189BTS NS3R189BTS 7926 2660 4.6177139
## AAACCCAGTCCATACA-NS3R189BTS NS3R189BTS 2239 1283 5.0468959
## AAACCCAGTGAGTTTC-NS3R189BTS NS3R189BTS 2339 1291 5.3014109
## AAACCCAGTTCTTAGG-NS3R189BTS NS3R189BTS 3537 1608 2.2900763
## percent.ribo percent.hb type
## AAACCCACAAGTCATC-NS3R189BTS 18.510901 0.00000000 HC
## AAACCCACACTGGAAG-NS3R189BTS 42.526158 0.00000000 HC
## AAACCCACAGCACACC-NS3R189BTS 19.151444 0.00000000 HC
## AAACCCACAGGAGGAG-NS3R189BTS 8.842593 0.00000000 HC
## AAACCCACATGACGTT-NS3R189BTS 18.978102 0.04293688 HC
## AAACCCAGTAGTTACC-NS3R189BTS 18.985453 0.01864976 HC
## AAACCCAGTCCAAATC-NS3R189BTS 30.052990 0.02523341 HC
## AAACCCAGTCCATACA-NS3R189BTS 19.830281 0.04466280 HC
## AAACCCAGTGAGTTTC-NS3R189BTS 21.675930 0.00000000 HC
## AAACCCAGTTCTTAGG-NS3R189BTS 23.353124 0.00000000 HC1.5. Generate before and after QC plots
# 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,"before-qc-violineplot.pdf",sep=""), plot = qcbefore, width = 15, height = 20, units = "cm")
# 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,"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.
# Ask SCassist to recommend normalization method
normalization_recommendation<-SCassist_recommend_normalization("allsamplesgood", llm_server="google", api_key_file = api_key_file)## ## Recommended Normalization Method: SCTransform
##
## Based on the provided characteristics of your 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 34,967 cells, computational efficiency becomes a significant factor, and SCTransform excels in this regard.
##
## **2. High Variability in Gene Expression:** The high standard deviation of 6.34 in gene expression values indicates significant variability across cells. SCTransform accounts for this variability by modeling both technical and biological sources of variation, leading to more accurate normalization.
##
## **3. Moderate Library Size Variation:** The coefficient of variation of 0.58 for library sizes suggests moderate variation in sequencing depth. SCTransform effectively addresses this by regressing out library size effects during normalization.
##
## **Why SCTransform is Preferred:**
##
## * **Comprehensive Normalization:** SCTransform combines several normalization steps, including library size correction, removal of unwanted variation (e.g., cell cycle effects), and gene-specific scaling. This comprehensive approach leads to more accurate and robust results.
## * **Model-Based Approach:** SCTransform utilizes a statistical model to estimate and correct for technical biases, resulting in more reliable normalization compared to simple scaling methods.
## * **Improved Downstream Analysis:** SCTransform-normalized data often leads to better performance in downstream analyses like dimensionality reduction, clustering, and differential gene expression analysis.
##
## **Potential Alternatives:**
##
## * **LogNormalize:** This method is a simple scaling method that divides each gene count by the total number of counts in a cell and then takes the logarithm. While it is computationally efficient, it does not account for technical biases or gene-specific variability.
## * **CLR (Centered Log Ratio):** This method is based on the log ratio of gene counts to the geometric mean of all genes in a cell. It is effective for removing library size effects but does not address other sources of variation.
## * **RC (Relative Counts):** This method scales gene counts by the total number of counts in a cell. It is a simple method but does not account for technical biases or gene-specific variability.
##
## **Conclusion:**
##
## While other normalization methods might be suitable for smaller datasets or datasets with less variability, SCTransform is the most appropriate choice for this dataset due to its ability to handle large datasets, account for technical and biological variation, and improve 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
# List top 10 variable genes
top10 <- head(VariableFeatures(allsamplesgood), 10)
top10## [1] "IGKC" "IGLC1" "IGLC2" "IGHA1" "IGLC3" "PPBP" "S100A9" "IGHM"
## [9] "GNLY" "CCL4"# 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,"variablefeatureplot.pdf",sep=""), plot = vfp1, width = 20, height = 15, units = "cm")3.1. Analyze Variable Features
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 = "pbmcs of one healthy control and one uveitis patient"
# Ask SCassist to analyze variable features
variable_feature_analysis <- SCassist_analyze_variable_features("allsamplesgood", experimental_design = experimental_design, api_key_file = api_key_file, llm_server="google")## The provided list of genes suggests an enrichment in pathways related to **immune response, inflammation, and leukocyte activation**.
##
## Here's a breakdown:
##
## **Immune Response & Inflammation:**
##
## * **Immunoglobulin genes (IGKC, IGLC1, IGLC2, IGHA1, IGLC3, IGHM):** These genes encode components of antibodies, crucial for adaptive immune responses. Their presence suggests activation of B cells and antibody production.
## * **Cytokines (IL1B, CCL4, CCL3, CCL4L2, CXCL8):** These genes encode signaling molecules that regulate immune cell activity and inflammation. Their presence indicates a pro-inflammatory environment.
## * **HLA genes (HLA-DRA, HLA-DRB1):** These genes encode major histocompatibility complex (MHC) proteins, crucial for antigen presentation to T cells. Their presence suggests activation of T cell responses.
## * **S100 proteins (S100A8, S100A9):** These proteins are involved in inflammation and immune cell activation. Their presence suggests a pro-inflammatory state.
## * **JCHAIN:** This gene encodes a protein essential for the assembly of multimeric immunoglobulins, further supporting antibody production.
## * **CD74:** This gene encodes a receptor for the cytokine MIF, which is involved in inflammation and immune regulation.
## * **IFIT2:** This gene encodes an interferon-induced protein involved in antiviral responses and inflammation.
##
## **Leukocyte Activation & Function:**
##
## * **GNLY:** This gene encodes a cytotoxic granule protein found in NK cells and cytotoxic T cells, suggesting activation of these cell types.
## * **PPBP:** This gene encodes a platelet factor that promotes platelet aggregation and activation.
## * **PF4:** This gene encodes a chemokine that attracts and activates platelets.
## * **GP1BB:** This gene encodes a platelet receptor involved in platelet aggregation.
## * **LYZ:** This gene encodes lysozyme, an enzyme found in neutrophils and macrophages that plays a role in bacterial killing.
## * **CST3:** This gene encodes cystatin C, a cysteine protease inhibitor involved in immune regulation.
## * **CD79A:** This gene encodes a protein involved in B cell receptor signaling, suggesting B cell activation.
##
## **Other Relevant Genes:**
##
## * **NRGN:** This gene encodes a protein involved in neuronal development, but its presence in this context could suggest potential neuro-inflammatory aspects.
## * **PTGDS:** This gene encodes a prostaglandin synthase involved in inflammation and pain.
## * **EREG:** This gene encodes a growth factor involved in wound healing and inflammation.
## * **HSPA1A:** This gene encodes a heat shock protein involved in cellular stress response.
##
## **Relevance to the Experimental Design:**
##
## The presence of these genes in PBMCs from a healthy control and a uveitis patient suggests that the patient's immune system is activated and potentially dysregulated. The specific genes involved point towards a pro-inflammatory environment with activation of both innate and adaptive immune responses. This is consistent with the known pathology of uveitis, an inflammatory disease affecting the eye.
##
## The comparison between the healthy control and the uveitis patient allows for the identification of genes that are differentially expressed in the patient, potentially providing insights into the specific mechanisms underlying the disease.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: S100A9, LYZ, S100A8, IL1B, HLA-DRA
## Negative: RPS26, RPS27, RPS12, EEF1A1, RPL10
## PC_ 2
## Positive: IGKC, IGHA1, IGLC2, IGLC1, JCHAIN
## Negative: S100A9, S100A8, PPBP, LYZ, IL1B
## PC_ 3
## Positive: PPBP, PF4, NRGN, GP1BB, GNG11
## Negative: IGLC1, GNLY, IL1B, CXCL8, CCL3
## PC_ 4
## Positive: IGKC, RPS26, HLA-DRA, CD74, RPS27
## Negative: IGLC1, IGLC3, IGHM, PPBP, JCHAIN
## PC_ 5
## Positive: IGLC1, HLA-DRA, IGLC3, CD74, RPS26
## Negative: GNLY, NKG7, IGKC, GZMB, CCL44.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
pc_recommendation=SCassist_recommend_pcs("allsamplesgood", experimental_design = experimental_design, llm_server="google", api_key_file = api_key_file)## 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 the 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 65% of the total variance. This is a good balance between capturing a significant portion of the variation while avoiding excessive complexity.
## * **Balance between complexity and interpretability:** Using 10 PCs provides a reasonable level of detail without making the analysis overly complex. Including more PCs might capture subtle variations but could lead to overfitting and difficulty in interpreting the results.
##
## Therefore, using the first 10 PCs strikes a good balance between capturing a significant portion of the variance and maintaining interpretability.# 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 HC vs BCR-UV
Idents(allsamplesgood) <- "orig.ident"
# Visualize the cells after pca
DimPlot(object = allsamplesgood, reduction = "pca")
# 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)
4.2. SCassist analyze PCs
Here we ask SCassist_analyze_pcs to analzye the SCassist recommended PCs and provide insights on the gene sets identified
pc_analyzed=SCassist_analyze_pcs("allsamplesgood", num_pcs = 5, experimental_design = experimental_design, llm_server="google", api_key_file = api_key_file)##
## **PC1 Summary:**
## The top contributing genes for PC1 are primarily associated with immune response and inflammation. Genes like S100A9, LYZ, S100A8, IL1B, HLA-DRA, CST3, CXCL8, CD74, and IFI30 are all known to be involved in the activation and recruitment of immune cells, particularly neutrophils and macrophages. Additionally, genes like KYNU, TYROBP, and FCN1 are involved in antigen presentation and immune signaling. The presence of these genes suggests that PC1 might be capturing variations in the immune response, potentially reflecting differences in the activation state of immune cells between the healthy control and the uveitis patient. This is further supported by the presence of genes like ZEB2, SPI1, and NRGN, which are involved in the regulation of immune cell differentiation and function. Therefore, PC1 likely reflects differences in the immune response and inflammatory processes between the healthy control and the uveitis patient.
##
##
## **PC2 Summary:**
## The top contributing genes for PC2 are primarily involved in immune response and inflammation, particularly those related to B cell activation and antibody production. Genes like IGKC, IGHA1, IGLC2, IGLC1, JCHAIN, IGHM, and CD79A are components of immunoglobulin molecules and B cell signaling pathways. Additionally, genes like HLA-DRA, HLA-DQA1, HLA-DRB1, and MS4A1 are involved in antigen presentation and T cell activation. The presence of genes like IL1B, CXCL8, and S100A8 suggests an inflammatory response. These findings suggest that PC2 might be capturing variations in B cell activation and antibody production, potentially reflecting differences in immune response between the healthy control and the uveitis patient. This could be related to the pathogenesis of uveitis, an inflammatory condition affecting the eye.
##
##
## **PC3 Summary:**
## The top contributing genes for PC3 include several chemokines (CCL3, CCL4, CXCL8), inflammatory mediators (IL1B, PTGS1), and genes involved in immune cell activation and differentiation (GNLY, ITGA2B, CD9, TSPAN33). These genes are strongly associated with immune responses, particularly those related to inflammation and leukocyte recruitment. The presence of genes like PPBP, PF4, and SPARC suggests potential involvement of platelet activation and coagulation. Furthermore, genes like TUBB1, CAVIN2, and TLN1 point towards cytoskeletal remodeling and cell migration. Taken together, these findings suggest that PC3 might be capturing variations related to **immune cell activation, inflammation, and leukocyte trafficking**, potentially reflecting differences in immune responses between the healthy control and the uveitis patient.
##
##
## **PC4 Summary:**
## The top contributing genes for PC4 are primarily involved in immune response and cell signaling. Genes like IGLC1, IGKC, IGHM, and JCHAIN are components of immunoglobulin production, suggesting a strong B cell signature. Other genes like GNLY, PF4, GZMB, and NKG7 are associated with cytotoxic T cell activity and NK cell function. Furthermore, genes like CCL4 and CCL5 are chemokines involved in attracting immune cells to sites of inflammation. The presence of these genes suggests that PC4 might be capturing variations in the immune cell composition and activation state, potentially reflecting differences in the immune response between the healthy control and the uveitis patient.
##
##
## **PC5 Summary:**
## The top contributing genes for PC5 are primarily involved in immune response and inflammation, particularly those related to cytotoxic T lymphocytes (CTLs) and NK cells. Genes like GNLY, NKG7, GZMB, PRF1, and GZMA are known effector molecules of CTLs and NK cells, responsible for killing target cells. Additionally, genes like CCL4, CCL5, and IFNG are chemokines and cytokines that play a role in attracting and activating immune cells. The presence of HLA genes (HLA-DRA, HLA-DRB1, HLA-DPA1, HLA-DPB1, HLA-DQA1) suggests differences in antigen presentation and immune recognition between the healthy control and uveitis patient. Therefore, PC5 likely captures variations in the immune response, potentially reflecting the activation of CTLs and NK cells, and the associated inflammatory processes, which could be relevant to the pathogenesis of uveitis.
##
##
## ## Overall Summary of PCs:
##
## The five principal components (PCs) identified in this analysis all highlight significant variations in immune response and inflammation between healthy controls and uveitis patients.
##
## **Common Themes:**
##
## * **Immune Response:** All PCs show a strong enrichment of genes involved in immune response, particularly those related to B cell activation, antibody production, T cell activation, and cytotoxic T lymphocyte (CTL) and NK cell activity.
## * **Inflammation:** Genes associated with inflammation, including chemokines, inflammatory mediators, and genes involved in leukocyte recruitment, are consistently present across all PCs.
##
## **Specific Contributions of Each PC:**
##
## * **PC1:** Focuses on the activation and recruitment of immune cells, particularly neutrophils and macrophages.
## * **PC2:** Highlights variations in B cell activation and antibody production, potentially reflecting differences in immune response related to uveitis pathogenesis.
## * **PC3:** Captures variations in immune cell activation, inflammation, and leukocyte trafficking, suggesting differences in immune responses between groups.
## * **PC4:** Shows a strong B cell signature, along with cytotoxic T cell activity and NK cell function, indicating variations in immune cell composition and activation state.
## * **PC5:** Emphasizes the activation of CTLs and NK cells, and associated inflammatory processes, potentially relevant to uveitis pathogenesis.
##
## **Overall, the analysis suggests that uveitis is associated with distinct patterns of immune response and inflammation, potentially involving different immune cell populations and activation states.** Further investigation is needed to understand the specific roles of each PC in the pathogenesis of uveitis and to identify potential therapeutic targets.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
recommended_k <- SCassist_recommend_k("allsamplesgood", num_pcs = 10, llm_server="google", api_key_file = api_key_file)## ## Recommended K: 10-50
##
## ## Reasoning:
##
## The `k.param` value in `FindNeighbors()` determines the number of nearest neighbors considered for each cell when constructing the k-nearest neighbor graph. This graph is then used for clustering. A higher `k.param` value leads to a more connected graph, potentially merging distinct populations. Conversely, a lower value might result in overly fragmented clusters.
##
## Given your dataset size of 34,967 cells and the use of 10 principal components, a range of `k.param` values between 10 and 50 is recommended. This range allows for exploration of different levels of connectivity while considering the dimensionality of your data. Starting with a lower value like 10 might reveal finer subpopulations, while increasing the value to 50 could capture broader cell types.# Run findneighbors, with above identified pcs and k
allsamplesgood <- FindNeighbors(allsamplesgood, dims = 1:10, k.param = 30, 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", llm_server="google", api_key_file = api_key_file)## Based on the data characteristics, I recommend:
##
## **Recommended Resolution:** seq(0.2, 1.2, 0.1)
##
## **Reasoning:**
##
## The mean expression variability of 0.868717336264306 suggests a moderate level of heterogeneity in your dataset. This indicates that there are likely distinct cell populations, but they may not be drastically different from each other. The median neighbor distance of 2.50502443313599 in the k-nearest neighbor graph further supports this notion, as it implies a moderate level of separation between cells.
##
## Therefore, a resolution range of 0.2 to 1.2 with increments of 0.1 is recommended. This range allows for the identification of both distinct and subtle cell populations. Lower resolutions (closer to 0.2) will capture broader, more general cell types, while higher resolutions (closer to 1.2) will reveal finer distinctions within those populations. By exploring this range, you can effectively identify the optimal resolution for your specific analysis goals.allsamplesgood <- FindNeighbors(allsamplesgood, dims = 1:10, k.param = 30)
# Perform clustering using the SCassist recommended resolution
allsamplesgood <- FindClusters(allsamplesgood, resolution = seq(0.2,1.2,0.1))## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9408
## Number of communities: 10
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9242
## Number of communities: 12
## Elapsed time: 10 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9123
## Number of communities: 13
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9039
## Number of communities: 15
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8956
## Number of communities: 17
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8892
## Number of communities: 17
## Elapsed time: 7 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8819
## Number of communities: 20
## Elapsed time: 8 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8761
## Number of communities: 21
## Elapsed time: 7 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8707
## Number of communities: 23
## Elapsed time: 6 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8667
## Number of communities: 25
## Elapsed time: 7 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 34967
## Number of edges: 1730396
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8610
## Number of communities: 28
## Elapsed time: 7 secondshead(allsamplesgood)## orig.ident nCount_RNA nFeature_RNA percent.mito
## AAACCCACAAGTCATC-NS3R189BTS NS3R189BTS 2431 1431 4.1135335
## AAACCCACACTGGAAG-NS3R189BTS NS3R189BTS 2676 1078 0.8221226
## AAACCCACAGCACACC-NS3R189BTS NS3R189BTS 3394 1760 5.1856217
## AAACCCACAGGAGGAG-NS3R189BTS NS3R189BTS 2160 1279 6.5277778
## AAACCCACATGACGTT-NS3R189BTS NS3R189BTS 2329 1271 4.3795620
## AAACCCAGTAGTTACC-NS3R189BTS NS3R189BTS 5362 2261 3.8418501
## AAACCCAGTCCAAATC-NS3R189BTS NS3R189BTS 7926 2660 4.6177139
## AAACCCAGTCCATACA-NS3R189BTS NS3R189BTS 2239 1283 5.0468959
## AAACCCAGTGAGTTTC-NS3R189BTS NS3R189BTS 2339 1291 5.3014109
## AAACCCAGTTCTTAGG-NS3R189BTS NS3R189BTS 3537 1608 2.2900763
## percent.ribo percent.hb type nCount_SCT
## AAACCCACAAGTCATC-NS3R189BTS 18.510901 0.00000000 HC 2476
## AAACCCACACTGGAAG-NS3R189BTS 42.526158 0.00000000 HC 2660
## AAACCCACAGCACACC-NS3R189BTS 19.151444 0.00000000 HC 3016
## AAACCCACAGGAGGAG-NS3R189BTS 8.842593 0.00000000 HC 2321
## AAACCCACATGACGTT-NS3R189BTS 18.978102 0.04293688 HC 2431
## AAACCCAGTAGTTACC-NS3R189BTS 18.985453 0.01864976 HC 3350
## AAACCCAGTCCAAATC-NS3R189BTS 30.052990 0.02523341 HC 2781
## AAACCCAGTCCATACA-NS3R189BTS 19.830281 0.04466280 HC 2379
## AAACCCAGTGAGTTTC-NS3R189BTS 21.675930 0.00000000 HC 2440
## AAACCCAGTTCTTAGG-NS3R189BTS 23.353124 0.00000000 HC 3036
## nFeature_SCT SCT_snn_res.0.2 SCT_snn_res.0.3
## AAACCCACAAGTCATC-NS3R189BTS 1429 0 0
## AAACCCACACTGGAAG-NS3R189BTS 1078 1 6
## AAACCCACAGCACACC-NS3R189BTS 1758 0 0
## AAACCCACAGGAGGAG-NS3R189BTS 1278 3 3
## AAACCCACATGACGTT-NS3R189BTS 1271 6 8
## AAACCCAGTAGTTACC-NS3R189BTS 2104 4 4
## AAACCCAGTCCAAATC-NS3R189BTS 1373 2 5
## AAACCCAGTCCATACA-NS3R189BTS 1283 0 0
## AAACCCAGTGAGTTTC-NS3R189BTS 1291 0 0
## AAACCCAGTTCTTAGG-NS3R189BTS 1608 2 5
## SCT_snn_res.0.4 SCT_snn_res.0.5 SCT_snn_res.0.6
## AAACCCACAAGTCATC-NS3R189BTS 0 0 0
## AAACCCACACTGGAAG-NS3R189BTS 2 2 1
## AAACCCACAGCACACC-NS3R189BTS 0 0 0
## AAACCCACAGGAGGAG-NS3R189BTS 6 4 4
## AAACCCACATGACGTT-NS3R189BTS 9 11 13
## AAACCCAGTAGTTACC-NS3R189BTS 5 6 7
## AAACCCAGTCCAAATC-NS3R189BTS 4 7 3
## AAACCCAGTCCATACA-NS3R189BTS 0 0 0
## AAACCCAGTGAGTTTC-NS3R189BTS 0 0 7
## AAACCCAGTTCTTAGG-NS3R189BTS 4 7 3
## SCT_snn_res.0.7 SCT_snn_res.0.8 SCT_snn_res.0.9
## AAACCCACAAGTCATC-NS3R189BTS 0 0 5
## AAACCCACACTGGAAG-NS3R189BTS 9 13 4
## AAACCCACAGCACACC-NS3R189BTS 0 0 5
## AAACCCACAGGAGGAG-NS3R189BTS 6 4 8
## AAACCCACATGACGTT-NS3R189BTS 13 15 15
## AAACCCAGTAGTTACC-NS3R189BTS 4 11 17
## AAACCCAGTCCAAATC-NS3R189BTS 3 3 2
## AAACCCAGTCCATACA-NS3R189BTS 0 0 0
## AAACCCAGTGAGTTTC-NS3R189BTS 5 11 0
## AAACCCAGTTCTTAGG-NS3R189BTS 3 3 2
## SCT_snn_res.1 SCT_snn_res.1.1 SCT_snn_res.1.2
## AAACCCACAAGTCATC-NS3R189BTS 7 3 3
## AAACCCACACTGGAAG-NS3R189BTS 2 4 4
## AAACCCACAGCACACC-NS3R189BTS 7 3 3
## AAACCCACAGGAGGAG-NS3R189BTS 15 15 19
## AAACCCACATGACGTT-NS3R189BTS 16 18 17
## AAACCCAGTAGTTACC-NS3R189BTS 12 20 21
## AAACCCAGTCCAAATC-NS3R189BTS 5 17 16
## AAACCCAGTCCATACA-NS3R189BTS 0 6 2
## AAACCCAGTGAGTTTC-NS3R189BTS 12 6 2
## AAACCCAGTTCTTAGG-NS3R189BTS 5 17 16
## seurat_clusters
## AAACCCACAAGTCATC-NS3R189BTS 3
## AAACCCACACTGGAAG-NS3R189BTS 4
## AAACCCACAGCACACC-NS3R189BTS 3
## AAACCCACAGGAGGAG-NS3R189BTS 19
## AAACCCACATGACGTT-NS3R189BTS 17
## AAACCCAGTAGTTACC-NS3R189BTS 21
## AAACCCAGTCCAAATC-NS3R189BTS 16
## AAACCCAGTCCATACA-NS3R189BTS 2
## AAACCCAGTGAGTTTC-NS3R189BTS 2
## AAACCCAGTTCTTAGG-NS3R189BTS 16# 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.2 SCT_snn_res.0.3 SCT_snn_res.0.4 SCT_snn_res.0.5 SCT_snn_res.0.6
## 10 12 13 15 17
## SCT_snn_res.0.7 SCT_snn_res.0.8 SCT_snn_res.0.9 SCT_snn_res.1 SCT_snn_res.1.1
## 17 20 21 23 25
## SCT_snn_res.1.2
## 28# change default identity
Idents(allsamplesgood) <- "SCT_snn_res.0.5"
# list cell number in each cluster for HC vs UV
table(Idents(allsamplesgood),allsamplesgood$type)##
## BCR-UV HC
## 0 58 6404
## 1 3953 433
## 2 481 3506
## 3 3029 29
## 4 1467 1333
## 5 521 2193
## 6 887 1710
## 7 17 2129
## 8 890 885
## 9 43 1308
## 10 1020 68
## 11 3 923
## 12 9 817
## 13 389 84
## 14 5 373# 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 HC vs BCR-UV
DimPlot(object = allsamplesgood, pt.size=0.5, reduction = "umap", group.by = "type")
# Plot umap
dimplot1<-DimPlot(object = allsamplesgood, pt.size=0.5, reduction = "umap", label = T)
# Color by HC vs BCR-UV
dimplot2<-DimPlot(object = allsamplesgood, pt.size=0.5, reduction = "umap", group.by = "type")
# Save umaps
ggsave(paste(outputpath,"umap-clusters.pdf",sep=""), plot = dimplot1, width = 20, height = 15, units = "cm")
ggsave(paste(outputpath,"umap-clusters-type.pdf",sep=""), plot = dimplot2, width = 20, height = 15, units = "cm")
# Custom list of genes to visualize
mygenes = c("NKG7","CCL5")
RidgePlot(allsamplesgood, features = mygenes)
VlnPlot(allsamplesgood, features = mygenes, pt.size=0)
FeaturePlot(allsamplesgood, features = mygenes)