pre-process

1 Functions

tpm <- function(counts, lengths) {
  rpk <- counts / lengths
  coef <- sum(rpk) / 1e6
  rpk/coef
}

rpk <- function(counts, lengths) {
  rpk <- counts / lengths
  rpk
}

2 Read data

acc_immune_data_part1 <- read.delim("./Data/raw_data/SN0276364",skip=1,header=T, 
                       sep="\t",stringsAsFactors=F,row.names=1,check.names = F)

acc_immune_data_part2 <- read.delim("./Data/raw_data/SN0276892",skip=1,header=T, 
                       sep="\t",stringsAsFactors=F,row.names=1,check.names = F)

3 Pre-process

#save before omitting
genes_names = acc_immune_data_part1$gene_name
genes_lengths = acc_immune_data_part1[,5]

#omit non relevant cols
acc_immune_data_part1=acc_immune_data_part1[,7:ncol(acc_immune_data_part1)]
acc_immune_data_part2=acc_immune_data_part2[,7:ncol(acc_immune_data_part2)]

#combine datasets
acc_immune_data = cbind(acc_immune_data_part1,acc_immune_data_part2)

#create genes names
genes_names=make.unique(genes_names) %>% replace_na('NA')
rownames(acc_immune_data) = genes_names

#omit non relevant genes
omitgenes= startsWith(rownames(acc_immune_data),"NA")
acc_immune_data=acc_immune_data[!omitgenes,]
genes_lengths = genes_lengths[!omitgenes] #update genes_lengths

#calcualte gene length and MT genes
mt_genes = startsWith(rownames(acc_immune_data),"MT-")| startsWith(rownames(acc_immune_data),"ERCC-")

#get colnames
cell.labels <- gsub(pattern = "/.*$",replacement = "",colnames(acc_immune_data))

#change colnames
colnames(acc_immune_data) <- cell.labels


acc_immune_counts <- CreateSeuratObject(counts = acc_immune_data, project = "acc_immune_counts", min.cells = 3, min.features = 1000)

Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')

4 QC

acc_immune_counts@meta.data[["percent.mt"]] <- PercentageFeatureSet(acc_immune_counts, pattern = "^MT-")
print_tab(plt = 
            FeatureScatter(acc_immune_counts, feature1 = "nCount_RNA", feature2 = "percent.mt") + 
            theme(legend.position="none", axis.text.x = element_text(size=8)) + 
              geom_point(color='darkblue')
  
            ,title = "MT percentages")

MT percentages

print_tab(plt = 
            VlnPlot(acc_immune_counts, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
,title = "violin plots")

violin plots

NA

acc_immune_data.rpk=apply(acc_immune_data[!mt_genes,], 2, function(x) rpk(x, genes_lengths[!mt_genes]))
acc_immune <- CreateSeuratObject(counts = acc_immune_data.rpk, project = "acc_immune", min.cells = 3, min.features = 1000)

Warning: Feature names cannot have underscores ('_'), replacing with dashes ('-')

acc_immune = NormalizeData(object = acc_immune,scale.factor = 1e6) #create TPM

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

acc_immune@assays$RNA@data = log2(exp(1)) * acc_immune@assays$RNA@data #convert from e base to 2 base log

# new_genes = rownames(acc_immune_data) %in% rownames(acc_immune)
# acc_immune_data.tpm=apply(acc_immune_data[new_genes,], 2, function(x) tpm(x, genes_lengths[new_genes])) %>%as.data.frame()
# acc_immune_data.tpm = log2(acc_immune_data.tpm+1)

5 Filtering

nFeature_RNA_threshold = 1000
percent.mt_threshold = 40
print ("nFeature_RNA threshold = " %>% paste(nFeature_RNA_threshold))

[1] "nFeature_RNA threshold =  1000"

print ("percent.mt threshold = " %>% paste(percent.mt_threshold))

[1] "percent.mt threshold =  40"

acc_immune <- subset(acc_immune, subset = nFeature_RNA > nFeature_RNA_threshold & percent.mt < percent.mt_threshold)
acc_immune

An object of class Seurat 
22647 features across 621 samples within 1 assay 
Active assay: RNA (22647 features, 0 variable features)

6 PCA

# Identification of highly variable features
acc_immune <- FindVariableFeatures(acc_immune, selection.method = "vst", nfeatures = 15000) 

# Scaling the data
acc_immune <- ScaleData(acc_immune, vars.to.regress = c("percent.mt","nCount_RNA"))

# Perform linear dimensional reduction (PCA)
acc_immune <- RunPCA(acc_immune, features = VariableFeatures(object = acc_immune))

ElbowPlot(acc_immune, ndims = 50) # checking the dimensionality 

pc2use=1:10
clus_res=1
print("PCA dims = " %>% paste(max(pc2use)))

[1] "PCA dims =  10"

7 UMAP and clustering

acc_immune <- FindNeighbors(acc_immune, dims = pc2use,verbose = F)
acc_immune <- FindClusters(acc_immune, resolution = clus_res,verbose = F)

# Run non-linear dimensional reduction (UMAP)
acc_immune <- RunUMAP(acc_immune, dims = pc2use,verbose = F)
DimPlot(object = acc_immune, reduction = "umap", pt.size = 1, label = F)

8 UMAPS

#get metedata
plate = str_extract(colnames(acc_immune), "^.*-P[0-9]*")
patient.ident = str_extract(colnames(acc_immune), "ACC[0-9]*")
origin = str_extract(colnames(acc_immune), "LN")
origin = origin %>% replace_na('Primary')

acc_immune <- AddMetaData(object = acc_immune, metadata = as.factor(patient.ident), col.name = "patient.ident")
acc_immune <- AddMetaData(object = acc_immune, metadata = as.factor(plate), col.name = "plate")
acc_immune <- AddMetaData(object = acc_immune, metadata = as.factor(origin), col.name = "origin")

print_tab(DimPlot(acc_immune, reduction = "umap", label = F, pt.size = 1,group.by = "patient.ident") 
,title = "by patient")

by patient

print_tab(plt = 
            DimPlot(acc_immune, reduction = "umap", label = F, pt.size = 1,group.by = "plate") 
        ,title =  "by plate")

by plate

print_tab(plt =
            FeaturePlot(object = acc_immune, features = "PTPRC")
          ,
          title = "CD45")

CD45

print_tab(plt = 
            DimPlot(acc_immune, reduction = "umap", label = F, pt.size = 1,group.by = "origin") 
        ,title =  "by origin")

by origin

NA

library(SeuratDisk)
SaveH5Seurat(object = acc_immune,filename = "./Data/acc_immune")

9 Primary

acc_immune_pri  = subset(acc_immune, subset = origin == "Primary")

9.1 PCA

# Identification of highly variable features
acc_immune_pri <- FindVariableFeatures(acc_immune_pri, selection.method = "vst", nfeatures = 15000) 

# Scaling the data
acc_immune_pri <- ScaleData(acc_immune_pri, vars.to.regress = c("percent.mt","nCount_RNA"))

# Perform linear dimensional reduction (PCA)
acc_immune_pri <- RunPCA(acc_immune_pri, features = VariableFeatures(object = acc_immune_pri))

ElbowPlot(acc_immune_pri, ndims = 50) # checking the dimensionality 

pc2use=1:10
clus_res=1
print("PCA dims = " %>% paste(max(pc2use)))

[1] "PCA dims =  10"

9.2 UMAP and clustering

acc_immune_pri <- FindNeighbors(acc_immune_pri, dims = pc2use,verbose = F)
acc_immune_pri <- FindClusters(acc_immune_pri, resolution = clus_res,verbose = F)

# Run non-linear dimensional reduction (UMAP)
acc_immune_pri <- RunUMAP(acc_immune_pri, dims = pc2use,verbose = F,  metric = "euclidean")
DimPlot(object = acc_immune_pri, reduction = "umap", pt.size = 1, label = F)

9.3 Assignment

es.max = sctype_score(scRNAseqData = acc_immune_pri[["RNA"]]@scale.data, scaled = TRUE, 
                      gs = gs_list$gs_positive, gs2 = gs_list$gs_negative) 
# merge by cluster
cL_resutls = do.call("rbind", lapply(unique(acc_immune_pri@meta.data$seurat_clusters), function(cl){
    es.max.cl = sort(rowSums(es.max[ ,rownames(acc_immune_pri@meta.data[acc_immune_pri@meta.data$seurat_clusters==cl, ])]), decreasing = !0)
    head(data.frame(cluster = cl, type = names(es.max.cl), scores = es.max.cl, ncells = sum(acc_immune_pri@meta.data$seurat_clusters==cl)), 10)
}))
sctype_scores = cL_resutls %>% group_by(cluster) %>% top_n(n = 1, wt = scores)  

# set low-confident (low ScType score) clusters to "unknown"
sctype_scores$type[as.numeric(as.character(sctype_scores$scores)) < sctype_scores$ncells/4] = "Unknown"
print(sctype_scores[,1:3])

acc_immune_pri@meta.data$cell_identity = ""
for(j in unique(sctype_scores$cluster)){
  cl_type = sctype_scores[sctype_scores$cluster==j,]; 
  acc_immune_pri@meta.data$cell_identity[acc_immune_pri@meta.data$seurat_clusters == j] = as.character(cl_type$type[1])
}

DimPlot(acc_immune_pri, reduction = "umap", label = TRUE, repel = TRUE, group.by = 'cell_identity')  

cL_resutls

10 LN

acc_immune_ln  = subset(acc_immune, subset = origin == "LN")

10.1 PCA

# Identification of highly variable features
acc_immune_ln <- FindVariableFeatures(acc_immune_ln, selection.method = "vst", nfeatures = 15000) 

# Scaling the data
acc_immune_ln <- ScaleData(acc_immune_ln, vars.to.regress = c("percent.mt","nCount_RNA"))

# Perform linear dimensional reduction (PCA)
acc_immune_ln <- RunPCA(acc_immune_ln, features = VariableFeatures(object = acc_immune_ln))

ElbowPlot(acc_immune_ln, ndims = 50) # checking the dimensionality 

pc2use=1:10
clus_res=1
print("PCA dims = " %>% paste(max(pc2use)))

[1] "PCA dims =  10"

10.2 UMAP and clustering

acc_immune_ln <- FindNeighbors(acc_immune_ln, dims = pc2use,verbose = F)
acc_immune_ln <- FindClusters(acc_immune_ln, resolution = clus_res,verbose = F)

# Run non-linear dimensional reduction (UMAP)
acc_immune_ln <- RunUMAP(acc_immune_ln, dims = pc2use,verbose = F,  metric = "euclidean")
DimPlot(object = acc_immune_ln, reduction = "umap", pt.size = 1, label = F)

10.3 Assignment

es.max = sctype_score(scRNAseqData = acc_immune_ln[["RNA"]]@scale.data, scaled = TRUE, 
                      gs = gs_list$gs_positive, gs2 = gs_list$gs_negative) 
# merge by cluster
cL_resutls = do.call("rbind", lapply(unique(acc_immune_ln@meta.data$seurat_clusters), function(cl){
    es.max.cl = sort(rowSums(es.max[ ,rownames(acc_immune_ln@meta.data[acc_immune_ln@meta.data$seurat_clusters==cl, ])]), decreasing = !0)
    head(data.frame(cluster = cl, type = names(es.max.cl), scores = es.max.cl, ncells = sum(acc_immune_ln@meta.data$seurat_clusters==cl)), 10)
}))
sctype_scores = cL_resutls %>% group_by(cluster) %>% top_n(n = 1, wt = scores)  

# set low-confident (low ScType score) clusters to "unknown"
sctype_scores$type[as.numeric(as.character(sctype_scores$scores)) < sctype_scores$ncells/4] = "Unknown"
print(sctype_scores[,1:3])

acc_immune_ln@meta.data$cell_identity = ""
for(j in unique(sctype_scores$cluster)){
  cl_type = sctype_scores[sctype_scores$cluster==j,]; 
  acc_immune_ln@meta.data$cell_identity[acc_immune_ln@meta.data$seurat_clusters == j] = as.character(cl_type$type[1])
}

DimPlot(acc_immune_ln, reduction = "umap", label = TRUE, repel = TRUE, group.by = 'cell_identity')  

cL_resutls