Preprocess
#use all pc gene
all_pc_gene=read.table("/Users/niehu/Documents/IGDB/lab_works/SMY/genes/mm10.all_pc.gene_id2gene_name.txt", header=FALSE, row.names=1, stringsAsFactors=FALSE)
colnames(all_pc_gene)=c("GeneName")
all_pc_gene=all_pc_gene[!duplicated(all_pc_gene$GeneName),,drop=FALSE]
pos = rownames(data) %in% rownames(all_pc_gene)
data=data[pos,]
rownames(data)= all_pc_gene[rownames(data),1]
data=data[,data["Nes",] >10]
head(data[,c(1:10)])
## X0dpi.P0.A2 X0dpi.P0.B8 X0dpi.P0.B9 X0dpi.P0.C11 X0dpi.P0.C4
## mt-Nd6 712.7080 553.829 924.844 1012.370 794.0990
## mt-Nd4 9738.9700 3897.930 12857.200 8432.000 7849.8600
## mt-Nd3 78.6259 118.362 201.452 266.859 17.5377
## mt-Co3 4708.0000 3427.060 9657.460 3658.390 8860.1700
## mt-Atp8 2139.4300 1270.890 3218.900 2151.670 1998.0000
## mt-Co2 2822.1700 2588.940 3547.280 989.302 1013.2300
## X0dpi.P0.F3 X0dpi.P0.F4 X0dpi.P0.F5 X0dpi.P0.F7 X0dpi.P0.F9
## mt-Nd6 563.855 730.344 290.851 443.5580 246.1650
## mt-Nd4 8288.080 7916.270 5947.040 5544.0200 5474.8800
## mt-Nd3 261.026 211.954 0.000 37.9959 29.7024
## mt-Co3 9319.630 10388.000 3534.500 3513.0000 4685.5100
## mt-Atp8 2406.880 1537.750 903.896 1184.9100 1536.0000
## mt-Co2 2632.830 2831.790 740.061 1993.9100 754.4320
sample_name = colnames(data)
colnames(data) = 1:length(colnames(data))
name = colnames(data)
mapping =cbind(sample_name, name)
mapping
## sample_name name
## [1,] "X0dpi.P0.A2" "1"
## [2,] "X0dpi.P0.B8" "2"
## [3,] "X0dpi.P0.B9" "3"
## [4,] "X0dpi.P0.C11" "4"
## [5,] "X0dpi.P0.C4" "5"
## [6,] "X0dpi.P0.F3" "6"
## [7,] "X0dpi.P0.F4" "7"
## [8,] "X0dpi.P0.F5" "8"
## [9,] "X0dpi.P0.F7" "9"
## [10,] "X0dpi.P0.F9" "10"
## [11,] "X0dpi.P0.G5" "11"
## [12,] "X0dpi.P0.G9" "12"
## [13,] "X0dpi.P1.A1" "13"
## [14,] "X0dpi.P1.A3" "14"
## [15,] "X0dpi.P1.A6" "15"
## [16,] "X0dpi.P1.A8" "16"
## [17,] "X0dpi.P1.B10" "17"
## [18,] "X0dpi.P1.B11" "18"
## [19,] "X0dpi.P1.B4" "19"
## [20,] "X0dpi.P1.B5" "20"
## [21,] "X0dpi.P1.C11" "21"
## [22,] "X0dpi.P1.C1" "22"
## [23,] "X0dpi.P1.C3" "23"
## [24,] "X0dpi.P1.C4" "24"
## [25,] "X0dpi.P1.C6" "25"
## [26,] "X0dpi.P1.D12" "26"
## [27,] "X0dpi.P1.D6" "27"
## [28,] "X0dpi.P1.D8" "28"
## [29,] "X0dpi.P1.E12" "29"
## [30,] "X0dpi.P1.E2" "30"
## [31,] "X0dpi.P1.E5" "31"
## [32,] "X0dpi.P1.E8" "32"
## [33,] "X0dpi.P1.E9" "33"
## [34,] "X0dpi.P1.F11" "34"
## [35,] "X0dpi.P1.F12" "35"
## [36,] "X0dpi.P1.F1" "36"
## [37,] "X0dpi.P1.F3" "37"
## [38,] "X0dpi.P1.F6" "38"
## [39,] "X0dpi.P1.G10" "39"
## [40,] "X0dpi.P1.G3" "40"
## [41,] "X0dpi.P1.G9" "41"
## [42,] "X0dpi.P1.H11" "42"
## [43,] "X0dpi.P1.H12" "43"
## [44,] "X0dpi.P1.H2" "44"
## [45,] "X0dpi.P1.H3" "45"
## [46,] "X0dpi.P1.H6" "46"
Plot expression
sub = data[,c(1:6)]
summary(sub)
## 1 2 3
## Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.00
## Median : 0.00 Median : 0.00 Median : 0.00
## Mean : 54.38 Mean : 32.33 Mean : 62.16
## 3rd Qu.: 0.00 3rd Qu.: 0.00 3rd Qu.: 0.00
## Max. :22625.60 Max. :15010.10 Max. :36232.40
## 4 5 6
## Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 0.00 1st Qu.: 0.00 1st Qu.: 0.00
## Median : 0.00 Median : 0.00 Median : 0.00
## Mean : 29.67 Mean : 36.72 Mean : 35.09
## 3rd Qu.: 0.00 3rd Qu.: 0.00 3rd Qu.: 0.00
## Max. :14183.50 Max. :21192.10 Max. :14648.00
par(mfrow=c(1,2))
hist(log2(data +1))
par(mfrow=c(1,2))
num_cell = length(colnames(data))
num_cell
## [1] 46
for(index in 1:num_cell){
sub_index = data[,index]
sub_index_filter = sub_index[which(sub_index >0)]
hist(log2(sub_index_filter + 1), main=colnames(data)[index])
}
par(mfrow=c(1,1), cex.axis=0.7,cex.lab=0.7)
boxplot(log2(data + 1), names=c(1:num_cell), las=2)
filter_data = data[rowSums(data) > num_cell, ]
par(mfrow=c(1,1), cex.axis=0.7,cex.lab=0.7)
boxplot(log2(filter_data + 1), names=c(1:num_cell), las=2)
Create Seurat object
N0dpi<-CreateSeuratObject(raw.data = data, min.cells=2, min.genes=1000,project="SCI-0dpi")
QC and selecting cells for further analysis
mito.genes <- grep(pattern = "^mt-", x = rownames(x = N0dpi@data), value = TRUE)
percent.mito <- Matrix::colSums(N0dpi@raw.data[mito.genes, ])/Matrix::colSums(N0dpi@raw.data)
#add metadatal
N0dpi <- AddMetaData(object = N0dpi, metadata = percent.mito, col.name = "percent.mito")
VlnPlot(object = N0dpi,features.plot = c("nGene","percent.mito"),nCol=2)
##create gene plot
GenePlot(object = N0dpi, gene1 = "nGene", gene2 = "percent.mito")
#filter cells by gene numbers and mito percentage
N0dpi=FilterCells(object =N0dpi,subset.names = c("nGene","percent.mito"), low.thresholds = c(2500,0.02),high.thresholds = c(6000,Inf))
Normalizing the data
N0dpi=NormalizeData(object = N0dpi,normalization.method = "LogNormalize",scale.factor = 10000)
Detection of variable genes across the single cells
N0dpi=FindVariableGenes(object = N0dpi,mean.function = ExpMean,dispersion.function = LogVMR, x.low.cutoff = 0.25, y.cutoff = 0.5)
abline(h=0.5,col="red")
abline(v=0.25,col="red")
abline(v=10,col="red")
length(x = N0dpi@var.genes)
## [1] 1978
Scaling the data and removing unwanted sources of variation
N0dpi=ScaleData(object = N0dpi, vars.to.regress = c("percent.mito"),display.progress=FALSE)
par(mfrow=c(1,1), cex.axis=0.7,cex.lab=0.7)
boxplot(N0dpi@scale.data, las=2)
Determine statistically significant principal components
N0dpi <- JackStraw(object = N0dpi, num.replicate = 100, do.print = FALSE)
JackStrawPlot(object = N0dpi, PCs = 1:15)
## Warning: Removed 25429 rows containing missing values (geom_point).
PCElbowPlot(object = N0dpi)
Cluster the cells
N0dpi <- FindClusters(object = N0dpi, reduction.type = "pca", dims.use = 1:4, resolution = 0.6, print.output = 0, save.SNN = TRUE)
PrintFindClustersParams(object = N0dpi)
## Parameters used in latest FindClusters calculation run on: 2018-02-04 14:16:18
## =============================================================================
## Resolution: 0.6
## -----------------------------------------------------------------------------
## Modularity Function Algorithm n.start n.iter
## 1 1 100 10
## -----------------------------------------------------------------------------
## Reduction used k.param k.scale prune.SNN
## pca 30 25 0.0667
## -----------------------------------------------------------------------------
## Dims used in calculation
## =============================================================================
## 1 2 3 4
#Run Non-linear dimensional reduction (tSNE)
N0dpi <- RunTSNE(object = N0dpi, dims.use = 1:4, do.fast = TRUE, perplexity=2)
TSNEPlot(object = N0dpi)
save(N0dpi, file = "N0dpi.Robj")
plot(N0dpi@dr$tsne@cell.embeddings) +text(N0dpi@dr$tsne@cell.embeddings, labels=rownames(N0dpi@dr$tsne@cell.embeddings), cex=0.7)
## integer(0)
Plot features
#ependymal
VlnPlot(object = N0dpi, nCol=3, size.x.use = 8, size.y.use = 8, size.title.use = 10, features.plot = c("Vim","Gfap","Nes","S100b","Ncam1","Ncam2","Sox2","Sox9"))
VlnPlot(object = N0dpi, nCol=3, size.x.use = 8, size.y.use = 8, size.title.use = 10, features.plot = c("Fabp7","Pax6","Foxj1","Slc1a3","Hes1","Notch1","Egfr","Aldh1l1"))
#neronal
VlnPlot(object = N0dpi, nCol=3, size.x.use = 8, size.y.use = 8, size.title.use = 10, features.plot = c("Dcx","Tubb3","Prom1","Prom2"))
#oligodendrocytes
VlnPlot(object = N0dpi, nCol=3, size.x.use = 8, size.y.use = 8, size.title.use = 10, features.plot = c("Olig1","Olig2","Mbp","Sox10","Cspg4"))
