T-cell receptor selection for TCR-T studies from antigen specific culture and scRNA/VDJ sequencing
library(ggplot2)
library(dplyr)
library(Seurat)
library(patchwork)
so_data <- Read10X(data.dir = "filtered_feature_bc_matrix", gene.column=1)
so <- CreateSeuratObject(counts = so_data, project = "H37-H37004-P13-Vin_5GEX_P13", min.cells = 3, min.features = 200)
so
tcells <- subset(so, (ENSG00000153563.16 > 0 | ENSG00000172116.23 > 0) & ENSG00000010610.10 > 0)
cd8s <- subset(so, (ENSG00000153563.16 > 0 | ENSG00000172116.23 > 0) & ENSG00000010610.10 == 0)
cd4s <- subset(so, (ENSG00000153563.16 == 0 & ENSG00000172116.23 == 0) & ENSG00000010610.10 > 0)
tcells_cellNames <- rownames(tcells@meta.data)
cd8s_cellNames <- rownames(cd8s@meta.data)
cd4s_cellNames <- rownames(cd4s@meta.data)
cd8s_cd8a_raw_counts <- FetchData(object = cd8s, vars="ENSG00000153563.16", layer = "counts")
quantile(cd8s_cd8a_raw_counts$ENSG00000153563.16)
cd8s_cd8b_raw_counts <- FetchData(object = cd8s, vars="ENSG00000172116.23", layer = "counts")
quantile(cd8s_cd8b_raw_counts$ENSG00000172116.23)
cd4s_cd4_raw_counts <- FetchData(object = cd4s, vars="ENSG00000010610.10", layer = "counts")
quantile(cd4s_cd4_raw_counts$ENSG00000010610.10)
so[['exp_annot']] = rep('NA', ncol(so))
so@meta.data[['exp_annot']] <- so@meta.data %>% mutate(exp_annot = case_when(rownames(so@meta.data) %in% tcells_cellNames ~ 'Tcell', rownames(so@meta.data) %in% cd8s_cellNames ~ 'CD8+', rownames(so@meta.data) %in% cd4s_cellNames ~ 'CD4+')) %>% select(exp_annot)
so@meta.data[['CD8A_reads']] <- FetchData(so, vars='ENSG00000153563.16')
so@meta.data[['CD8B_reads']] <- FetchData(so, vars='ENSG00000172116.23')
so@meta.data[['CD4_reads']] <- FetchData(so, vars='ENSG00000010610.10')
so@meta.data[['CD3D_reads']] <- FetchData(so, vars='ENSG00000167286.10')
so@meta.data[['CD3E_reads']] <- FetchData(so, vars='ENSG00000198851.10')
so@meta.data[['CD3G_reads']] <- FetchData(so, vars='ENSG00000160654.11')
so
so[['percent.mt']] <- PercentageFeatureSet(so, pattern = 'ENSG00000210049.1|ENSG00000211459.2|ENSG00000210077.1|ENSG00000210082.2|ENSG00000209082.1|ENSG00000198888.2|ENSG00000210100.1|ENSG00000210107.1|ENSG00000210112.1|ENSG00000198763.3|ENSG00000210117.1|ENSG00000210127.1|ENSG00000210135.1|ENSG00000210140.1|ENSG00000210144.1|ENSG00000198804.2|ENSG00000210151.2|ENSG00000210154.1|ENSG00000198712.1|ENSG00000210156.1|ENSG00000228253.1|ENSG00000198899.2|ENSG00000198938.2|ENSG00000210164.1|ENSG00000198840.2|ENSG00000210174.1|ENSG00000212907.2|ENSG00000198886.2|ENSG00000210176.1|ENSG00000210184.1|ENSG00000210191.1|ENSG00000198786.2|ENSG00000198695.2|ENSG00000210194.1|ENSG00000198727.2|ENSG00000210195.2|ENSG00000210196.2')
p <- VlnPlot(so, features = c('nFeature_RNA', 'nCount_RNA', 'percent.mt'), ncol = 3)
ggsave("H37-H37004-P13-Vin_5GEX_P13.vln_plot.pdf", p)
plot1 <- FeatureScatter(so, feature1 = 'nCount_RNA', feature2 = 'percent.mt')
plot2 <- FeatureScatter(so, feature1 = 'nCount_RNA', feature2 = 'nFeature_RNA')
jp <- plot1 + plot2
ggsave("H37-H37004-P13-Vin_5GEX_P13.scatter_plot.pdf", jp)
so <- subset(so, subset = nFeature_RNA > 200 & percent.mt < 10)
so <- NormalizeData(so, assay="RNA", normalization.method = "LogNormalize", scale.factor = 10000)
so <- FindVariableFeatures(so, selection.method = 'vst', nfeatures = 2000)
top10 <- head(VariableFeatures(so), 10)
plot1 <- FeatureScatter(so, feature1 = 'nCount_RNA', feature2 = 'percent.mt')
ggsave("H37-H37004-P13-Vin_5GEX_P13.var_plot.pdf", plot1)
all.genes <- rownames(so)
so <- ScaleData(so, features = all.genes)
so <- RunPCA(so, assay="RNA", npcs = 50, features = VariableFeatures(object = so))
print(so[['pca']], dims = 1:5, nfeatures = 5)
dim_loads <- VizDimLoadings(so, dims = 1:10, reduction = 'pca')
ggsave("H37-H37004-P13-Vin_5GEX_P13.dim_loads.pdf", dim_loads)
dim_heatmap <- DimHeatmap(so, dims = 1:15, cells = 500, balanced = TRUE)
ggsave("H37-H37004-P13-Vin_5GEX_P13.dim_heatmap.pdf", dim_heatmap)
so <- JackStraw(so, num.replicate = 100)
so <- ScoreJackStraw(so, dims = 1:20)
jackstraw_plot <- JackStrawPlot(so, dims = 1:15)
ggsave("H37-H37004-P13-Vin_5GEX_P13.jackstraw.pdf", jackstraw_plot)
elbow_plot <- ElbowPlot(so, ndims = 30)
ggsave("H37-H37004-P13-Vin_5GEX_P13.elbow_plot.pdf", elbow_plot)
ndims = length(which(so@reductions$pca@stdev > 2))
so <- FindNeighbors(so, dims = 1:ndims)
so <- FindClusters(so, resolution = 1.2)
so <- RunUMAP(so, dims = 1:ndims)
so.markers <- FindAllMarkers(so, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
so.markers %>% group_by(cluster) %>% slice_max(n = 2, order_by = avg_log2FC)
so <- RunUMAP(so, dims = 1:10)
umap_plot <- DimPlot(so, reduction = 'umap')
ggsave("H37-H37004-P13-Vin_5GEX_P13.umap.pdf", umap_plot)
ad <- GetAssayData(so, assay = "RNA", slot = "data")
saveRDS(so, file = "H37-H37004-P13-Vin_5GEX_P13.seurat.rds")
saveRDS(ad, file = "H37-H37004-P13-Vin_5GEX_P13.seurat.assay_data.rds")