4.0 安裝包
install.packages('Seurat')
加載所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
清除環(huán)境静汤,
rm(list=ls())
讀取10x的數(shù)據(jù)
BC10 文件夾內(nèi)容
4640_0.png
scRNA.counts=Read10X("F:/006_data/BC10")
class(scRNA.counts)
scRNA.counts是一個S4對象逆甜。S4對象就像一顆大樹一樣教藻,有主枝干篡悟,有分枝谜叹。
4642_0.png
創(chuàng)建Seurat對象
?CreateSeuratObject
scRNA = CreateSeuratObject(scRNA.counts ,min.cells = 3,project="os", min.features = 300)
view(scRNA)
min.cells 每個基因至少要在3個細胞中表達
min.features 每個細胞至少有多少基因表達/被檢測到
查看樣本的細胞數(shù)量
第一種s4對象的提取方法 點擊白框
第二提取s4對象的方法 @ $交替使用
一般來說S4第一個用@
先輸入枝干scRNA@assays$RNA@counts@Dim
table(scRNA@meta.data$orig.ident)
計算質(zhì)控指標匾寝,去除低質(zhì)量的細胞
計算細胞中線粒體基因比例,線粒體基因公認是一定要去除的荷腊。
1 線粒體是獨立遺傳的艳悔,不是染色體上基因控制的。
2 單細胞測序一般都是新鮮的組織女仰,防止mRNA降解猜年,臨床新鮮的。防止線粒體快速擴增疾忍。 不新鮮乔外,線粒體不會占比太高,線粒體占比太高一罩,細胞可能凋亡或壞死杨幼。
scRNA[[]]
scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
x <- PercentageFeatureSet(scRNA, pattern = "^MT-")
x列名是BULK細胞標簽,行名是線粒體在每個基因的占比
4646_1.png
線粒體基因都是用MT-開頭命名的聂渊。
用正則表達式^符號表示從首字母開始MT進行匹配提取差购。
計算紅細胞比例,紅細胞沒有細胞核歧沪,沒有轉(zhuǎn)錄組
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)
col.num <- length(levels(scRNA@active.ident))
Feature歹撒、count莲组、線粒體基因诊胞、紅細胞基因占比可視化。
violin <- VlnPlot(scRNA,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.HB"),
cols =rainbow(col.num),
pt.size = 0.01, #不需要顯示點锹杈,可以設置pt.size = 0
ncol = 4) +
theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())
把圖片畫到畫板上面
violin
4644_0.png
保存
ggsave("vlnplot_before_qc.pdf", plot = violin, width = 12, height = 6)
ggsave("vlnplot_before_qc.png", plot = violin, width = 12, height = 6)
這幾個指標之間的相關性撵孤。 把圖畫到畫板上,然后手動保存
plot1=FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2=FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot3=FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.HB")
pearplot <- CombinePlots(plots = list(plot1, plot2, plot3), nrow=1, legend="none")
看畫板
plot1
plot2
plot3
pearplot
我們可以看到竭望,nFeature_RNA的范圍在0到8000之內(nèi)邪码,每個細胞基因表達量
percent.mt代表線粒體含量
percent.HB代表紅細胞含量
我們默認線粒體含量至少要小于10%,這是根據(jù)生物學知識得出的默認閾值咬清。紅細胞的數(shù)目要至少小于3%
至于nFeature_RNA和nCount_RNA的閾值怎么確定闭专,這個要結合 pearplot的圖來判斷。我們質(zhì)控的目標就是刪除離異值旧烧。而且注意閾值盡可能取的寬松一下影钉,防止后面分析想要的細胞得不到。
接下來從pearplot的圖片來做質(zhì)控---剔除離異值
nFeature_RNA選擇大于300 小于7000的 nFeature_RNA選擇小于100000掘剪,percent.mt小于10平委,percent.HB小于3
scRNA1 <- subset(scRNA, subset = nFeature_RNA > 300& nFeature_RNA < 7000 & percent.mt < 10 & percent.HB < 3 & nCount_RNA < 100000)
scRNA
scRNA1
在控制臺中我們可以看到有500多細胞過濾了
過濾完之后 我們就要對數(shù)據(jù)進行均一化,使用NormalizeData這個函數(shù)夺谁。
注意均一化是用NormalizeData廉赔,標準化是用ScaleData
scRNA1 <- NormalizeData(scRNA1, normalization.method = "LogNormalize", scale.factor = 10000)
好了肉微,這一節(jié)數(shù)據(jù)加載、質(zhì)控的內(nèi)容就算是做完了蜡塌。
在我們關閉rstudio之前 先把環(huán)境中運行好的數(shù)據(jù)保存一下
數(shù)據(jù)將保存在之前設定好的路徑中碉纳。還有保存的scRNA1,不是scRNA馏艾,因為scRNA1才是過濾好的數(shù)據(jù)村象。
save(scRNA1,file='scRNA1.Rdata')