第一步 數(shù)據(jù)獲取與前期準(zhǔn)備
從NCBI上下載GEO數(shù)據(jù)后解壓拉庶，在Rstudio中加載需要的package
##系統(tǒng)報錯改為英文
Sys.setenv(LANGUAGE = "en")
##禁止轉(zhuǎn)化為因子
options(stringsAsFactors = FALSE)
##清空環(huán)境
rm(list=ls())
##安裝Seurat包
install.packages('Seurat')
##改變工作路徑，將相關(guān)結(jié)果儲存至對應(yīng)文件夾
setwd("C:/Users/86269/Desktop/shun.C/single_cell")
##加載需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
第二步 讀取數(shù)據(jù)并進行質(zhì)控
讀取數(shù)據(jù)
##讀取10×數(shù)據(jù)，文件夾中需要包括features.tsv,barcodes.tsv,matrix.mtx,分別代表基因名祷愉，標(biāo)記序列，表達量矩陣
scRNA_counts = Read10X("C:/Users/86269/Desktop/shun.C/single_cell/BC21")
class(scRNA_counts)
##創(chuàng)建Seurat對象
##min.cells表示基因至少要在多少個細(xì)胞中被檢測到表達才算有效
##min.features至少要檢測到多少個基因才算有效
scRNA = CreateSeuratObject(scRNA_counts,project = "sample_21",
min.cells = 3, min.features = 300)
進行質(zhì)控
過濾線粒體DNA含量過高的基因执解，因為線粒體基因含量過高表示細(xì)胞即將凋亡或狀態(tài)不佳帝洪，過濾標(biāo)準(zhǔn)需要根據(jù)實際情況考慮，例如肝臟細(xì)胞與肌肉細(xì)胞本身含線粒體量較多做修，過濾閾值也應(yīng)當(dāng)適當(dāng)增大霍狰。
過濾紅細(xì)胞，將紅細(xì)胞表達的標(biāo)志性基因與樣本基因匹配饰及，如果一個細(xì)胞的這些基因表達量過高蔗坯，就將該細(xì)胞視作紅細(xì)胞。
#計算線粒體基因比例
scRNA[["percent.mt"]]=PercentageFeatureSet(scRNA,pattern = "^MT-")
#計算紅細(xì)胞比例
HB.genes = c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
HB_m <- match(HB.genes,rownames(scRNA@assays$RNA))
HB.genes = rownames(scRNA@assays$RNA)[HB_m]
HB.genes = HB.genes[!is.na(HB.genes)]
scRNA[["percent.HB"]] <- PercentageFeatureSet(scRNA,features = HB.genes)
#過濾低質(zhì)量細(xì)胞
scRNA1 <- subset(scRNA,subset=nFeature_RNA>500 &
nCount_RNA>1000 & percent.mt<20 & percent.HB<1)
第三步 歸一化燎含，降維與聚類
1.歸一化宾濒，排除測序深度的影響
2.降維與聚類
- 尋找高變基因，即在各細(xì)胞中表達量差異最大的基因屏箍，一般挑選3000個绘梦，這一步將維度降到3000
- 中心化處理，PCA降維前的必要步驟
- 對高變基因通過PCA降維赴魁，一般降至50個維度
- 制無向有權(quán)圖并切圖卸奉，對細(xì)胞進行聚類
- 進行降維，降至二維進行展示颖御，有tsne與umap法
#歸一化處理
scRNA1 <- NormalizeData(scRNA1,normalization.method="LogNormalize",scale.factor=10000)
#挑選高變基因
scRNA1 <- FindVariableFeatures(scRNA1,selection.method="vst",nfeatures=3000)
#對基因進行中心化處理榄棵，將每個細(xì)胞中高變基因的表達量改變，使細(xì)胞的平均表達為0潘拱，細(xì)胞間差異為1疹鳄，成為標(biāo)準(zhǔn)數(shù)據(jù)，賦予每個基因相同的權(quán)重泽铛，因此高表達基因不占優(yōu)勢
scale_gene <- VariableFeatures(scRNA1)
scRNA1 <- ScaleData(scRNA1,features=scale_gene)
#將高變基因通過PCA降維,RunPCA默認(rèn)選取50個PC
scRNA1 <- RunPCA(scRNA1,features=VariableFeatures(scRNA1))
##聚類處理
pc.num=1:20
scRNA1 <- FindNeighbors(scRNA1,dims=pc.num)
scRNA1 <- FindClusters(scRNA1,resolution=1.0)
##可視化聚類尚辑，將高維數(shù)據(jù)降低到二維or三維進行展示
scRNA1 <- RunTSNE(scRNA1,dims=pc.num)
embed_tsne <- Embeddings(scRNA1,"tsne")
#tsne法
DimPlot(scRNA1,reduction="tsne",label=TRUE)
#umap法
DimPlot(scRNA1,reduction="umap",label=TRUE)
第四步 細(xì)胞周期評分
scRNA測序中，單個樣品里有許多細(xì)胞盔腔，每個細(xì)胞的表達的基因種類與表達量都不相同杠茬，根據(jù)基因表達的不同將細(xì)胞進行分類月褥，但細(xì)胞基因表達的差距除了是因為細(xì)胞種類的差別，也有可能是因為細(xì)胞所處周期的差別瓢喉，細(xì)胞聚類可能會被細(xì)胞周期影響宁赤，所以需要檢查細(xì)胞周期的影響，如果細(xì)胞周期影響較大則需要排除細(xì)胞周期因素
#細(xì)胞周期評分
#cc.genes是Seurat包自帶的一個list栓票，含不同細(xì)胞周期的marker基因
CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA1))
g2m_genes=cc.genes$g2m.genes
g2m_genes=CaseMatch(search = g2m_genes,match=rownames(scRNA1))
s_genes=cc.genes$s.genes
s_genes=CaseMatch(search=s_genes,match=rownames(scRNA1))
scRNA1 <- CellCycleScoring(object=scRNA1,g2m.features = g2m_genes,s.features = s_genes)
#查看細(xì)胞周期基因?qū)?xì)胞聚類的影響
scRNAa <- RunPCA(scRNA1,features=c(s_genes,g2m_genes))
p <- DimPlot(scRNAa,reduction = "pca",group.by="Phase")
p
#如果細(xì)胞周期有影響則消除細(xì)胞周期的影響
#scRNAb <- ScaleData(scRNA1,vars.to.regress="S.Score","G2M.Score")
細(xì)胞周期影響的判定需要在中心化處理后决左，排除細(xì)胞周期影響后就可以繼續(xù)進行PCA降維等后續(xù)步驟
