SNL CMV Tumor Cells Control vs CXCR2i 010721
Abbie Ireland
1/7/2022
Load packages
suppressPackageStartupMessages({
library(BiocFileCache)
library(scran)
library(scater)
library(tidyverse)
library(SingleCellExperiment)
library(DropletUtils)
library(ggplot2)
library(cowplot)
library(Seurat)
library(monocle3)
library(SeuratObject)
})Read in all SNL seurat object from previous analysis
SNL_All<-readRDS("/scratch/Aireland/SNL for R01 Analysis/SNL_All tumor seurat object QC.rds")Subset just for CMV-Cre initiated samples
SNL_CMV<-subset(SNL_All, idents=c("SNL_CMV1","SNL_CMV2","SNL_CMV3","SNL_CMV_Tx"))
# Check cell number
table(SNL_CMV@meta.data$tx)##
## Control CXCR2i
## 8540 3008Normalize data now
SNL_CMV <- NormalizeData(SNL_CMV, normalization.method = "LogNormalize", scale.factor = 10000)
SNL_CMV## An object of class Seurat
## 32454 features across 11548 samples within 1 assay
## Active assay: RNA (32454 features, 0 variable features)Identify Highly variable features
SNL_CMV <- FindVariableFeatures(SNL_CMV, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
# features that are highly expressed in some cells, and lowly expressed in others.. often carrying biological
# significance in R
top10 <- head(VariableFeatures(SNL_CMV), 50)
top10## [1] "Scgb1a1" "Psca" "Scgb3a2" "Gkn2"
## [5] "Gkn1" "Ccl6" "Sprr2f" "Spp1"
## [9] "Reg3g" "Chil3" "Cnfn" "Wfdc17"
## [13] "Lyz1" "Ifit1bl1" "Ecm1" "Retnla"
## [17] "Krt15" "Scgb3a1" "Cyp2f2" "Clca3b"
## [21] "Tyrobp" "Cym" "Ambp" "Ctss"
## [25] "Ccdc153" "Cldn5" "Sprr2d" "Fam183b"
## [29] "Fcer1g" "Prap1" "Wfdc12" "Chil4"
## [33] "Duoxa2" "Ly6g6c" "Cxcl5" "9130204L05Rik"
## [37] "Ltf" "Prss27" "Rgs13" "Sprr1a"
## [41] "2300002M23Rik" "Reg3b" "Plet1" "Tuba1a"
## [45] "Ifit1" "AW112010" "Isg15" "Reg1"
## [49] "AU040972" "Tm4sf4"Scale data now
all.genes <- rownames(SNL_CMV)
SNL_CMV <- ScaleData(SNL_CMV, features = all.genes)## Centering and scaling data matrixPerform linear dimension reduction/ PCA
SNL_CMV <- RunPCA(SNL_CMV, features = VariableFeatures(object = SNL_CMV))## PC_ 1
## Positive: Sftpa1, Lyz2, Sftpb, Cd74, Slc34a2, H2-Aa, Cxcl15, Lamp3, Napsa, Scd1
## Hc, Sfta2, Lpcat1, Chil1, H2-Ab1, Ager, Lcn2, Tceal9, Apoe, Npc2
## Dram1, Cd36, H2-Eb1, Lgi3, Ppp1r14c, Tmem243, Selenow, Egfl6, Fasn, S100g
## Negative: S100a6, Gsto1, Plac8, Fxyd3, Lgals4, Krt19, Oit1, Wfdc1, Pglyrp1, Gsta4
## Krt18, Gnb2l1, Tff2, Pigr, Krt8, Agr2, Tff1, Anxa10, AA467197, Dpcr1
## Qsox1, Selk, Prom1, Sepp1, Atp5o, Ctse, Tmsb10, Tspan1, Foxq1, Il18
## PC_ 2
## Positive: Ramp2, Calcrl, Cldn5, Cdh5, Tmem100, Aqp1, Icam2, BC028528, Ptprb, Pecam1
## Thbd, Gpihbp1, Ace, Foxf1, Clec14a, Acvrl1, Clec1a, Ctla2a, S1pr1, Pltp
## Cav1, Hpgd, Cd93, Tspan7, Cavin2, Scn7a, Eng, Flt1, Cyyr1, Itga1
## Negative: Sftpd, Chchd10, Sfta2, Napsa, Slc34a2, Hc, Lamp3, Cxcl15, Chil1, Dram1
## Rnase4, Fasn, Npc2, Lgi3, Sftpb, Egfl6, Ppp1r14c, S100g, Scd1, Apoe
## Lcn2, Apoc1, Ptprf, Ager, Lpcat1, Emb, Cbr2, Nkx2-1, Ctsc, Lrp2
## PC_ 3
## Positive: Sftpd, Fabp5, Il33, Cd36, Ctsc, Lrg1, H2-Aa, Napsa, Hc, Fasn
## Slco2a1, H2-Eb1, Lamp3, Apoe, S100g, Cd200, Lpcat1, Lgi3, Slc34a2, Egfl7
## Nrp1, Sfta2, Chil1, Egfl6, Sparc, Aqp5, Selenop, Scd1, Pdzk1ip1, Etv5
## Negative: Ccdc153, Fam183b, 1110017D15Rik, Tmem212, 1700016K19Rik, 1700007K13Rik, Mlf1, Meig1, Efcab10, Dnah12
## Rsph1, Ccdc113, Sntn, Tctex1d4, 1700001C02Rik, Riiad1, Dynlrb2, Gm867, Foxj1, Cfap126
## Tekt1, Nme5, Stmnd1, Gm29538, Gm1673, Pifo, Sec14l3, 2410004P03Rik, BC051019, Iqca
## PC_ 4
## Positive: Ecm1, Clca3b, Prss27, Gm9573, Tmprss11g, Psca, Ifit1bl1, Nccrp1, Prss22, Irf7
## Fth1, Il1a, Duoxa2, Cnfn, 2300002M23Rik, Ly6a, Calml3, Arg1, Plpp1, Cysrt1
## Gpr87, Plat, Isg15, Ccng2, Mxd1, Pmaip1, Cdkn2b, Cpe, Pyhin1, Il1r2
## Negative: Car2, Tff2, Gstm1, Fabp5, Sepp1, Agr2, Gas6, Ctse, Lgals4, Gm3336
## Mgst2, Cfi, Sult1c2, Id1, Gsta4, Pla2g10, Tmsb10, Aqp5, Oit1, Prom1
## Pam, Vsig1, Id3, Ftl1, Csrp2, Sox9, Mgst3, Crip1, Cela1, Creb3l1
## PC_ 5
## Positive: Alox5ap, Ltc4s, Rgs13, Ly6g6f, Rac2, Dclk1, Sh2d6, Hck, Fyb, Hepacam2
## Gng13, Sh2d7, Trpm5, Matk, Cd300lf, Rgs2, Kctd12, Etv1, Pik3r5, Hmx2
## Vav1, Lrmp, Nrgn, Spib, Bmx, Hpgds, Espn, Cdhr2, Alox5, Avil
## Negative: Sftpd, Krt19, Gsta4, Il33, Mal, Krt7, S100a6, Wfdc1, Dpcr1, Oit1
## Pdzk1ip1, Fabp5, Lgals4, Phgr1, Areg, Fos, Gsto1, Agr2, Ctse, Rgs17
## Ccnd2, Phlda1, Sfn, Serpinb6b, 2210407C18Rik, Gprc5a, Krt20, Nupr1, Vsig1, Tacstd2DimPlot(SNL_CMV, reduction = "pca", group.by='tx')+ggtitle("PCA of CMV-Cre Initiated SNL Tumors")
Looks like we have some outlying cells that are likely non-tumor.
Determine the dimensionality for clustering
# Use elbow plot to visualize dimensions that encompass most of variance
ElbowPlot(SNL_CMV, ndims=50)
# shows ~25 dims captures most of variance
# Choosing 25 PCs for dimensionality
SNL_CMV <- FindNeighbors(SNL_CMV, dims = 1:25)## Computing nearest neighbor graph## Computing SNNSNL_CMV<- FindClusters(SNL_CMV, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 11548
## Number of edges: 398145
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9046
## Number of communities: 17
## Elapsed time: 1 secondsRun Non-Linear Dimension reduction (tSNE/UMAP)
SNL_CMV <- RunUMAP(SNL_CMV, dims = 1:25)## 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## 15:55:45 UMAP embedding parameters a = 0.9922 b = 1.112## 15:55:45 Read 11548 rows and found 25 numeric columns## 15:55:45 Using Annoy for neighbor search, n_neighbors = 30## 15:55:45 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 15:55:48 Writing NN index file to temp file /tmp/Rtmp3tK3hy/file97dd27f939aa
## 15:55:48 Searching Annoy index using 1 thread, search_k = 3000
## 15:55:52 Annoy recall = 100%
## 15:55:53 Commencing smooth kNN distance calibration using 1 thread
## 15:55:55 Initializing from normalized Laplacian + noise
## 15:55:57 Commencing optimization for 200 epochs, with 480098 positive edges
## 15:56:15 Optimization finishedPlot UMAP
x<-DimPlot(SNL_CMV, reduction = "umap",label=TRUE,pt.size=.1)+ggtitle("All cells captured by Cluster ID")
y<-DimPlot(SNL_CMV, reduction = "umap", group.by = "tx",label=FALSE,pt.size=0.1)+ggtitle("All cells captured by Treatment Status")
multiplot(x,y,cols=2)
There is a lot of variability. It does seem like treatment status is driving clustering, which is good..
Is clustering driven by batch effect?
DimPlot(SNL_CMV, reduction = "umap", group.by = "gnomex",label=FALSE,pt.size=0.1)+ggtitle("All cells captured by Experiment")
Batch effect is driving clustering too much. Going to just focus on 15242R experiment. See “CMV SNL tumor analysis 15739R” for analysis of additional CMV tumors separately, to see if captured more squamous cells. #####################################################################################
Subset and re-cluster
SNL_CMV<-subset(SNL_All, idents=c("SNL_CMV1","SNL_CMV_Tx"))
# Cells by treatment
table(SNL_CMV@meta.data$tx)##
## Control CXCR2i
## 3893 3008# Total cells
table(SNL_CMV@meta.data$Cre)##
## CMV
## 6901Normalize Data
SNL_CMV <- NormalizeData(SNL_CMV, normalization.method = "LogNormalize", scale.factor = 10000)
SNL_CMV## An object of class Seurat
## 32454 features across 6901 samples within 1 assay
## Active assay: RNA (32454 features, 0 variable features)Identify Highly variable features
SNL_CMV <- FindVariableFeatures(SNL_CMV, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
# features that are highly expressed in some cells, and lowly expressed in others.. often carrying biological
# significance in R
top10 <- head(VariableFeatures(SNL_CMV), 50)
top10## [1] "Lyz2" "Sftpc" "Sftpa1" "Ccl6"
## [5] "Sprr2f" "Sftpb" "Chil3" "Psca"
## [9] "Gkn2" "Wfdc17" "Gkn1" "Cnfn"
## [13] "Spp1" "Ifit1bl1" "Lcn2" "Reg3g"
## [17] "Tyrobp" "Ecm1" "Scd1" "Ctss"
## [21] "Sprr2d" "Krt15" "Clca3b" "Fcer1g"
## [25] "Cd74" "Chil1" "Prap1" "Apoe"
## [29] "Scgb1a1" "Hp" "9130204L05Rik" "Wfdc12"
## [33] "Cym" "Chil4" "Cxcl5" "Ifit1"
## [37] "Duoxa2" "Prss27" "Isg15" "Rgs13"
## [41] "Ambp" "2300002M23Rik" "S100g" "Scgb3a1"
## [45] "Ly6g6c" "Arg1" "Slc34a2" "Cxcl2"
## [49] "AW112010" "Ifit3b"Scale data
all.genes <- rownames(SNL_CMV)
SNL_CMV <- ScaleData(SNL_CMV, features = all.genes)## Centering and scaling data matrix##Perform linear dimension reduction/ PCA
SNL_CMV <- RunPCA(SNL_CMV, features = VariableFeatures(object = SNL_CMV))## PC_ 1
## Positive: Ecm1, Clca3b, Fth1, Prss27, Ly6a, Plpp1, Gm9573, Psca, Tmprss11g, Irf7
## Ifit1bl1, Il1a, Prss22, Nccrp1, Arg1, Isg15, Duoxa2, Pim1, Cnfn, 2300002M23Rik
## Calml3, Lgals3, Ccng2, Ier5, Plat, Mxd1, Sat1, Pmaip1, Gpr87, Abcg1
## Negative: Tff2, Wfdc2, Car2, Gsta4, Fabp5, Agr2, Gstm1, Krt8, Gm3336, Sult1c2
## Tff1, Gas6, Cfi, Pam, Cldn18, Id1, Sftpd, Aqp5, 2210407C18Rik, Dpcr1
## Krt18, AU021092, Id3, Prdx6, 5330417C22Rik, Dmbt1, Lrg1, Cbr2, Ptgr1, Sox9
## PC_ 2
## Positive: Fcer1g, Tyrobp, Ccl6, Wfdc17, Ctss, Cybb, Atp6v0d2, Vim, Mpeg1, Chil3
## Gpnmb, Trem2, Lyz2, Lilrb4a, Bcl2a1b, Plek, Ctsk, Cd52, Lgals1, Lpl
## Wfdc21, Fabp4, Cd84, Mcemp1, Mrc1, Sirpa, Cd53, Gngt2, Ptprc, Cd302
## Negative: Gsto1, Plac8, Krt7, Krt23, AA467197, Ltf, Clca3b, Psca, Ecm1, Prss27
## Ly6a, Gm9573, S100a14, Duoxa2, Ifit1bl1, Tmprss11g, Irf7, Plat, Mal, Rdh10
## Pmaip1, Tspan1, Nccrp1, Prss22, Sdcbp2, Plpp1, Cnfn, Pglyrp1, Lmo7, 2300002M23Rik
## PC_ 3
## Positive: Krt19, H2afz, Il33, Anxa10, Krt20, Dpcr1, Il18, Anxa2, Fabp5, Areg
## Ezr, Actg1, Ccnd2, Sftpd, Ftl1, Lmna, Wfdc2, Gapdh, Tnfrsf12a, Phgr1
## Msln, Nap1l1, Krt7, Stmn1, Ran, Birc5, Sfn, Tpm4, Serpinb6b, Far1
## Negative: Rgs13, Dclk1, Ly6g6f, Hepacam2, Sh2d6, Fyb, Hck, Gng13, Sh2d7, Trpm5
## Matk, Ltc4s, Hmx2, Etv1, Alox5ap, Nkd1, Lrmp, Rac2, Nrgn, Spib
## Espn, Selm, Kctd12, Cdhr2, Bmx, Itpr2, Agt, Hpgds, Vav1, 1810046K07Rik
## PC_ 4
## Positive: Tff2, Txnip, Cbr2, Gas6, 5330417C22Rik, Gstm1, Lrg1, Cfi, Sgk1, Aqp5
## Car2, Dmbt1, Pla2g1b, Fth1, AU021092, Clca3b, Igfbp5, Cela1, Prss27, Gm3336
## Chad, Pglyrp1, Ces1f, Gm9573, 2300002M23Rik, Fbxo32, Tmprss11g, Nccrp1, Klk1, 2010007H06Rik
## Negative: H2afz, Krt8, Gapdh, Anxa2, Rgs13, Dstn, Avil, Fyb, Dclk1, Ly6g6f
## Hepacam2, Il33, Ran, Sh2d6, Krt20, Actg1, Hck, Gng13, Krt18, Etv1
## Tnfrsf12a, Areg, Basp1, Sh2d7, Hsp90aa1, Ltc4s, Tm4sf1, Matk, Anxa3, Alox5ap
## PC_ 5
## Positive: Dynlrb2, Ccdc153, AU040972, 1110017D15Rik, Tmem212, Sntn, Fam183b, Vpreb3, Mlf1, Sec14l3
## Cfap126, 1700007K13Rik, Dnah12, Efcab10, 1700016K19Rik, Meig1, Ccdc113, 1700026L06Rik, Riiad1, Pcp4l1
## Tctex1d4, Hdc, Foxj1, Prr29, Gm867, 1700001C02Rik, Rsph1, Gm29538, Nme5, Cfap45
## Negative: Plac8, Fabp5, Tff1, AA467197, Dpcr1, Cldn18, Ltc4s, Alox5ap, Pglyrp1, Car2
## Rgs13, Sult1c2, Ly6g6f, Hck, Dclk1, Sh2d6, Fyb, Hepacam2, Avil, Phgr1
## Rac2, Gm3336, Sh2d7, Etv1, Matk, Trpm5, Cd300lf, Gng13, Pik3r5, Hmx2DimPlot(SNL_CMV, reduction = "pca", group.by='tx')+ggtitle("PCA of CMV-Cre Initiated SNL Tumors")
This looks better. Moving forward.
Determine dimensionality for clustering.
# Use elbow plot to visualize dimensions that encompass most of variance
ElbowPlot(SNL_CMV, ndims=50)
# shows ~25 dims captures most of variance
# Choosing 25 PCs for dimensionality
SNL_CMV <- FindNeighbors(SNL_CMV, dims = 1:25)## Computing nearest neighbor graph## Computing SNNSNL_CMV<- FindClusters(SNL_CMV, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 6901
## Number of edges: 236207
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8682
## Number of communities: 14
## Elapsed time: 0 secondsRun Non-Linear Dimension reduction (tSNE/UMAP)
SNL_CMV <- RunUMAP(SNL_CMV, dims = 1:25)## 15:56:54 UMAP embedding parameters a = 0.9922 b = 1.112## 15:56:54 Read 6901 rows and found 25 numeric columns## 15:56:54 Using Annoy for neighbor search, n_neighbors = 30## 15:56:54 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 15:56:56 Writing NN index file to temp file /tmp/Rtmp3tK3hy/file97dd7c192fe1
## 15:56:56 Searching Annoy index using 1 thread, search_k = 3000
## 15:56:58 Annoy recall = 100%
## 15:56:59 Commencing smooth kNN distance calibration using 1 thread
## 15:57:01 Initializing from normalized Laplacian + noise
## 15:57:01 Commencing optimization for 500 epochs, with 286310 positive edges
## 15:57:28 Optimization finishedPlot UMAP
x<-DimPlot(SNL_CMV, reduction = "umap",label=TRUE,pt.size=.1)+ggtitle("All cells captured by Cluster ID")
y<-DimPlot(SNL_CMV, reduction = "umap", group.by = "tx",label=FALSE,pt.size=0.1)+ggtitle("All cells captured by Treatment Status")
multiplot(x,y,cols=2)
Ensure we are analyzing only tumor cells
Assign cell types per cluster using Dotplots of common cell markers
more_types<-c("GFP","Sox2", "Nkx2-1", "Trp63", "Krt5","Krt13","Krt4",
"Myc","Yap1","Sox9","Runx2","Sox6","Sp7",
"Vcam1","Fgf10","Dcn","Clec3b","Fbln1", # rpma
"Col14a1", # matrix fibroblast/osteoclastic
"Acta2","Myh11","Tagln","Mustn1", #myofibroblast
"Lpl","Lipa","Pparg","Plin2","Ear1","Fabp1","Spp1", #lipofibroblast/osteoblastic
"Mertk","Marco","Mrc1","Ly75","Adgre1","Itgax","Cd68","Csf1r","Mafb","Msr1","Arg1","Adgre4","Clec4a1", #Macs/Myeloid
"Cx3cr1","Itgam","Cd14", #Monocytes
"S100a9","S100a8","Mmp9","Csf3r","Cxcr2","Ly6g", #Neuts
"Batf3","Xcr1","Clec9a","Ccl17","Ccl22", #DC
"Cd3d","Cd3e","Cd3g","Cd28","Cd8a","Cd4","Foxp3", # Tcell
"Gzma","Ncr1","Gzmb", #NK
"Fcmr","Cd19","Fcer2a","Pax5","Cd22","Cd79b","Cd79a", #B cells
"Slamf7", "Prdm1","Mzb1", #Plasma
"Mcam","Pecam1","Icam2","Cd36","Cd93", #Endothelial
"Scd1","Lamp3","Sftpb","Sftpa1","Aqp5", #Alveolar
"Scgb1a1","Cyp2f2","Scgb3a2","Scgb3a1","Lypd2", #club
"Foxj1","Rfx2","Rfx3","Trp73", #ciliated
"Muc5b","Muc5ac", #goblet
"Krt17","Perp", #basal
"Pou2f3","Ascl2","Rgs13","Il25","Tslp", #tuft
"Foxi1","Cftr","Ascl3", #ionocyte
"Epgn","Slc16a14","Grhl1", #Eosinophils
"Cpa3","Fcer1a","Ms4a2", #Basophils
"Cma1","Tpsb2", #Mast
"Cd200r3" #Combined with Mast, Baso, Eosinophil
)
DotPlot(SNL_CMV,features = more_types, group.by="seurat_clusters")+theme(axis.text.x=element_text(angle=90, hjust=1,size=8), axis.text.y=element_text(hjust=1,size=10))
Also perform unbiased marker analysis to confirm cell type ID
DEGs between treatments
Cluster.Markers <- FindAllMarkers(SNL_CMV, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)## Calculating cluster 0## Calculating cluster 1## Calculating cluster 2## Calculating cluster 3## Calculating cluster 4## Calculating cluster 5## Calculating cluster 6## Calculating cluster 7## Calculating cluster 8## Calculating cluster 9## Calculating cluster 10## Calculating cluster 11## Calculating cluster 12## Calculating cluster 13test<-as.data.frame(Cluster.Markers)
write.csv(test, "~/Cluster.Markers SNL CMV All cells captured.csv")
top100 <- Cluster.Markers %>% group_by(cluster) %>% top_n(n = 50, wt = avg_log2FC)
top100_names<-top100$gene
top100_names## [1] "Tff2" "Pigr" "mt-Co3" "Sepp1"
## [5] "mt-Nd1" "mt-Cytb" "Atp1b1" "mt-Nd4"
## [9] "Aqp5" "mt-Co1" "Fxyd3" "mt-Nd2"
## [13] "mt-Atp6" "mt-Co2" "Pde4d" "Car2"
## [17] "Cela1" "Xist" "Pam" "Ern2"
## [21] "Sgk1" "Dmbt1" "Sox9" "Galnt7"
## [25] "2010007H06Rik" "Ctse" "Mecom" "Cbr2"
## [29] "5330417C22Rik" "Stard10" "Tff1" "Txnip"
## [33] "AU021092" "Gm3336" "Tcf4" "C2cd4b"
## [37] "Lrg1" "PISD" "Thbs1" "Muc5b"
## [41] "Wnt4" "Pde2a" "Asph" "Galnt6"
## [45] "Itih2" "Creb3l4" "Pik3r1" "Id3"
## [49] "Irf2bp2" "Gkn3" "Fabp5" "Sptssb"
## [53] "Phgr1" "S100a6" "Dpcr1" "Tff1"
## [57] "Lgals4" "Oit1" "Tmsb10" "2210407C18Rik"
## [61] "Prss32" "Mgst3" "Pla2g10" "Gsta1"
## [65] "Aldh3a1" "Agr2" "Lypd8" "Sprr2a3"
## [69] "S100a1" "Bcas1" "Ucp2" "Glul"
## [73] "Muc16" "Gm3336" "Ces1f" "Fam132a"
## [77] "Gsta4" "Adh7" "Il18" "Sult1b1"
## [81] "Crip1" "Prr13" "Galnt6" "Id1"
## [85] "Spink4" "Dbi" "Gkn2" "Ctse"
## [89] "Mettl7b" "Prdx6" "Paqr5" "Far1"
## [93] "Ccnd2" "Rgs17" "Fut9" "Gcnt3"
## [97] "Rbp1" "Gkn1" "Rbp2" "Reg3g"
## [101] "Nupr1" "H2afz" "Ly6g6c" "Eif4ebp1"
## [105] "Pkib" "Upk3bl" "Tnfrsf12a" "Cxcl16"
## [109] "Fst" "Prss12" "Mthfd2" "Asns"
## [113] "Anxa3" "Dhrs1" "Anxa2" "Dstn"
## [117] "Il33" "Rhoc" "Tgif1" "Areg"
## [121] "Krt7" "Myl12a" "Serpinb9" "Fkbp1a"
## [125] "Serpinb6b" "Ffar4" "Pls3" "Tpm1"
## [129] "Krt23" "Anxa10" "Phlda1" "8430408G22Rik"
## [133] "Cd9" "Krt20" "Fam213a" "Krt18"
## [137] "Lgals3" "Nudt4" "Cd55" "Ugp2"
## [141] "F3" "Tm4sf1" "Actg1" "Azin1"
## [145] "Atf3" "S100g" "Hspb1" "Klf2"
## [149] "Anxa1" "AW112010" "Ttr" "Pla2g1b"
## [153] "Gas6" "Scgb3a1" "Gsta4" "Kcne3"
## [157] "Bex4" "Cyp4b1" "Lurap1l" "Sult1c2"
## [161] "Car8" "Gsta3" "Pdzk1ip1" "Sox4"
## [165] "Gstm1" "Rbp2" "Gsta1" "Clu"
## [169] "Mbip" "B4galnt1" "Ihh" "2210407C18Rik"
## [173] "Tmsb4x" "Cyb5a" "Fam3c" "Cfi"
## [177] "Ifitm3" "Thbs1" "Ly6d" "Car2"
## [181] "Sepp1" "Spink4" "Jun" "Smim5"
## [185] "Itm2c" "Lrpap1" "Tst" "Zfp36l1"
## [189] "Scara3" "Ptgr1" "Atp1a1" "Foxq1"
## [193] "Sptssb" "Hspa1a" "Kctd14" "Pgc"
## [197] "Ier2" "Egr1" "Socs3" "Junb"
## [201] "Ccnb2" "Cenpf" "Cdca3" "Cenpa"
## [205] "Birc5" "Knstrn" "Cdkn3" "Ube2c"
## [209] "Mki67" "Stmn1" "2810417H13Rik" "Pttg1"
## [213] "Lockd" "Spc24" "Hmgn2" "Cdca8"
## [217] "Ppia" "Ptma" "Ccdc34" "Rbm3"
## [221] "H2afz" "Hmgb1" "2700094K13Rik" "Ran"
## [225] "Slc25a5" "Txn1" "Gapdh" "Krt8"
## [229] "Actg1" "Hnrnpa1" "Atp5b" "Krt20"
## [233] "Eif5a" "Fabp5" "Hn1" "Hnrnpa3"
## [237] "Fdps" "Tubb5" "Khdc1a" "Nucks1"
## [241] "Lmna" "Lsm4" "Prdx4" "Vdac1"
## [245] "Hsp90b1" "Ptms" "Krt19" "Ptgr1"
## [249] "Arl6ip1" "Reg3g" "Hp" "Cyp2f2"
## [253] "Ifitm1" "Gp2" "Gabrp" "Aldh1a1"
## [257] "Cckar" "Cp" "Ifitm3" "Sox2"
## [261] "Hspb1" "Tmem176a" "Cbr2" "Wfdc2"
## [265] "Gpx8" "Scgb1a1" "Tmem176b" "Cldn3"
## [269] "Rdx" "Mgst1" "Wfdc3" "Chad"
## [273] "Epas1" "Ldhb" "Rbp4" "Bsg"
## [277] "Itm2b" "Krt4" "Scgb3a1" "Adh7"
## [281] "Crip2" "Them5" "Krt19" "Ifitm2"
## [285] "Cxcl17" "Reg1" "Ces1d" "Tmed3"
## [289] "Jun" "Sdpr" "Wfdc1" "Iyd"
## [293] "Reg3g" "Ier2" "Ly6d" "Btg2"
## [297] "Muc5b" "Fos" "Junb" "Atf3"
## [301] "Ecm1" "Ifit1bl1" "Cnfn" "Clca3b"
## [305] "Duoxa2" "Prss27" "Plpp1" "Irf7"
## [309] "Isg15" "2300002M23Rik" "Arg1" "Nccrp1"
## [313] "Ifit1" "Prss22" "Il1a" "Gm9573"
## [317] "Calml3" "Cysrt1" "Mxd1" "Ccng2"
## [321] "Tmprss11g" "Ly6a" "Sprr2f" "Ifit3"
## [325] "Plat" "Psca" "Pmaip1" "Ier5"
## [329] "S100a14" "Ltf" "Fth1" "Plet1"
## [333] "Emp1" "AA467197" "Tmbim4" "Gsto1"
## [337] "Malat1" "Map1lc3b" "Slc39a4" "Dusp1"
## [341] "Krt23" "Ctsb" "Rdh10" "Zfand5"
## [345] "Cited2" "Sat1" "Lgals3" "Glrx"
## [349] "Pglyrp1" "Hist1h1c" "Socs3" "Fosb"
## [353] "Egr1" "Jun" "Fos" "Junb"
## [357] "Zfp36" "Ier2" "Gadd45g" "Btg2"
## [361] "Klk1" "Hspa1a" "Pgc" "Tff2"
## [365] "Nr4a1" "Jund" "Dusp1" "Atf3"
## [369] "Creb3l4" "AU021092" "Klf2" "Cfi"
## [373] "Aqp5" "Klf6" "Dnajb1" "Gkn3"
## [377] "Car2" "Ctse" "Clu" "Pdia3"
## [381] "Klf5" "Lrg1" "Hgfac" "Bace2"
## [385] "Pdia6" "Gas6" "Vsig2" "Thbs1"
## [389] "Cela1" "C2cd4b" "Sepp1" "Cyp2b10"
## [393] "Cbr2" "Pla2g1b" "Id3" "Pam"
## [397] "Gmds" "Dmbt1" "Sgk1" "Muc5b"
## [401] "Trpv6" "1700011H14Rik" "Prap1" "Sult1d1"
## [405] "Spp1" "Cited4" "2310043M15Rik" "Ces1d"
## [409] "Ppp2r3a" "Ces1f" "Prom1" "Rbp1"
## [413] "Paqr5" "Pdzk1ip1" "Cystm1" "Aplp2"
## [417] "Ndufb8" "Socs2" "Adora1" "Lgals4"
## [421] "Rpl30" "Pmm1" "Ascc1" "Epb41l2"
## [425] "Gpx2" "Oit1" "Rpl32" "Smim6"
## [429] "Smim1" "Slc35a3" "Gm3336" "Pigr"
## [433] "Them5" "Fermt1" "Gclc" "Cd2ap"
## [437] "Noxo1" "Zfp703" "Ccnd2" "Igfbp5"
## [441] "Rbbp7" "Ugdh" "Vopp1" "Chd3"
## [445] "Cox17" "Ucp2" "Tmem54" "Sult1b1"
## [449] "C2cd4b" "Eps8" "Rgs13" "Dclk1"
## [453] "Alox5ap" "Ltc4s" "Fyb" "Gng13"
## [457] "Ly6g6f" "Etv1" "Hck" "Rgs2"
## [461] "Sh2d6" "Selm" "Rac2" "Lrmp"
## [465] "Hepacam2" "Avil" "Matk" "Trpm5"
## [469] "Pate4" "Kctd12" "Nrgn" "Sh2d7"
## [473] "Nkd1" "Itpr2" "Hmx2" "Hpgds"
## [477] "Bmx" "Dpysl2" "Strip2" "Spib"
## [481] "Zfp428" "Tuba1a" "Espn" "Ahnak2"
## [485] "Ptpn6" "Bcl2l14" "Fam221a" "Marcksl1"
## [489] "Bpgm" "Adh1" "Ptpn18" "Aldh2"
## [493] "Tspan6" "Myo1b" "Ethe1" "Skap2"
## [497] "Reep5" "Cnn3" "Calm2" "Crip1"
## [501] "2810417H13Rik" "Top2a" "Hist1h2ap" "Prc1"
## [505] "Nusap1" "Cdk1" "Spc25" "Ccna2"
## [509] "Tpx2" "Rrm2" "Pbk" "Esco2"
## [513] "Mki67" "Smc2" "Birc5" "Ccnb1"
## [517] "Cdc20" "Cdca8" "Cenpe" "Spc24"
## [521] "Ube2c" "Lockd" "Cenpa" "Cks1b"
## [525] "Cdca3" "Cks2" "Cenpf" "Knstrn"
## [529] "Lig1" "Stmn1" "Dut" "H2afx"
## [533] "Smc4" "Hmgb2" "Hmgn2" "Ccdc34"
## [537] "Tuba1b" "Tmpo" "Ptma" "Tubb5"
## [541] "2700094K13Rik" "Ube2s" "Ran" "Nucks1"
## [545] "Dek" "Ranbp1" "Hn1" "H2afv"
## [549] "Arl6ip1" "Gkn2" "Chil3" "Ccl6"
## [553] "Wfdc17" "Ctss" "Tyrobp" "Fcer1g"
## [557] "Wfdc21" "S100a8" "Cybb" "Ctsk"
## [561] "Fabp4" "Cxcl2" "Lpl" "Vim"
## [565] "Gpnmb" "Atp6v0d2" "Mpeg1" "Fxyd5"
## [569] "Car4" "Bcl2a1b" "Cd52" "Trem2"
## [573] "Fn1" "Gngt2" "Il1b" "Cd53"
## [577] "Mmp12" "Hebp1" "Mrc1" "Ear2"
## [581] "C1qb" "Lgals1" "Cd68" "S100a9"
## [585] "Lyz2" "Slpi" "Cd44" "Apoe"
## [589] "Plin2" "Spp1" "Laptm5" "Lgmn"
## [593] "Msrb1" "Ctsd" "Psap" "Ftl1"
## [597] "Ctsz" "B2m" "Ctsb" "Cd74"
## [601] "Sftpb" "Slc34a2" "Napsa" "Apoc1"
## [605] "Hc" "H2-Aa" "Ager" "Lamp3"
## [609] "Itih4" "Dram1" "Egfl6" "Fgg"
## [613] "Rgcc" "Cd36" "Lrp2" "Lgi3"
## [617] "Ppp1r14c" "Apoe" "Sftpa1" "Lpcat1"
## [621] "H2-Eb1" "H2-Ab1" "Cd302" "Lcn2"
## [625] "Sfta2" "Cxcl15" "Scd1" "Chil1"
## [629] "Cd74" "Lyz2" "Fasn" "S100g"
## [633] "Hp" "Atp11a" "Acsl4" "Sftpc"
## [637] "Prnp" "Lyz1" "Selenbp1" "Elovl1"
## [641] "Npc2" "Sftpd" "Rnase4" "Icam1"
## [645] "Ctsc" "Scd2" "Chchd10" "Abca3"
## [649] "Lrg1" "Arl6ip1" "AU040972" "Ccdc153"
## [653] "Tmem212" "Fam183b" "Dynlrb2" "Sec14l3"
## [657] "1110017D15Rik" "1700016K19Rik" "Pcp4l1" "Mlf1"
## [661] "Vpreb3" "1700007K13Rik" "Sntn" "Efcab10"
## [665] "Aoc1" "Meig1" "Dnah12" "Foxj1"
## [669] "Ccdc113" "Gm867" "Cdhr4" "Tekt1"
## [673] "Stmnd1" "Dnah5" "Cfap126" "Hdc"
## [677] "Rsph1" "Slc23a2" "Riiad1" "Tuba1a"
## [681] "BC048546" "Aldh3b1" "Traf3ip1" "Lrrc51"
## [685] "Ccdc181" "Ifitm1" "Aldh1a1" "Cyp2f2"
## [689] "Tppp3" "Acot1" "Hspa4l" "Gsn"
## [693] "1110004E09Rik" "Chchd10" "4933434E20Rik" "Igfbp5"
## [697] "Hsp90aa1" "Elof1" "Scgb1a1" "Tubb4b"DoHeatmap(subset(SNL_CMV, downsample = 100), features = top100_names, size = 3)
#group.colors=c("orchid4","orange"))Top 50 marker genes by treatment status
DT::datatable(
Cluster.Markers %>%
mutate(gene=rownames(.)) %>%
group_by(cluster) %>%
top_n(n = 50, wt = avg_log2FC) %>%
select(gene, everything())
)Cluster 11 is Macs
Cluster 12 is normal AT2/Club lung cells
Cluster 13 is ciliated cells
Cluster 9 are tuft cells! :)
Cluster 6 is macs/myeloid, but has tumor markers, tricky but keep
Now subset Seurat Object to just take ID’D tumor cells
SNL_CMV_tumor<-subset(SNL_CMV, idents=c("0","1","2","3","4","5","6","7", "8","10"))
# New cell number
table(SNL_CMV_tumor@meta.data$tx)##
## Control CXCR2i
## 3723 2862Re-clustering (PCA, UMAP)
all.genes <- rownames(SNL_CMV_tumor)
SNL_CMV_tumor <- ScaleData(SNL_CMV_tumor, features = all.genes)## Centering and scaling data matrix#Perform linear dimension reduction/ PCA
SNL_CMV_tumor <- RunPCA(SNL_CMV_tumor, features = VariableFeatures(object = SNL_CMV_tumor))## PC_ 1
## Positive: Ecm1, Clca3b, Prss27, Ly6a, Plpp1, Psca, Gm9573, Tmprss11g, Irf7, Fth1
## Ifit1bl1, Prss22, Nccrp1, Il1a, Duoxa2, Isg15, Arg1, 2300002M23Rik, Cnfn, Calml3
## Ltf, Ccng2, Plat, Krt23, Pmaip1, Mxd1, Gpr87, Cysrt1, Cdkn2b, Ier5
## Negative: Tff2, Fabp5, Wfdc2, Gstm1, Car2, Gsta4, Agr2, Ftl1, Gm3336, Sult1c2
## Gas6, Cfi, Tff1, Id1, Pam, Krt8, Aqp5, Lrg1, 2210407C18Rik, Sftpd
## Crip1, AU021092, Id3, Mgst1, 5330417C22Rik, Ucp2, Cldn18, Dpcr1, Cbr2, Hsp90b1
## PC_ 2
## Positive: Tff2, Txnip, Cbr2, 5330417C22Rik, Gas6, Fth1, Lrg1, Clca3b, Cfi, Gstm1
## Prss27, Sgk1, Gm9573, 2300002M23Rik, Nccrp1, Pla2g1b, Tmprss11g, Ecm1, Dmbt1, Aqp5
## Trp53inp2, Igfbp5, Car2, Irf7, Sox9, Isg15, Calml3, Chad, Arg1, Creb3l4
## Negative: H2afz, Anxa2, Il33, Krt8, Gapdh, Krt20, Actg1, Areg, Tnfrsf12a, Dstn
## Anxa10, Ran, Anxa3, Nap1l1, Serpinb6b, Tpm4, Tm4sf1, Eif4ebp1, Serpinb9, Nupr1
## Dhrs1, Lmna, Ezr, Stmn1, Pkib, Birc5, Hsp90aa1, Cxcl16, Gldc, Msln
## PC_ 3
## Positive: Rhoc, Nupr1, Krt18, Asns, Ffar4, Eif4ebp1, Atf3, Dstn, Areg, 8430408G22Rik
## Jund, Hspb1, Pkib, Mthfd2, Klf2, Tgif1, Prss12, Anxa3, Ftl1, F3
## Cxcl16, Plaur, Ly6g6c, Tnfrsf12a, Fst, Junb, Tm4sf1, Serpinb6b, Phgdh, Dhrs1
## Negative: Birc5, 2810417H13Rik, Cenpa, Mki67, Tpx2, Ccna2, Ube2c, Top2a, Nusap1, Cdc20
## Cenpf, Pbk, Ccnb1, Prc1, Cdca3, Cdca8, Cenpe, Hmmr, Racgap1, Ckap2l
## Cdkn3, Cks2, Ckap2, Casc5, Spc24, Cdk1, Spc25, Lockd, Ccnb2, Plk1
## PC_ 4
## Positive: Egfl7, Sparc, Ramp2, Aqp1, Epas1, Calcrl, Ly6c1, Cldn5, Vim, Ctla2a
## Clec14a, Cav1, BC028528, Ace, Cd36, Gng11, Nrp1, Pecam1, Igfbp7, Slc43a3
## Icam2, Ptprb, Tspan7, Cd200, Acvrl1, Tmem100, Thbd, Timp3, Kdr, Clec1a
## Negative: Plac8, Dpcr1, Tff1, Fabp5, Phgr1, Anxa10, Gsto1, Cldn18, Far1, AA467197
## Il33, Il18, Krt20, Sptssb, 2210407C18Rik, Azin1, Basp1, Gkn1, Krt8, Lypd8
## Gkn2, Mal, Ezr, Agr2, Sdcbp2, Gsta1, Neat1, Ces1f, Gcnt3, Mettl7b
## PC_ 5
## Positive: Hp, Cyp2f2, Cp, Jun, Gabrp, Tmem176b, Scgb3a1, Mgst1, Ifitm1, Gp2
## Cckar, Ier2, Fosb, Tmem176a, Junb, Egr1, Wfdc2, Fos, Chad, Zfp36
## Ldhb, Gpx8, Socs3, Cbr2, Hspb1, Atf3, Cldn3, Sftpd, Gsta3, Ifitm3
## Negative: Tff1, Dpcr1, Calcrl, Ramp2, Aqp1, Cldn5, Il33, BC028528, Clec14a, Egfl7
## Ly6c1, Thbd, Icam2, Vim, Pecam1, Far1, Hpgd, Ptprb, Ace, Ctla2a
## Phgr1, Sparc, Tmem100, Cd36, Stmn2, Clec1a, Kdr, Cav1, Gpihbp1, Gng11DimPlot(SNL_CMV_tumor, reduction = "pca", group.by='tx')+ggtitle("PCA of Tumor Cells Initated by CMV-Cre")
# Determine the dimensionality for clustering
# Use elbow plot to visualize dimensions that encompass most of variance
ElbowPlot(SNL_CMV_tumor, ndims=50)
# shows ~20 dims captures most of variance
# Choosing 15 PCs for dimensionality
SNL_CMV_tumor <- FindNeighbors(SNL_CMV_tumor, dims = 1:15)## Computing nearest neighbor graph## Computing SNNSNL_CMV_tumor<- FindClusters(SNL_CMV_tumor, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 6585
## Number of edges: 218933
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8574
## Number of communities: 10
## Elapsed time: 0 seconds# Run Non-Linear Dimension reduction (tSNE/UMAP)
SNL_CMV_tumor <- RunUMAP(SNL_CMV_tumor, dims = 1:15)## 16:01:24 UMAP embedding parameters a = 0.9922 b = 1.112## 16:01:24 Read 6585 rows and found 15 numeric columns## 16:01:24 Using Annoy for neighbor search, n_neighbors = 30## 16:01:24 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 16:01:25 Writing NN index file to temp file /tmp/Rtmp3tK3hy/file97dd60f146a
## 16:01:25 Searching Annoy index using 1 thread, search_k = 3000
## 16:01:28 Annoy recall = 100%
## 16:01:29 Commencing smooth kNN distance calibration using 1 thread
## 16:01:31 Initializing from normalized Laplacian + noise
## 16:01:31 Commencing optimization for 500 epochs, with 271004 positive edges
## 16:01:56 Optimization finished# Plot UMAP
x<-DimPlot(SNL_CMV_tumor, reduction = "umap",label=TRUE,pt.size=.1)+ggtitle("Tumor cells by cluster ID")
y<-DimPlot(SNL_CMV_tumor, reduction = "umap", group.by = "tx",label=FALSE,pt.size=0.1)+ggtitle("Tumor cells by treatment")
multiplot(x,y,cols=2)
CMV treated tumors don’t respond very much to neut depletion hmm…. There are still outyling pops, check cell types again.
DotPlot(SNL_CMV_tumor,features = more_types, group.by="seurat_clusters")+ theme(axis.text.x=element_text(angle=90, hjust=1,size=8), axis.text.y=element_text(hjust=1,size=10))
Cluster 9 is endothelial cells, re-subset and cluster
SNL_CMV_tumor<-subset(SNL_CMV, idents=c("0","1","2","3","4","5","6","7","8"))
table(SNL_CMV_tumor@meta.data$tx)##
## Control CXCR2i
## 3658 2806Re-clustering, PCA, UMAP
all.genes <- rownames(SNL_CMV_tumor)
SNL_CMV_tumor <- ScaleData(SNL_CMV_tumor, features = all.genes)## Centering and scaling data matrix#Perform linear dimension reduction/ PCA
SNL_CMV_tumor <- RunPCA(SNL_CMV_tumor, features = VariableFeatures(object = SNL_CMV_tumor))## PC_ 1
## Positive: Ecm1, Clca3b, Prss27, Ly6a, Plpp1, Psca, Gm9573, Tmprss11g, Irf7, Fth1
## Ifit1bl1, Prss22, Nccrp1, Duoxa2, Il1a, Isg15, Cnfn, Arg1, 2300002M23Rik, Calml3
## Ltf, Ccng2, Plat, Krt23, Pmaip1, Mxd1, Cysrt1, Gpr87, Cdkn2b, Ier5
## Negative: Tff2, Fabp5, Wfdc2, Gstm1, Car2, Gsta4, Agr2, Ftl1, Gm3336, Gas6
## Sult1c2, Cfi, Tff1, Pam, Id1, Aqp5, Lrg1, Krt8, 2210407C18Rik, AU021092
## Sftpd, Mgst1, 5330417C22Rik, Id3, Cbr2, Crip1, Cldn18, Dpcr1, Dmbt1, Ucp2
## PC_ 2
## Positive: Tff2, Txnip, Fth1, Clca3b, Prss27, Gm9573, Gas6, Nccrp1, Cbr2, 2300002M23Rik
## Ecm1, Tmprss11g, 5330417C22Rik, Lrg1, Gstm1, Car2, Sgk1, Cfi, Irf7, Trp53inp2
## Isg15, Calml3, Arg1, Il1r2, Il1a, Cnfn, Gpr87, Dmbt1, Spns2, Psca
## Negative: H2afz, Krt8, Areg, Anxa2, Il33, Dstn, Tnfrsf12a, Gapdh, Krt20, Anxa3
## Anxa10, Nupr1, Actg1, Eif4ebp1, Krt18, Tm4sf1, Serpinb6b, Serpinb9, Ffar4, Dhrs1
## Pkib, Rhoc, Cxcl16, Ezr, Nap1l1, Mthfd2, Basp1, Msln, Tpm4, Ftl1
## PC_ 3
## Positive: Egfl7, Sparc, Ramp2, Aqp1, Epas1, Calcrl, Ly6c1, Cldn5, Vim, Clec14a
## Ctla2a, Cav1, BC028528, Ace, Gng11, Nrp1, Cd36, Pecam1, Slc43a3, Igfbp7
## Icam2, Ptprb, Cd200, Acvrl1, Thbd, Tspan7, Tmem100, Kdr, Timp3, Clec1a
## Negative: Plac8, Dpcr1, Tff1, Fabp5, Gsto1, AA467197, Phgr1, Cldn18, Sptssb, Anxa10
## Far1, Il33, Krt20, Mal, 2210407C18Rik, Gsta1, Gkn1, Il18, Agr2, Lypd8
## Gkn2, Pglyrp1, Gm3336, Krt8, Glrx, Aldh3a1, Sult1b1, Gcnt3, Ces1f, Sprr2a3
## PC_ 4
## Positive: Tff1, Dpcr1, Ramp2, Calcrl, Aqp1, Cldn5, Egfl7, Phgr1, BC028528, Clec14a
## Hpgd, Fabp5, Ly6c1, Vim, Icam2, Pecam1, Ptprb, Ace, Ctla2a, Sparc
## Thbd, Stmn2, Tmem100, Cd36, Clec1a, Kdr, Cav1, Gng11, Gpihbp1, Slc43a3
## Negative: Hp, Cyp2f2, Cp, Jun, Tmem176b, Mgst1, Scgb3a1, Wfdc2, Junb, Fosb
## Ifitm1, Ier2, Gabrp, Gp2, Cckar, Chad, Tmem176a, Fos, Egr1, Atf3
## Hspb1, Zfp36, Ldhb, Cbr2, Socs3, Gpx8, Jund, Cldn3, Klf6, Kcne3
## PC_ 5
## Positive: Mthfd2, Prss12, Fst, Asns, Pkib, Eif4ebp1, Azin1, Rhoc, Nupr1, Phgdh
## Kit, Ptrh1, Ly6g6c, Far1, Areg, Dhrs1, E130012A19Rik, Tgif1, Abcc4, Upk3bl
## Atf5, Basp1, Ptgs2, Ereg, Trib3, Gcnt1, Fxyd5, Rgs6, Pax8, Maff
## Negative: Birc5, Ccnb2, Cenpa, Cdca3, Stmn1, Cenpf, 2810417H13Rik, Knstrn, Cdkn3, Mki67
## Ube2c, Pttg1, Hmgn2, Cdca8, Hmmr, Cenpe, Lockd, Cdc20, Ccnb1, Tpx2
## Spc24, Racgap1, Cenpm, Ccna2, Ccdc34, Mad2l1, C330027C09Rik, Cep55, Cenpw, NrmDimPlot(SNL_CMV_tumor, reduction = "pca", group.by='tx')+ggtitle("PCA of Tumor Cells Initated by CMV-Cre")
# Determine the dimensionality for clustering
# Use elbow plot to visualize dimensions that encompass most of variance
ElbowPlot(SNL_CMV_tumor, ndims=50)
# Choosing 15 PCs for dimensionality
SNL_CMV_tumor <- FindNeighbors(SNL_CMV_tumor, dims = 1:15)## Computing nearest neighbor graph## Computing SNNSNL_CMV_tumor<- FindClusters(SNL_CMV_tumor, resolution = 0.5)## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 6464
## Number of edges: 215442
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8546
## Number of communities: 10
## Elapsed time: 0 seconds# Run Non-Linear Dimension reduction (tSNE/UMAP)
SNL_CMV_tumor <- RunUMAP(SNL_CMV_tumor, dims = 1:15)## 16:02:32 UMAP embedding parameters a = 0.9922 b = 1.112## 16:02:32 Read 6464 rows and found 15 numeric columns## 16:02:32 Using Annoy for neighbor search, n_neighbors = 30## 16:02:32 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 16:02:33 Writing NN index file to temp file /tmp/Rtmp3tK3hy/file97dd17e5b665
## 16:02:33 Searching Annoy index using 1 thread, search_k = 3000
## 16:02:35 Annoy recall = 100%
## 16:02:37 Commencing smooth kNN distance calibration using 1 thread
## 16:02:39 Initializing from normalized Laplacian + noise
## 16:02:39 Commencing optimization for 500 epochs, with 265092 positive edges
## 16:03:04 Optimization finished# Plot UMAP
x<-DimPlot(SNL_CMV_tumor, reduction = "umap",label=TRUE,pt.size=.1)+ggtitle("Tumor cells by cluster ID")
y<-DimPlot(SNL_CMV_tumor, reduction = "umap", group.by = "tx",label=FALSE,pt.size=0.1)+ggtitle("Tumor cells by treatment")
multiplot(x,y,cols=2)
Is cluster 6 tumor?
DotPlot(SNL_CMV_tumor,features = more_types, group.by="seurat_clusters")+ theme(axis.text.x=element_text(angle=90, hjust=1,size=8), axis.text.y=element_text(hjust=1,size=10))
No reason to exclude, proceeding.
Is clustering driven by cell cycle phase?
library(hciR)
library(hciRdata)
s.genes<-cc.genes$s.genes
g2m.genes<-cc.genes$g2m.genes
s.genes.list<-subset(mouse94, mouse94$human_homolog %in% s.genes)
s.genes.mouse<-(s.genes.list$gene_name)
s.genes.mouse## [1] "Cdc45" "Uhrf1" "Mcm2" "Slbp" "Mcm5" "Pola1"
## [7] "Gmnn" "Cdc6" "Rrm2" "Atad2" "Dscc1" "Mcm4"
## [13] "Chaf1b" "Rfc2" "Msh2" "Fen1" "Hells" "Prim1"
## [19] "Tyms" "Mcm6" "Wdr76" "Rad51" "Pcna" "Ccne2"
## [25] "Casp8ap2" "Usp1" "Nasp" "Rpa2" "Ung" "Rad51ap1"
## [31] "Blm" "Pold3" "Rrm1" "Gins2" "Tipin" "Brip1"
## [37] "Dtl" "Exo1" "Ubr7" "Clspn" "E2f8" "Cdca7"g2m.genes.list<-subset(mouse94, mouse94$human_homolog %in% g2m.genes)
g2m.genes.mouse<-(g2m.genes.list$gene_name)
g2m.genes.mouse## [1] "Ube2c" "Lbr" "Ctcf" "Cdc20" "Cbx5" "Kif11" "Anp32e"
## [8] "Birc5" "Cdk1" "Tmpo" "Hmmr" "Aurkb" "Top2a" "Gtse1"
## [15] "Rangap1" "Cdca3" "Ndc80" "Kif20b" "Cenpf" "Nek2" "Nuf2"
## [22] "Nusap1" "Bub1" "Tpx2" "Aurka" "Ect2" "Cks1b" "Kif2c"
## [29] "Cdca8" "Cenpa" "Mki67" "Ccnb2" "Kif23" "Smc4" "G2e3"
## [36] "Tubb4b" "Anln" "Tacc3" "Dlgap5" "Ckap2" "Ncapd2" "Ttk"
## [43] "Ckap5" "Cdc25c" "Hjurp" "Cenpe" "Ckap2l" "Cdca2" "Hmgb2"
## [50] "Cks2" "Psrc1" "Gas2l3"SNL_CMV_tumor <- CellCycleScoring(SNL_CMV_tumor, s.features = s.genes.mouse, g2m.features = g2m.genes.mouse)
DimPlot(SNL_CMV_tumor, reduction = "umap", group.by = "Phase", pt.size=.2)+ggtitle("Tumor cells by cell cycle phase")
Looks okay.
What % of cells in each cell cycle phase occupy each cluster?
proportions <- as.data.frame(100*prop.table(x = table(SNL_CMV_tumor$seurat_clusters, SNL_CMV_tumor$Phase), margin = 2))
colnames(proportions)<-c("Cluster", "Sample", "Frequency")
library(ggpubr)##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legendggbarplot(proportions, x="Sample", y="Frequency", fill = "Sample", color="black",
group = "Sample", ylab = "Frequency (percent)", xlab="",
# palette = sample_cols
)+
theme_bw()+
facet_wrap(facets = "Cluster", scales="free_y", ncol =4)+
theme(axis.text.y = element_text(size=12)) +
rotate_x_text(angle = 45)
Cluster 7 is highly proliferative.
What % of cells in each treatment type occupy each cluster?
proportions <- as.data.frame(100*prop.table(x = table(SNL_CMV_tumor$seurat_clusters, SNL_CMV_tumor$tx), margin = 2))
colnames(proportions)<-c("Cluster", "Sample", "Frequency")
library(ggpubr)
ggbarplot(proportions, x="Sample", y="Frequency", fill = "Sample", color="black",
group = "Sample", ylab = "Frequency (percent)", xlab="",
# palette = sample_cols
)+
theme_bw()+
facet_wrap(facets = "Cluster", scales="free_y", ncol =4)+
theme(axis.text.y = element_text(size=12)) +
rotate_x_text(angle = 45)
Seems like CMV-Cre initiated tumors don’t respond as much to neutrophil depletion (or neut depletion was just less effective in this experiment…?)
Differential gene expression analysis
DEGs between seurat clusters
Cluster.Markers <- FindAllMarkers(SNL_CMV_tumor, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)## Calculating cluster 0## Calculating cluster 1## Calculating cluster 2## Calculating cluster 3## Calculating cluster 4## Calculating cluster 5## Calculating cluster 6## Calculating cluster 7## Calculating cluster 8## Calculating cluster 9test<-as.data.frame(Cluster.Markers)
write.csv(test, "~/SNL CMV Tumor Cluster Marker Genes.csv")
top100 <- Cluster.Markers %>% group_by(cluster) %>% top_n(n = 50, wt = avg_log2FC)
top100_names<-top100$gene
top100_names## [1] "Sptssb" "Fabp5" "Phgr1" "S100a6"
## [5] "Dpcr1" "Tff1" "Lgals4" "Tmsb10"
## [9] "Mgst3" "2210407C18Rik" "Oit1" "Ucp2"
## [13] "Gsta1" "S100a1" "Prss32" "Aldh3a1"
## [17] "Sprr2a3" "Crip1" "Bcas1" "Pla2g10"
## [21] "Lypd8" "Muc16" "Agr2" "Sult1b1"
## [25] "Adh7" "mt-Nd1" "Glul" "Ces1f"
## [29] "Dbi" "Gm3336" "Laptm5" "Il18"
## [33] "Tmsb4x" "Basp1" "Gsta4" "Fam132a"
## [37] "Prdx6" "Galnt6" "Id1" "Paqr5"
## [41] "Fut9" "Spink4" "Mettl7b" "Rbp1"
## [45] "Far1" "Ccnd2" "Tppp3" "Gkn2"
## [49] "Reg3g" "Gkn1" "Tff2" "Pigr"
## [53] "mt-Co3" "mt-Nd1" "Sepp1" "mt-Cytb"
## [57] "mt-Nd4" "mt-Co1" "mt-Nd2" "Aqp5"
## [61] "mt-Atp6" "mt-Co2" "Car2" "Cela1"
## [65] "Pam" "Pde4d" "H2afj" "Sox9"
## [69] "Rgs6" "Galnt7" "Ern2" "Xist"
## [73] "Sgk1" "2010007H06Rik" "Dmbt1" "Ctse"
## [77] "Cbr2" "Mecom" "Klk1" "Lrg1"
## [81] "Tff1" "5330417C22Rik" "Tcf4" "Thbs1"
## [85] "Pde2a" "AU021092" "Gm3336" "Galnt6"
## [89] "C2cd4b" "Gkn3" "PISD" "Creb3l1"
## [93] "Muc5b" "Creb3l4" "Itih2" "Pik3r1"
## [97] "Asph" "Id3" "Spdef" "Wnt4"
## [101] "H2afz" "Nupr1" "Eif4ebp1" "Pkib"
## [105] "Upk3bl" "Fst" "Prss12" "Mthfd2"
## [109] "Ly6g6c" "Cxcl16" "Asns" "Tnfrsf12a"
## [113] "Anxa3" "Dhrs1" "Anxa2" "Dstn"
## [117] "Il33" "Tgif1" "Rhoc" "Myl12a"
## [121] "Fkbp1a" "Areg" "Krt7" "Serpinb9"
## [125] "Serpinb6b" "Tpm1" "Pls3" "Hsp90aa1"
## [129] "Krt23" "Anxa10" "Phlda1" "8430408G22Rik"
## [133] "Cd9" "Fam213a" "Krt18" "Krt20"
## [137] "Lgals3" "Nudt4" "Cd55" "Ugp2"
## [141] "Tm4sf1" "F3" "Actg1" "Azin1"
## [145] "Atf3" "S100g" "Hspb1" "Klf2"
## [149] "Anxa1" "AW112010" "Ttr" "Pla2g1b"
## [153] "Gas6" "Kcne3" "Scgb3a1" "Bex4"
## [157] "Gsta4" "Lurap1l" "Sult1c2" "Sox4"
## [161] "Cyp4b1" "Gsta3" "Gstm1" "Mbip"
## [165] "Car8" "Ptma" "Pdzk1ip1" "B4galnt1"
## [169] "Clu" "Spp1" "Rbp2" "Cyb5a"
## [173] "Tmsb4x" "Ihh" "Gsta1" "Car2"
## [177] "Thbs1" "Lrpap1" "Cfi" "Itm2c"
## [181] "Fam3c" "Sepp1" "Ifitm3" "2210407C18Rik"
## [185] "Tgm2" "Tst" "Smim5" "Jun"
## [189] "Pgc" "Scara3" "Zfp36l1" "Spink4"
## [193] "Ptgr1" "Idh2" "Cbr1" "Cd14"
## [197] "Ly6d" "Ier2" "Junb" "Socs3"
## [201] "Fosb" "Socs3" "Egr1" "Jun"
## [205] "Fos" "Zfp36" "Junb" "Ier2"
## [209] "Hspa1a" "Btg2" "Gadd45g" "Atf3"
## [213] "Jund" "Dusp1" "Nr4a1" "Klf2"
## [217] "Dnajb1" "Klf6" "Tff2" "Creb3l4"
## [221] "Klk1" "AU021092" "Pdia3" "Pgc"
## [225] "Ctse" "Hgfac" "Klf5" "Pdia6"
## [229] "Cfi" "Oit1" "Slc25a5" "Bace2"
## [233] "Car2" "Atp5b" "Rhob" "Vsig1"
## [237] "Ubc" "Ckmt1" "Vsig2" "Agr2"
## [241] "Aqp5" "Cela1" "Cyp2b10" "Gmds"
## [245] "Gkn3" "Lrg1" "C2cd4b" "Id3"
## [249] "Klf4" "Sfn" "Hp" "Cyp2f2"
## [253] "Ifitm1" "Gp2" "Gabrp" "Aldh1a1"
## [257] "Cckar" "Cp" "Sox2" "Wfdc2"
## [261] "Ifitm3" "Tmem176a" "Cbr2" "Gpx8"
## [265] "Hspb1" "Cldn3" "Scgb1a1" "Tmem176b"
## [269] "Mgst1" "Wfdc3" "Ldhb" "Lcn2"
## [273] "Chad" "Rbp4" "Krt4" "Bsg"
## [277] "Adh7" "Them5" "Scgb3a1" "Fxyd6"
## [281] "Krt19" "Cxcl17" "Crip2" "Anxa1"
## [285] "Fam3c" "Ces1d" "Reg1" "Tmed3"
## [289] "Wfdc1" "Iyd" "Reg3g" "Tst"
## [293] "Jun" "Ly6d" "Ier2" "Muc5b"
## [297] "Btg2" "Fos" "Junb" "Atf3"
## [301] "Ecm1" "Ifit1bl1" "Cnfn" "Clca3b"
## [305] "Duoxa2" "Plpp1" "Prss27" "Irf7"
## [309] "Isg15" "2300002M23Rik" "Arg1" "Nccrp1"
## [313] "Ifit1" "Prss22" "Il1a" "Gm9573"
## [317] "Calml3" "Cysrt1" "Mxd1" "Ccng2"
## [321] "Tmprss11g" "Ly6a" "Sprr2f" "Ifit3"
## [325] "Plat" "Psca" "Ier5" "Pmaip1"
## [329] "Fth1" "Ltf" "S100a14" "Plet1"
## [333] "AA467197" "Tmbim4" "Emp1" "Gsto1"
## [337] "Malat1" "Map1lc3b" "Ctsb" "Slc39a4"
## [341] "Dusp1" "Krt23" "Zfand5" "Rdh10"
## [345] "Sat1" "Cited2" "Lgals3" "Glrx"
## [349] "Pglyrp1" "Hist1h1c" "Stmn1" "Cenpf"
## [353] "Ccnb2" "Cdca3" "Cenpa" "Birc5"
## [357] "Pttg1" "2810417H13Rik" "Mki67" "Knstrn"
## [361] "Ube2c" "Lockd" "Spc24" "Cdkn3"
## [365] "Cenpe" "Cdca8" "Cdc20" "Hmgn2"
## [369] "Ptma" "Ccdc34" "Ppia" "Hmgb1"
## [373] "2700094K13Rik" "H2afz" "Ran" "Rbm3"
## [377] "Nucks1" "Lsm5" "Txn1" "Tubb5"
## [381] "Ptms" "Cks2" "Eif5a" "Lsm4"
## [385] "Hn1" "Dynll1" "Hnrnpa3" "Khdc1a"
## [389] "Gapdh" "Hnrnpa1" "Slc25a5" "Vdac1"
## [393] "Fdps" "Erh" "Hsp90b1" "Arl6ip1"
## [397] "Krt20" "Hmgb2" "Gkn2" "Gkn1"
## [401] "Trpv6" "Spp1" "1700011H14Rik" "Prap1"
## [405] "Sult1d1" "Cited4" "Ces1d" "Ces1f"
## [409] "2310043M15Rik" "Ppp2r3a" "Prom1" "Rbp1"
## [413] "Pdzk1ip1" "Cystm1" "Rps29" "Aplp2"
## [417] "Ndufb8" "Paqr5" "Rpl37a" "Lgals4"
## [421] "Socs2" "Ascc1" "Pmm1" "Rpl32"
## [425] "Epb41l2" "Xist" "Smim6" "Smim1"
## [429] "Slc35a3" "Pigr" "Gpx2" "Chd3"
## [433] "Them5" "Gclc" "Noxo1" "Igfbp5"
## [437] "Fermt1" "Cnn3" "Gm3336" "Cox17"
## [441] "Vopp1" "Zfp703" "Rbbp7" "C2cd4b"
## [445] "Ccnd2" "Sorl1" "Eps8" "Ugdh"
## [449] "Sult1b1" "Id1" "Cldn5" "Ramp2"
## [453] "Sparc" "Calcrl" "Tmem100" "Gng11"
## [457] "Cd36" "Aqp1" "Cav1" "Gpihbp1"
## [461] "Icam2" "Ace" "Clec14a" "Ptprb"
## [465] "Cd200" "Ctla2a" "Pecam1" "Vim"
## [469] "Car4" "Clec1a" "BC028528" "Cdh5"
## [473] "Kdr" "Cxcl12" "Scn7a" "Igfbp4"
## [477] "Egfl7" "Ly6c1" "Stmn2" "Thbd"
## [481] "Nrp1" "Gpx3" "Igfbp7" "H2-Aa"
## [485] "Tuba1a" "Hpgd" "Tspan7" "Timp3"
## [489] "Cd74" "Epas1" "Ecscr" "Klf2"
## [493] "Ehd4" "Adgrf5" "Slco2a1" "Sdpr"
## [497] "Sptbn1" "Ifitm3" "Slc9a3r2" "Hilpda"DoHeatmap(subset(SNL_CMV_tumor, downsample = 100), features = top100_names, size = 3) + NoLegend()
Top 50 marker genes per seurat cluster
DT::datatable(
Cluster.Markers %>%
mutate(gene=rownames(.)) %>%
group_by(cluster) %>%
top_n(n = 50, wt = avg_log2FC) %>%
select(gene, everything())
)DEGs between treatments
Idents(SNL_CMV_tumor)<-"tx"
Tx.Markers <- FindAllMarkers(SNL_CMV_tumor, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)## Calculating cluster Control## Calculating cluster CXCR2itest<-as.data.frame(Tx.Markers)
write.csv(test, "~/Tx.Markers SNL CMV Tumor cells.csv")
top100 <- Tx.Markers %>% group_by(cluster) %>% top_n(n = 100, wt = avg_log2FC)
top100_names<-top100$gene
top100_names## [1] "Tpt1" "Rps12" "Rpl12" "Rplp1"
## [5] "Rpl10" "Xist" "Rpl21" "Rpl5"
## [9] "Rpl8" "Rpl23" "Gnb2l1" "Rps20"
## [13] "Eef1a1" "Rps10" "Rps24" "Rps8"
## [17] "Rpl7" "Rps2" "Stard10" "Rpl30"
## [21] "Csrp2" "Rpl18" "Rpl11" "Rpsa"
## [25] "Rps27a" "Rpl28" "Gstm1" "Gkn1"
## [29] "Tff2" "Rpl7a" "Gkn2" "Atp1b1"
## [33] "Smim6" "Ern2" "P4hb" "Rpl10a"
## [37] "Gnas" "Spdef" "Eef1b2" "Gas6"
## [41] "5330417C22Rik" "mt-Nd2" "Sox9" "AU021092"
## [45] "Fabp5" "Txnip" "Rab3d" "Thbs1"
## [49] "Chad" "B2m" "Cbr2" "Creb3l1"
## [53] "Car2" "Ifi30" "Ptprf" "Dpcr1"
## [57] "Cldn18" "Lypd8" "Aqp5" "Igfbp5"
## [61] "Wnt4" "Wfdc2" "Pdgfa" "Ihh"
## [65] "Galnt6" "Pde4d" "Lrpap1" "Cyp2s1"
## [69] "mt-Nd3" "Nr4a1" "Spink4" "Tff1"
## [73] "Bace2" "C2cd4b" "Tmem176b" "Egr1"
## [77] "Atf3" "Tmem176a" "Fos" "Hspa1a"
## [81] "Tppp3" "Jun" "Muc5b" "AW112010"
## [85] "Ltf" "Ecm1" "Ly6a" "Psca"
## [89] "Sdcbp" "Malat1" "Ly6d" "Lgals3"
## [93] "2200002D01Rik" "Cd81" "Far1" "Slc16a11"
## [97] "S100a11" "Fam213a" "Anxa2" "Krt7"
## [101] "Sdcbp2" "Pmaip1" "Sumo1" "Ttc9"
## [105] "Vmp1" "Sat1" "Plac8" "S100a14"
## [109] "Krt23" "Phlda1" "Sept7" "Anxa1"
## [113] "Zfand6" "Capza2" "Cd9" "Atp6v1e1"
## [117] "Krt20" "Rdh10" "Ggh" "Psme2"
## [121] "9130008F23Rik" "Capg" "Chic2" "Fkbp1a"
## [125] "Urah" "Nupr1" "Zfand5" "Pls3"
## [129] "Plscr1" "Prdx5" "Map1lc3b" "Serpinb1a"
## [133] "H2afz" "Tmbim4" "Abhd2" "Cited2"
## [137] "Taldo1" "Hist1h2bc" "Hspa5" "Ugp2"
## [141] "Dhrs1" "Actb" "Plet1" "Cldn23"
## [145] "AA467197" "Glrx" "Hist1h1c" "Slc39a4"
## [149] "H1f0"DoHeatmap(subset(SNL_CMV_tumor, downsample = 2000), features = top100_names, size = 3)
#group.colors=c("orchid4","orange"))Top 100 marker genes by treatment status
DT::datatable(
Tx.Markers %>%
mutate(gene=rownames(.)) %>%
group_by(cluster) %>%
top_n(n = 100, wt = avg_log2FC) %>%
select(gene, everything())
)Save work for easy return later.
saveRDS(SNL_CMV_tumor,“/scratch/Aireland/SNL for R01 Analysis/SNL_CMV_Tumor seurat.rds”) SNL_CMV_tumor<-readRDS(“/scratch/Aireland/SNL for R01 Analysis/SNL_CMV_Tumor seurat.rds”)