plate_bias_fixing

1 Functions

source_from_github(repositoy = "HMSC_functions",version = "0.1.12",script_name = "functions.R")

ℹ SHA-1 hash of file is 2934dee5f6b9fee69635192f19b1bcc205e05b62

2 Data

acc1_cancer_cells = readRDS("./Data/acc1_cancer_cells_15KnCount_V3.RDS")
acc1_cancer_cells$plate = acc1_cancer_cells$orig.ident

3 Original UMAP

DimPlot(object = acc1_cancer_cells,pt.size = 2,group.by = "plate")

4 Seurat intergration

acc1_cancer_cells.list <- SplitObject(acc1_cancer_cells, split.by = "plate")

# normalize and identify variable features for each dataset independently
acc1_cancer_cells.list <- lapply(X = acc1_cancer_cells.list, FUN = function(x) {
    # x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = acc1_cancer_cells.list,nfeatures = 2000)
acc.anchors <- FindIntegrationAnchors(object.list = acc1_cancer_cells.list, anchor.features = features,k.filter = 50)

acc.combined <- IntegrateData(anchorset = acc.anchors,k.weight = 50)

Merging dataset 1 into 2
Extracting anchors for merged samples
Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Integrating data

DefaultAssay(acc.combined) <- "integrated"

acc.combined <- ScaleData(acc.combined, verbose = FALSE)
acc.combined <- RunPCA(acc.combined, npcs = 30, verbose = FALSE)
ElbowPlot(acc.combined)

acc.combined <- RunUMAP(acc.combined, reduction = "pca", dims = 1:10)
acc.combined <- FindNeighbors(acc.combined, reduction = "pca", dims = 1:10)
acc.combined <- FindClusters(acc.combined, resolution = 0.5)

4.1 UMAPS

DimPlot(acc.combined, reduction = "umap", group.by = "plate")

DimPlot(acc.combined, reduction = "umap")

4.2 clusters DEG

acc_deg <- FindMarkers(acc.combined, ident.1 = "0",ident.2 = "1",features = VariableFeatures(acc.combined),densify = T)

enrichment_analysis(acc_deg,background = VariableFeatures(acc.combined),fdr_Cutoff = 0.01,ident.1 = "0",ident.2 = "1",show_by = 1)

4.3 myo-lum score

original_myo_genes = c( "TP63", "TP73", "CAV1", "CDH3", "KRT5", "KRT14", "ACTA2", "TAGLN", "MYLK", "DKK3")
original_lum_genes = c("KIT", "EHF", "ELF5", "KRT7", "CLDN3", "CLDN4", "CD24", "LGALS3", "LCN2", "SLPI" )

calculate_score = function(dataset,myo_genes,lum_genes,lum_threshold =1 , myo_threshold = -1) {
  myoscore=FetchData(object =dataset,vars =  myo_genes,slot = "data") %>% rowMeans()
  lescore=FetchData(object =dataset,vars =  lum_genes,slot = "data") %>% rowMeans()
  correlation = cor(lescore,myoscore) %>% round(digits = 2)
  message("correlation of lum score and myo score:" %>% paste(correlation))
  




  dataset=AddMetaData(dataset,lescore-myoscore,"luminal_over_myo")
  print(
    FeaturePlot(object = dataset,features = "luminal_over_myo")
  )
  data = FetchData(object = dataset,vars = "luminal_over_myo")
  print(
    data %>% 
    ggplot(aes( x=luminal_over_myo)) + 
    geom_density() 
    )
  
lum_cells_num = subset(x = dataset,luminal_over_myo >(lum_threshold)) %>% ncol() /ncol(dataset)
myo_cells_num = subset(x = dataset,luminal_over_myo <(myo_threshold)) %>% ncol()/ncol(dataset)
df = data.frame(cell_type = c("myo_cells","lum_cells"),percentage = c(myo_cells_num,lum_cells_num))
ggplot(data=df, aes(x=cell_type, y=percentage)) +
  geom_bar(stat="identity") + ggtitle("ACC cell types")
}

calculate_score(dataset = acc.combined,myo_genes = original_myo_genes,lum_genes = original_lum_genes)

Warning: Could not find TP63 in the default search locations, found in RNA assay instead
Warning: Could not find TP73 in the default search locations, found in RNA assay instead
Warning: Could not find CDH3 in the default search locations, found in RNA assay instead
Warning: Could not find MYLK in the default search locations, found in RNA assay instead
Warning: Could not find KIT in the default search locations, found in RNA assay instead
Warning: Could not find EHF in the default search locations, found in RNA assay instead
Warning: Could not find ELF5 in the default search locations, found in RNA assay instead
Warning: Could not find CLDN3 in the default search locations, found in RNA assay instead
Warning: Could not find CD24 in the default search locations, found in RNA assay instead
correlation of lum score and myo score: -0.04

acc1_cancer_cells.list <- SplitObject(acc1_cancer_cells, split.by = "plate")

# normalize and identify variable features for each dataset independently
acc1_cancer_cells.list <- lapply(X = acc1_cancer_cells.list, FUN = function(x) {
    # x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = nrow(acc1_cancer_cells))
})

# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = acc1_cancer_cells.list,nfeatures = nrow(acc1_cancer_cells))
acc.anchors <- FindIntegrationAnchors(object.list = acc1_cancer_cells.list, anchor.features = features,k.filter = 50)

acc.combined <- IntegrateData(anchorset = acc.anchors,k.weight = 50)

Merging dataset 1 into 2
Extracting anchors for merged samples
Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Integrating data

DefaultAssay(acc.combined) <- "integrated"

4.4 HPV-MYB

HPV33_P3 = fread("./Data/HPV33_P3.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P3.df = HPV33_P3 %>% mutate(
  plate = gsub(x =HPV33_P3$plate, replacement = "",pattern = "_.*$") 
  %>% gsub(pattern = "-P",replacement = ".P") 
  %>% gsub(pattern = "-",replacement = "_",)
  )
HPV33_P3.df = HPV33_P3.df %>% dplyr::filter(HPV33_P3.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P3.df)  <- HPV33_P3.df$plate
HPV33_P3.df$plate = NULL


HPV33_P2 = fread("./Data/HPV33_P2.txt",col.names = c("plate","reads")) %>% as.data.frame()
HPV33_P2.df = HPV33_P2 %>% mutate(
  plate = gsub(x =HPV33_P2$plate, replacement = "",pattern = "_.*$") 
  %>% gsub(pattern = "plate2-",replacement = "plate2_",)
  %>% gsub(pattern = "-",replacement = "\\.",)
  )
HPV33_P2.df = HPV33_P2.df %>% dplyr::filter(HPV33_P2.df$plate %in% colnames(acc1_cancer_cells))
rownames(HPV33_P2.df)  <- HPV33_P2.df$plate
HPV33_P2.df$plate = NULL

HPV33 = rbind(HPV33_P3.df,HPV33_P2.df)
acc.combined = AddMetaData(object = acc.combined,metadata = HPV33,col.name = "HPV33.reads")

hpv33_positive = HPV33 %>% dplyr::mutate(hpv33_positive = case_when(reads >= 10 ~ "positive",
                                                                    reads < 10 ~ "negative")
)



hpv33_positive$reads = NULL
acc.combined = AddMetaData(object = acc.combined,metadata = hpv33_positive)

myb_vs_hpv = FetchData(object = acc.combined,vars = c("hpv33_positive","MYB"))
myb_vs_hpv $hpv33_positive = as.character(myb_vs_hpv $hpv33_positive )

ggboxplot(myb_vs_hpv, x = "hpv33_positive", y = "MYB",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative")))+ stat_summary(fun.data = function(x) data.frame(y=15, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(MYB)")

4.5 HPV-MYB per plate

plate_1 = subset(acc.combined,subset = plate == "ACC.plate2")
myb_vs_hpv = FetchData(object = plate_1,vars = c("hpv33_positive","MYB"))
myb_vs_hpv $hpv33_positive = as.character(myb_vs_hpv $hpv33_positive )

ggboxplot(myb_vs_hpv, x = "hpv33_positive", y = "MYB",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative")))+ stat_summary(fun.data = function(x) data.frame(y=15, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(MYB)")+ggtitle("ACC.plate2")

plate_2 = subset(acc.combined,subset = plate == "ACC1.P3")
myb_vs_hpv = FetchData(object = plate_2,vars = c("hpv33_positive","MYB"))
myb_vs_hpv $hpv33_positive = as.character(myb_vs_hpv $hpv33_positive )

ggboxplot(myb_vs_hpv, x = "hpv33_positive", y = "MYB",
          palette = "jco",
          add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("positive","negative")))+ stat_summary(fun.data = function(x) data.frame(y=15, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(MYB)")+ggtitle("ACC1.P3")