scRepertorie+Seurat

ref: https://ncborcherding.github.io/vignettes/vignette.html#6_Interacting_with_Seurat
data: download.file("https://bioshare.bioinformatics.ucdavis.edu/bioshare/download/iimg5mz77whzzqc/vdj_v1_mm_balbc_pbmc.zip"）

#scRepertoire
rm(list=ls())
library(scRepertoire)
csv1 <- read.csv("vdj_v1_mm_balbc_pbmc_t_filtered_contig_annotations.csv", stringsAsFactors = F)
contig_list <- list(csv1)
combined <- combineTCR(contig_list, samples = c("balbc"), ID = c("pbmc"),cells = "T-AB")
#the total or relative numbers of unique clonotypes. 
p1 <- quantContig(combined, cloneCall="gene+nt", scale = T)
#the relative distribution of clonotypes
p2 <- abundanceContig(combined, cloneCall = "gene", scale = F)
#the length distribution of the CDR3 sequences
p3 <- lengthContig(combined, cloneCall="aa", chains = "combined") 
#the individual chains
p4 <- lengthContig(combined, cloneCall="nt", chains = "single") 
#the relative usage of vgenes of the TCR
p5 <- vizVgenes(combined, TCR="TCR1", facet.x = "sample", facet.y = "ID")
library(patchwork)
p1+p2+p3+p4
p5

#Clonal Space Homeostasis
p6 <- clonalHomeostasis(combined, cloneCall = "gene")
#Clonal Proportion
p7 <- clonalProportion(combined, cloneCall = "gene") 
p6+p7

#Interacting with Seurat
library(Seurat)
library(cowplot)
library(hdf5r)
#加载Cellranger output file
balbc_pbmc <- Read10X_h5("vdj_v1_mm_balbc_pbmc_5gex_filtered_feature_bc_matrix.h5")
s_balbc_pbmc <- CreateSeuratObject(counts = balbc_pbmc, min.cells = 3, min.features = 200, project = "cellranger")
#提取线粒体基因
s_balbc_pbmc$percent.mito <- PercentageFeatureSet(s_balbc_pbmc, pattern = "^mt-")
#seurat workflow
seurat_workflow <- function(s_balbc_pbmc){
s_balbc_pbmc <- subset(s_balbc_pbmc, percent.mito <= 10)
s_balbc_pbmc <- subset(s_balbc_pbmc, nCount_RNA >= 500 & nCount_RNA <= 40000)
s_balbc_pbmc <- NormalizeData(s_balbc_pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
s_balbc_pbmc <- FindVariableFeatures(s_balbc_pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(s_balbc_pbmc)
s_balbc_pbmc <- ScaleData(s_balbc_pbmc, features = all.genes)
s_balbc_pbmc <- RunPCA(s_balbc_pbmc, features = VariableFeatures(object = s_balbc_pbmc))
use.pcs = 1:30
s_balbc_pbmc <- FindNeighbors(s_balbc_pbmc, dims = use.pcs)
s_balbc_pbmc <- FindClusters(s_balbc_pbmc, resolution = c(0.5))
s_balbc_pbmc <- RunUMAP(s_balbc_pbmc, dims = use.pcs)}
s_balbc_pbmc <- seurat_workflow(s_balbc_pbmc)
seurat <- s_balbc_pbmc
#save(s_balbc_pbmc,file="seurat.rda")
#seurat <- get(load("seurat.rda"))
#Becouse the barcodes in Seurat have no prefix, we need remove the barcode prefix in combined
for (i in seq_along(combined)) {
  combined[[i]] <- stripBarcode(combined[[i]], column = 1, connector = "_", num_connects = 3)
}
seurat <- combineExpression(combined, seurat,cloneCall = "gene")
#visualization
colorblind_vector <- colorRampPalette(c("#FF4B20", "#FFB433", "#C6FDEC", "#7AC5FF", "#0348A6","#0348A6","#0348A6"))
#the distribution of the clonotype bins by first ordering the clonoType as a factor
slot(seurat, "meta.data")$cloneType <- factor(slot(seurat, "meta.data")$cloneType, 
levels = c("Hyperexpanded (100 < X <= 500)", 
"Large (20 < X <= 100)",  
"Medium (5 < X <= 20)", 
"Small (1 < X <= 5)",  
"Single (0 < X <= 1)", NA))
p1 <- DimPlot(seurat, group.by = "cloneType") + scale_color_manual(values = colorblind_vector(5), na.value="grey")
#look at the clonotypes by calling specific sequences
seurat <- highlightClonotypes(seurat, cloneCall= "aa", sequence = c("CAARDTGYQNFYF_CASSIRVNTEVFF"))
p2 <- DimPlot(seurat, group.by = "highlight")
 ##type1数量太少了，不是很明显
#look at count of cells by cluster assigned into specific frequency ranges
p3 <- occupiedscRepertoire(seurat, x.axis = "cluster")
p1/p2/p3

#clonaldiversity
combined2 <- expression2List(seurat, group = "cluster")
length(combined2) #now listed by cluster
p4 <- clonalDiversity(combined2, cloneCall = "nt")
#clonalhomeostasis
p5 <- clonalHomeostasis(combined2, cloneCall = "nt")
#clonalpropoetion
p6 <- clonalProportion(combined2, cloneCall = "nt")
#clonaloverlap
p7 <- clonalOverlap(combined2, cloneCall="aa", method="overlap")
p4/p5/p6
p7