cc_analysis

1 Functions

2 Data

3 CC data

cc_data  = Read10X(data.dir = "/sci/labs/yotamd/lab_share/avishai.wizel/R_projects/HMSC/Data/cervical_cancer_data/Tumor/")

cc_cancer = CreateSeuratObject(counts = cc_data, min.cells = 3, min.features = 200)
cc_cancer[["percent.mt"]] <- PercentageFeatureSet(cc_cancer, pattern = "^MT-")
cc_cancer <- subset(cc_cancer, subset = percent.mt < 10)

cc_cancer <- NormalizeData(cc_cancer)
cc_cancer <- FindVariableFeatures(cc_cancer,nfeatures = 2000)
cc_cancer <- ScaleData(cc_cancer,features = VariableFeatures(cc_cancer),vars.to.regress = c("percent.mt","nCount_RNA"))
cc_cancer  %<>%  RunPCA( features = VariableFeatures(object = cc_cancer))

ElbowPlot(cc_cancer,ndims = 50)

4 HPV

hpv_bam = read_tsv(file = "./Data/cervical_cancer_data/hpv_bam.tsv",col_select = 1,col_names = F)

Rows: 29481 Columns: 1
── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (1): X1

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

reads_per_cell = str_extract(pattern = "CR:Z:.*$",string = hpv_bam$X1) %>% str_sub( start = 6) %>% table() %>% as.data.frame() %>% column_to_rownames(".")
rownames(reads_per_cell) = paste0(rownames(reads_per_cell),"-1")
reads_per_cell

zero_cells = colnames(cc_cancer)[!colnames(cc_cancer) %in% rownames(reads_per_cell)]
zero_cells_freq = data.frame(row.names = zero_cells,Freq = rep(0,length(zero_cells)))
reads_per_cell_all_cells = rbind(reads_per_cell,zero_cells_freq)
cc_cancer %<>% AddMetaData(metadata = reads_per_cell_all_cells,col.name = "HPV_reads")
hpv_reads = FetchData(object = cc_cancer,vars = "HPV_reads")

data = hpv_reads %>% mutate("0 reads" = if_else(condition = HPV_reads == 0,true = 1,false = 0))
data = data %>% mutate("1 reads" = if_else(condition = HPV_reads == 1,true = 1,false = 0))
data = data %>% mutate("2 reads" = if_else(condition = HPV_reads == 2,true = 1,false = 0))
data = data %>% mutate("3-23 reads" = if_else(condition = HPV_reads >=3 &HPV_reads  <24,true = 1,false = 0))
data = data %>% mutate("24+ reads" = if_else(condition = HPV_reads >=24,true = 1,false = 0))
data = colSums(data[,2:ncol(data)]) %>% as.data.frame()
names(data)[1] = "count"
data = rownames_to_column(data,var = "bin")
ggplot(data=data, aes(x=factor(bin,levels = c("0 reads","1 reads","2 reads","3-23 reads","24+ reads")), y=count)) +
  geom_bar(stat="identity", fill="steelblue") + xlab("HPV Reads")+ theme_minimal()+
  geom_text(aes(label=count), vjust=-0.5, color="black", size=3.5)

5 Cell types

library(patchwork)

p1 = FeaturePlot(object = cc_cancer,features = c("CDH1","CDKN2A","EPCAM")) + plot_annotation(
  title = 'cancer markers')

p2 = FeaturePlot(object = cc_cancer,features = c("CD27", "PRF1","CD28"))+ plot_annotation(
  title = 'Lymphocytes markers')

p3 = FeaturePlot(object = cc_cancer,features = c("CD163", "FCGR2A"),ncol = 1)+ plot_annotation(
  title = 'Macrophages markers')

print_tab(plt = p1,title = "cancer markers")

cancer markers

print_tab(plt = p2,title = "Lymphocytes markers")

Lymphocytes markers

print_tab(plt = p3,title = "Macrophages markers")

Macrophages markers

NA

clusters 4 and 7 have been removed

cc_cancer <- subset(cc_cancer, subset = seurat_clusters != 4)
cc_cancer <- subset(cc_cancer, subset = seurat_clusters != 7)

6 MYB

FeaturePlot(object = cc_cancer,features = c("MYB"))

df = cc_cancer@assays$RNA@data["MYB",] %>% as.data.frame()
names(df) = "MYB"

data = df %>% mutate("MYB=0" = if_else(condition = MYB == 0,true = 1,false = 0))
data = data %>% mutate("0>MYB<1"  = if_else(condition = MYB >0 & MYB  < 1,true = 1,false = 0))
data = data %>% mutate("MYB >= 1"  = if_else(condition = MYB >=1,true = 1,false = 0))
data = colSums(data[,2:ncol(data)]) %>% as.data.frame()
names(data)[1] = "count"
data = rownames_to_column(data,var = "bin")
ggplot(data=data, aes(x=factor(bin,levels = c("MYB=0","0>MYB<1","MYB >= 1")), y=count)) +
  geom_bar(stat="identity", fill="steelblue") + xlab("log(TPM) levels")+ theme_minimal()+
  geom_text(aes(label=count), vjust=-0.5, color="black", size=3.5)

7 HPV positive

hpv_positive = hpv_reads %>% dplyr::mutate(hpv_positive = case_when(HPV_reads >= 1 ~ "positive",
                                                                    HPV_reads < 1 ~ "negative")
)


cc_cancer = AddMetaData(object = cc_cancer,metadata = hpv_positive)

DimPlot(object = cc_cancer,group.by  = c("hpv_positive"),pt.size = 0.5)

8 MYB positive

myb_data = FetchData(object = cc_cancer,vars = ("MYB"))%>% dplyr::mutate(myb_positive = case_when(MYB >= 0.2 ~ "positive",
                                                                    MYB < 0.2 ~ "negative")) %>% dplyr::select("myb_positive")



cc_cancer = AddMetaData(object = cc_cancer,metadata = myb_data,col.name = "MYB_positive")
DimPlot(object = cc_cancer,group.by  = c("MYB_positive"),pt.size = 0.5)

9 HPV percent / MYB levels

library(ggrepel)
data = FetchData(object = cc_cancer,vars =  c("MYB","hpv_positive","seurat_clusters"))
average_data1  = data %>% group_by(seurat_clusters) %>%
    dplyr::summarize(average_myb = mean(MYB, na.rm=TRUE))

average_data2  = data %>%  group_by(seurat_clusters,hpv_positive) %>%
  summarise(count = n(),) %>%  mutate(hpv_percent = count / sum(count))

`summarise()` has grouped output by 'seurat_clusters'. You can override using the `.groups` argument.

average_data2 = average_data2 %>% dplyr::filter(hpv_positive == "positive")


average_all = cbind(average_data1,average_data2) 
average_all = average_all[,c(-1)]

ggplot(average_all,
           aes(x = hpv_percent, y = average_myb, label=seurat_clusters)) +  geom_smooth(method = lm) +
  geom_point() + stat_cor(method = "pearson")+geom_text_repel()

`geom_smooth()` using formula 'y ~ x'

10 HPV percent / MYB percent

library(ggrepel)
data = FetchData(object = cc_cancer,vars =  c("MYB_positive","hpv_positive","seurat_clusters"))
data$MYB_positive = as.factor(data$MYB_positive) 
average_data1  = data %>%  group_by(seurat_clusters,MYB_positive, .drop = FALSE) %>%
  summarise(count = n())%>%  mutate(myb_percent = count / sum(count))

`summarise()` has grouped output by 'seurat_clusters'. You can override using the `.groups` argument.

average_data1 = average_data1 %>% dplyr::filter(MYB_positive == "positive" & seurat_clusters!=4 & seurat_clusters!=7 )


average_data2  = data %>%  group_by(seurat_clusters,hpv_positive) %>%
  summarise(count = n()) %>%  mutate(hpv_percent = count / sum(count))

`summarise()` has grouped output by 'seurat_clusters'. You can override using the `.groups` argument.

average_data2 = average_data2 %>% dplyr::filter(hpv_positive == "positive") %>% dplyr::filter(seurat_clusters!=4 & seurat_clusters!=7 )
average_data2 = average_data2[,c(-1)]


average_all = cbind(average_data1,average_data2) 

New names:
• `count` -> `count...3`
• `count` -> `count...6`

ggplot(average_all,
           aes(x = hpv_percent, y = myb_percent, label=seurat_clusters)) +  geom_smooth(method = lm) +
  geom_point() + stat_cor(method = "pearson")+geom_text_repel()

`geom_smooth()` using formula 'y ~ x'

11 Boxplot HPV_positive VS MYB

myb_vs_hpv = FetchData(object = cc_cancer, vars = c("hpv_positive", "MYB"))
myb_vs_hpv$hpv_positive = paste("HPV", myb_vs_hpv$hpv_positive)

p = ggboxplot(myb_vs_hpv, x = "hpv_positive", y = "MYB",
            palette = "jco",
            add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("HPV positive","HPV negative")))+ stat_summary(fun.data = function(x) data.frame(y=max(x)*1.2, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") +ylab("log2(gene)")+ggtitle(gene)
  print_tab(p,title = gene)

##   MYB {.unnumbered }  

NA

12 Boxplot MYB_positive VS HPV

myb_vs_hpv = FetchData(object = cc_cancer, vars = c("HPV_reads", "MYB_positive"))
myb_vs_hpv$MYB_positive = paste("MYB", myb_vs_hpv$MYB_positive)

p = ggboxplot(myb_vs_hpv, x = "MYB_positive", y = "HPV_reads",
            palette = "jco",
            add = "jitter")+ stat_compare_means(method = "wilcox.test",comparisons = list(c("MYB positive","MYB negative")))+ stat_summary(fun.data = function(x) data.frame(y=max(x)*1.2, label = paste("Mean=",round(mean(x),digits = 2))), geom="text") 
  print_tab(p,title = gene)

##   MYB {.unnumbered }  

NA

library(ggstatsplot)
df  = FetchData(object = cc_cancer,vars = c("hpv_positive","MYB_positive")) %>% droplevels() 
test = fisher.test(table(df))
plt = ggbarstats(
  df,
  hpv_positive,
  MYB_positive,
  results.subtitle = FALSE,
  subtitle = paste0("Fisher's exact test", ", p-value = ",
                    round(test$p.value, 13))
)

print_tab(plt = plt,title = "fisher")

##   fisher {.unnumbered }  

NA