Seurat應(yīng)用筆記 數(shù)據(jù)為Seurat官網(wǎng)示范數(shù)據(jù)
library(Seurat)
library(dplyr)
library(magrittr)
每次重新進(jìn)入R都要重新library包
每次得到什么中間結(jié)果可以
save.image("保存路徑")
或者直接點(diǎn)environment左上方的保存按鈕
（當(dāng)然可以一開始用setwd設(shè)定好工作位置足淆，然后調(diào)用save.image的時(shí)候就可以直接用“文件名+.RData”的形式），就可以對每次的數(shù)據(jù)分開保存。
load的作用就是把一個(gè)RData類型文件載入到environment中埃篓，這樣變量就全部在當(dāng)前環(huán)境下雕崩，可以直接使用
setwd("運(yùn)行環(huán)境位置")
單細(xì)胞數(shù)據(jù)文件分類：
由單細(xì)胞得的10x數(shù)據(jù)，怎么處理？
code1：
matrix <- Read10X_h5("filtered_gene_bc_matrices_h5.h5")
該函數(shù)直接讀取10x的h5文件枣抱，并且得到一個(gè)matrix矩陣
code2：
matrix <- Read10X("data/matrix_dir")
默認(rèn)讀取運(yùn)行環(huán)境中的 data/matrix_dir目錄下包含文件布疙，如果要讀取運(yùn)行環(huán)境之外的文件蚊惯，就要從頭什么盤開始輸入位置。
包括matrix.mtx, genes.tsv, and barcodes.tsv 把他們都放到matrix中
導(dǎo)入文件后用Seurat包處理
pbmc <-CreateSeuratObject(raw.data = matrix,min.cells = 3,
min.genes = 200, project = "NAME")
創(chuàng)建一個(gè)叫做"pbmc"的東西拐辽，是SeuratObject拣挪。
其中命名raw.data為以上matrix（有的代碼中命名不同，或者一開始導(dǎo)入的時(shí)候不以matrix形式導(dǎo)入俱诸，而以counts形式菠劝，總之這里raw.data = matrix可寫“你要?jiǎng)?chuàng)建的文件的名字”=“你前面導(dǎo)入東西的名字”）；基因至少在min.cells個(gè)細(xì)胞中表達(dá)睁搭；每個(gè)細(xì)胞中至少表達(dá)min.genes個(gè)基因赶诊，可以自己設(shè)定；project的命名自己定義就好园骆。
之后Seurat將數(shù)據(jù)保存在不同的slot中舔痪，如pbmc@raw.data, pbmc@data, pbmc@meta.data, pbmc@ident，其中raw.data存放的是每個(gè)細(xì)胞中每個(gè)gene的原始UMI數(shù)據(jù)锌唾，data存放的是gene的表達(dá)量锄码，meta.data存放的是每個(gè)細(xì)胞的統(tǒng)計(jì)數(shù)據(jù)如UMI數(shù)目，gene數(shù)目等晌涕，ident此時(shí)存放的是project信息滋捶。
由于技術(shù)原因，一個(gè)GEM（單細(xì)胞測試的油包水滴）中可能會(huì)包含2個(gè)或多個(gè)細(xì)胞余黎，也可能不包含細(xì)胞重窟，這時(shí)候可以通過觀察每個(gè)barcode中的基因數(shù)目或UMI數(shù)目來判斷。
處理篩掉異常GEM的方法：
- 可以計(jì)算每個(gè)barcode中的線粒體基因含量等惧财，從而更加仔細(xì)的觀察數(shù)據(jù)的質(zhì)量巡扇。
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^mt-")
mt-開頭的是線粒體基因扭仁，這里將其進(jìn)行標(biāo)記并統(tǒng)計(jì)其分布頻率
VlnPlot(object = pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")+ NoLegend()
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")+ NoLegend()
CombinePlots(plots = list(plot1, plot2))
以上對 pbmc 對象做小提琴圖，分別為基因數(shù)厅翔，細(xì)胞數(shù)和線粒體占比
以上根據(jù)圖片中基因數(shù)和線粒體數(shù)乖坠，分別設(shè)置過濾參數(shù)，這里基因數(shù)200-2500知给，線粒體百分比為小于 5%（常用的過濾參數(shù)）
作圖瓤帚，直觀地看一下數(shù)據(jù)的質(zhì)量。要輸入代碼plot或者Combineplot涩赢，才可以在右邊出現(xiàn)圖
標(biāo)準(zhǔn)化
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
鑒定高變基因戈次，之后就關(guān)注這些基因進(jìn)行研究
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
top10 <- head(VariableFeatures(pbmc), 10)
使用默認(rèn)參數(shù)，用vst方法選取2000個(gè)高變基因筒扒。
然后選取top10變化的
plot3 <- VariableFeaturePlot(pbmc)
plot4 <- LabelPoints(plot = plot3, points = top10, repel = TRUE,xnudge=0,ynudge=0)
plot3+plot4
得到圖怯邪，看到變化最大的前2000個(gè)和10個(gè)的圖
接下來進(jìn)行數(shù)據(jù)歸一化
1.線性姜維
有一堆細(xì)胞，但是你需要將這些細(xì)胞進(jìn)行分群花墩，但是你又不知道這些細(xì)胞如果進(jìn)行分群的話是按照什么依據(jù)進(jìn)行分群悬秉，就要先使用PCA 降維找到細(xì)胞分群的主要特點(diǎn)，這也就是我們常說的主成分分析
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc,features = all.genes,vars.to.regress = "percent.mt")
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
數(shù)據(jù)歸一化冰蘑，對所有基因進(jìn)行標(biāo)準(zhǔn)化和泌，默認(rèn)只是標(biāo)準(zhǔn)化高變基因（ 2000 ），速度更快祠肥，不影響 PCA 和分群武氓，但影響熱圖的繪制。
線性降維（PCA）仇箱，默認(rèn)用前面得到的2000個(gè)高變基因集县恕，但也可通過 features 參數(shù)自己指定 需要挺長時(shí)間的，從1%開始加載
定義可視化細(xì)胞和功能的幾種有用的方式PCA剂桥，包括VizDimReduction忠烛，DimPlot，和DimHeatmap 得到一堆positive negative基因
VizDimLoadings(object = pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")+ NoLegend()
DimHeatmap(pbmc, dims = 1:2, cells = 500, balanced = TRUE)
以上三行得到三個(gè)圖权逗，定義可視化細(xì)胞和功能美尸，還不知道有啥用？
鑒定數(shù)據(jù)集的可用維度斟薇，確定哪些主成分所代表的基因可以進(jìn)入下游分析用于后續(xù)細(xì)胞分類火惊。
這里可以使用JackStraw做重抽樣分析”伎眩可以用JackStrawPlot可視化看看哪些主成分可以進(jìn)行下游分析。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)
ElbowPlot(object = pbmc)
第一步需要一定計(jì)算時(shí)間最后得到圖JackStrawPlot和Elbowplot尸疆，后者比較直觀椿猎。
第二行惶岭，得到JackStrawPlot。是彩色的犯眠，虛線以上的為可用維度按灶，dim是展示多少個(gè)】疬郑可以調(diào)節(jié)dim鸯旁，畫出不同數(shù)量的 pca 查看。
也就是：基于每個(gè)主成分對方差解釋率的排名量蕊。建議嘗試選擇多個(gè)主成分個(gè)數(shù)做下游分析铺罢，對整體影響不大；在選擇此參數(shù)時(shí)残炮，建議選擇偏高的數(shù)字（ “寧濫勿缺”韭赘，為了獲取更多的稀有分群）；有些亞群很罕見势就，如果沒有先驗(yàn)知識(shí)泉瞻，很難將這種大小的數(shù)據(jù)集與背景噪聲區(qū)分開來。
第三行苞冯，得到Elbowplot袖牙，點(diǎn)圖。比較直觀地看到舅锄，前面的幾個(gè)點(diǎn)分得比較開鞭达。
之后就可以決定選取多少主成分用于細(xì)胞分類。從上到下巧娱，從p-value最小的開始選擇top few碉怔。
尊重官網(wǎng)的建議，這里選取top 10
- 非線性降維（ UMAP/tSNE)
TSNE降維與PCA降維并不相同禁添，TSNE降維主要是非線性降維撮胧。
TSNE優(yōu)勢：使得相似的對象有更高的概率被選擇，而不相似的對象有較低的概率被選擇老翘。它將多維數(shù)據(jù)映射到適合于人類觀察的兩個(gè)或多個(gè)維度芹啥，廣泛應(yīng)用于圖像處理。
一般情況下铺峭，在進(jìn)行細(xì)胞分群時(shí)：先進(jìn)行PCA主成分分析墓怀，再進(jìn)行TSNE降維分析進(jìn)行細(xì)胞分群，這樣的分群結(jié)果可靠度更高卫键。
Seurat 提供了幾種非線性降維的方法進(jìn)行數(shù)據(jù)可視化（在低維空間把相似的細(xì)胞聚在一起）傀履，比如 UMAP 和 t-SNE，運(yùn)行 UMAP 需要先安裝'umap-learn'包莉炉，兩種方法都可以使用钓账，但不要混用碴犬，這樣，后面的結(jié)算結(jié)果會(huì)將先前的聚類覆蓋掉梆暮，只能保留一個(gè)服协。
2.1 umap
umap優(yōu)于tsne：
UMAP 與 t-SNE 均是非線性降維，多用于數(shù)據(jù)可視化
UMAP 結(jié)構(gòu)與t-SNE一致
UMAP 計(jì)算更快
UMAP 能更好地反映高緯結(jié)構(gòu)啦粹，比t-SNE有著更好的連續(xù)性
這種連續(xù)性反映到單細(xì)胞分析中就是能更好滴可視化細(xì)胞的分化狀態(tài)（UMAP better represents the multi-branched continuous trajectory of hematopoietic development）
umap安裝包如下：
install(packages = 'umap-learn')
基于PCA 空間中的歐氏距離計(jì)算 nearest neighbor graph偿荷，優(yōu)化任意兩個(gè)細(xì)胞間的距離權(quán)重（輸入上一步得到的 PC 維數(shù)）
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
resolution 參數(shù)決定下游聚類分析得到的分群數(shù)，對于 3K 左右的細(xì)胞唠椭，設(shè)為 0.4-1.2 能得到較好的結(jié)果(官方說明)跳纳；如果數(shù)據(jù)量增大，該參數(shù)也應(yīng)該適當(dāng)增大泪蔫；增加的值會(huì)導(dǎo)致更多的群集棒旗。
使用 Idents（）函數(shù)可查看不同細(xì)胞的分群；查看每一類有多少個(gè)細(xì)胞
head(Idents(pbmc), 8)
table(pbmc@active.ident)
得到前8個(gè)細(xì)胞分群的table
pbmc <- RunUMAP(object = pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = "umap")
umapplot<-UMAPPlot(pbmc,label = TRUE, pt.size = 1.5)+ NoLegend()
這里用 DimPlot()函數(shù)繪制散點(diǎn)圖撩荣。如果只做tsne铣揉，參數(shù)reduction = "tsne"；如果不設(shè)定reduction餐曹，默認(rèn)先從搜索 umap逛拱， 然后 tsne, 再然后 pca；
也可以直接使用這 3 個(gè)函數(shù) PCAPlot()台猴、 TSNEPlot()朽合、UMAPPlot()； cols饱狂， pt.size 分別調(diào)整分組顏色和點(diǎn)的大胁懿健；
2.2 TSNE
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(object = pbmc, reduction = "tsne")
tsneplot<-TSNEPlot(pbmc,label = TRUE, pt.size = 1.5)+ NoLegend()
做完以上降維處理保存數(shù)據(jù)
saveRDS(pbmc, file = "pbmc_tutorial.rds")
save(pbmc,file="pbmc.RData")
load(file = "pbmc.RData")
保存在了工作路徑下文件夾中休讳，并且下次用的時(shí)候直接load就好
尋找差異表達(dá)的特征（聚類生物標(biāo)志物）
尋找某個(gè)聚類和其他所有聚類顯著表達(dá)的基因
cluster1.markers <- FindMarkers(object = pbmc, ident.1 = 1, min.pct = 0.25)
head(x = cluster1.markers, n = 5)
find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(object = pbmc, ident.1 = 5, ident.2 = c(0,3), min.pct = 0.25)
write.table(as.data.frame(cluster5.markers), file="cluster1vs2.xls", sep="\t", quote = F)
find markers for every cluster compared to all remaining cells, report only the positive ones
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)
pbmc.markers <- FindAllMarkers(object = pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25, return.thresh = 0.01)
library(dplyr)
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(object = pbmc, features = top10$gene) + NoLegend()
畫圖
VlnPlot(object = pbmc, features = c("MS4A1", "CD79A"))
FeaturePlot(object = pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
you can plot raw UMI counts as well
VlnPlot(object = pbmc, features.plot = c("NKG7", "PF4"), use.raw = TRUE, y.log = TRUE)
FeaturePlot(object = pbmc, features = c("hbaa1","hbaa2","ighv1-4","cd79a","cd79b","cd4-1","cd8a","cd8b","itga2b"), cols = c("grey", "blue"), reduction = "tsne")
將單元類型標(biāo)識(shí)分配給集群
new.cluster.ids <- c("1", "3", "4", "2", "5", "6", "7", "8")
names(x = new.cluster.ids) <- levels(x = pbmc)
pbmc <- RenameIdents(object = pbmc, new.cluster.ids)
DimPlot(object = pbmc, label = TRUE, pt.size = 1.5) + NoLegend()
氣泡圖
markers.to.plot <- c("cd74a", "cd74b")
DotPlot(pbmc, features = markers.to.plot) + RotatedAxis()
