cc_analysis

1 Functions

2 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)
cc_cancer <- NormalizeData(cc_cancer)
cc_cancer <- FindVariableFeatures(cc_cancer,nfeatures = 15000)
cc_cancer <- ScaleData(cc_cancer,features = VariableFeatures(cc_cancer))

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")
data

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)

cc_cancer  %<>%  RunPCA( features = VariableFeatures(object = cc_cancer))

PC_ 1 
Positive:  CD52, CD3D, VIM, SRGN, CCL5, PTPRC, CD2, RGS1, CD69, GZMA 
       TRAC, HCST, IL36G, CD3E, CD7, GZMB, NKG7, ALOX5AP, CD53, ARHGDIB 
       GNLY, SAMSN1, CD37, LAPTM5, CXCR4, GMFG, TRBC1, IGKC, LCP1, TRBC2 
Negative:  PTMA, CD9, RPS8, ATP1B3, RPL39, RPS5, PRDX1, IFITM3, RPS2, RPL8 
       SLC25A5, RPS24, RPS4X, RPS18, NDUFA4, RPL35A, RPS6, RPL13, LY6D, RPS27A 
       RPS12, RPL18A, RPL11, S100A6, KRT19, RPLP0, RPL7A, RPL26, RPL19, H2AFZ 
PC_ 2 
Positive:  HSPE1, PSME2, LDHB, HMGN2, NME1, ATP5G1, HSPD1, BIRC5, MRPL12, PSMB9 
       CD74, MRPS12, PTTG1, HMGB1, VIM, CDKN3, SRGN, POLR2L, HLA-DRB1, TMEM160 
       SNRPF, CENPF, CD52, PHB, SRM, CDC20, RANBP1, CKS2, HSP90AB1, CKS1B 
Negative:  ERO1A, MXD1, RHCG, NDRG1, IL1RN, CSTB, EGLN3, ANXA1, F3, SLC2A1 
       VEGFA, IL36G, ECM1, CEACAM1, FAM83A, KRT17, SLC6A8, GJB2, ACKR3, SERPINB1 
       PERP, BNIP3, PLAU, BMP2, DSC2, PGF, EMP1, NDUFA4L2, CDKN1A, CEACAM5 
PC_ 3 
Positive:  SRGN, VIM, LAPTM5, RGS1, FCER1G, TYROBP, C1QB, C1QA, CD53, ARHGDIB 
       ALOX5AP, C1QC, PTPRC, MS4A6A, HCST, FCGR3A, AIF1, ITGB2, LCP1, HLA-DQA1 
       LYZ, TMEM176B, COTL1, FCGR2A, CD37, HAVCR2, FGL2, SPP1, SAMSN1, CD14 
Negative:  LY6D, KRT5, KRT13, LYNX1, SERPINB3, DSP, SERPINB13, KRT19, KRT6A, GPX2 
       PERP, DSG3, KLHDC7B, ATP1B3, CRABP2, RHOV, CTHRC1, SFN, C10orf99, KRT1 
       CENPW, AQP3, CLCA2, S100A8, ID1, SLPI, CSTA, ADH7, S100A2, GABRP 
PC_ 4 
Positive:  MMP1, CAV1, MMP7, IGFBP7, MMP10, TNFRSF12A, S100A2, FGFBP1, COL17A1, LAMB3 
       FLNA, ANGPTL4, S100A10, BCAM, MT1X, RPS3, RPL18A, RPS2, MT2A, KRT14 
       RPS4X, CRIP2, RPL7A, GAPDH, RPS18, CRIP1, EFEMP1, MT1E, RPL8, RPS5 
Negative:  SPRR1B, TMPRSS11E, TOP2A, SPRR3, NUSAP1, MKI67, TMPRSS11D, TXNIP, NCCRP1, ASPM 
       IVL, CEACAM7, DUOX2, ANLN, ARL6IP1, KRT1, SPRR2D, UBE2C, CENPF, SPINK5 
       TPX2, AURKB, CDK1, CEACAM6, DLGAP5, KLK11, CLCA4, CXCL17, SPRR2A, PRC1 
PC_ 5 
Positive:  ANLN, TOP2A, UBE2C, ASPM, CENPF, MKI67, CDC20, DLGAP5, BIRC5, NUSAP1 
       CDKN3, CCNB2, TPX2, CENPA, CEP55, CCNA2, AURKB, CCNB1, UBE2S, KIF23 
       PRC1, CKS2, DEPDC1, PLK1, GTSE1, PTTG1, KIF20B, NUF2, KPNA2, CENPE 
Negative:  SPRR1B, CSTA, S100A8, CALML3, S100A9, SEPP1, LGALS7, AQP3, KRT1, BAALC 
       KRT13, LYNX1, ADH7, C10orf99, RAB25, SPINK5, SERPINB3, LYPD3, DEFB1, CTHRC1 
       GLUL, KRT16, SERPINB4, AKR1C3, CXCL17, CRABP2, LGALS7B, TMEM45A, TXNIP, LAGE3 

ElbowPlot(cc_cancer)

cc_cancer  %<>%  FindNeighbors(dims = 1:10) %>%  FindClusters(resolution = 0.5) %>%  RunUMAP(dims = 1:10)

Computing nearest neighbor graph
Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 14220
Number of edges: 423754

Running Louvain algorithm...

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

Maximum modularity in 10 random starts: 0.8478
Number of communities: 8
Elapsed time: 2 seconds

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
10:34:15 UMAP embedding parameters a = 0.9922 b = 1.112
10:34:15 Read 14220 rows and found 10 numeric columns
10:34:15 Using Annoy for neighbor search, n_neighbors = 30
10:34:15 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:34:17 Writing NN index file to temp file /tmp/RtmpC9FeRs/file11da44eabbbf1
10:34:17 Searching Annoy index using 1 thread, search_k = 3000
10:34:24 Annoy recall = 99.96%
10:34:25 Commencing smooth kNN distance calibration using 1 thread
10:34:28 Initializing from normalized Laplacian + noise
10:34:29 Commencing optimization for 200 epochs, with 569862 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
10:34:37 Optimization finished

DimPlot(cc_cancer)

FeaturePlot(object = cc_cancer,features = "HPV_reads")

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


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

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

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_all = cbind(average_data1,average_data2) 
average_all = average_all[,c(-1)]
average_all = average_all %>% dplyr::filter(hpv_positive == "positive")

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'

notch_genes = c("JAG1","JAG2","NOTCH3","NOTCH2","NOTCH1","DLL1","MYB","HES4","HEY1","HEY2","NRARP")
for (gene  in notch_genes) {
  myb_vs_hpv = FetchData(object = cc_cancer,vars = c("hpv_positive",gene)) 
  myb_vs_hpv$hpv_positive = paste("HPV",myb_vs_hpv$hpv_positive)

  p = ggboxplot(myb_vs_hpv, x = "hpv_positive", y = gene,
            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)
}

##   JAG1 {.unnumbered }  

 

##   JAG2 {.unnumbered }  

 

##   NOTCH3 {.unnumbered }  

 

##   NOTCH2 {.unnumbered }  

 

##   NOTCH1 {.unnumbered }  

 

##   DLL1 {.unnumbered }  

 

##   MYB {.unnumbered }  

 

##   HES4 {.unnumbered }  

 

##   HEY1 {.unnumbered }  

 

##   HEY2 {.unnumbered }  

 

##   NRARP {.unnumbered }  

NA