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