scuttle包对单细胞数据质控
2024-07-02 05:31:38
scuttle软件包为单细胞组学数据分析提供了一些实用功能。这包括一些计算和过滤质量控制的简单方法；基础数据转换涉及各种类型的缩放规范化；以及跨单元或功能的灵活聚合。
1. 载入SingleCellExperiment数据
# if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#
# BiocManager::install("scuttle")library(scRNAseq)
sce <- ZeiselBrainData()
sce
2. 增加PerCellQCMetrics数据到colData
library(scuttle)
is.mito <- grep("mt-", rownames(sce))##方法1# subsets:A named list containing one or more vectors (a character vector
# of feature names, a logical vector,
# or a numeric vector of indices), used to identify interesting feature subsets
# such as ERCC spike-in transcripts or mitochondrial genes.
per.cell <- perCellQCMetrics(sce, subsets=list(Mito=is.mito))
colnames(per.cell)
per.cell[1:3,]# 过滤
#discard <- perCellQCFilters(per.cell)
#filtered_sce <- sce[,!discard$discard]# 总counts
summary(per.cell$sum)
# 检测到的基因数量
summary(per.cell$detected)
# sum(assay(sce)[,"1772071015_C02"]>0)summary(per.cell$subsets_Mito_percent)
summary(per.cell$altexps_ERCC_percent)colData(sce) <- cbind(colData(sce), per.cell)
colnames(colData(sce))## 方法2
sce <- addPerCellQCMetrics(sce, subsets=list(Mito=is.mito))
colnames(colData(sce))# sce <- addPerCellQCMetrics(sce)
# colnames(colData(sce))3. 根据CellQC过滤
## 方法1
# An example with just mitochondrial filters.
qc.stats <- perCellQCFilters(per.cell, sub.fields="subsets_Mito_percent")
colSums(as.matrix(qc.stats))# Another example with mitochondrial + spike-in filters.
qc.stats2 <- perCellQCFilters(per.cell, sub.fields=c("subsets_Mito_percent", "altexps_ERCC_percent"))
colSums(as.matrix(qc.stats2))# 过滤
filtered_sce2 <- sce[,!qc.stats2$discard]# 方法2
filtered_sce3 <- quickPerCellQC(sce, subsets=list(Mito=is.mito), sub.fields=c("subsets_Mito_percent", "altexps_ERCC_percent"))4 . 过滤outliers
# Convenience function to determine which values in a numeric vector are
# outliers based on the median absolute deviation (MAD).
# 是否有低counts数（outlier）的单细胞
low.total <- isOutlier(per.cell$sum, type="lower", log=TRUE)
summary(low.total)attr(low.total, "threshold")### 去除总counts数outlier的细胞
sce[,!low.total]# 按tissue分别计算，log2 转化后有离群值。
# 按tissue分别计算，不进行log2 转化，没有离群值。
# 不按tissue分别计算，总体计算，log2 转化后没有离群值
low.total.batched <- isOutlier(per.cell$sum, type="lower", log=TRUE, batch=sce$tissue)
low.total.batched <- isOutlier(per.cell$sum, type="lower", batch=sce$tissue)
low.total.batched <- isOutlier(per.cell$sum, type="lower", log=TRUE)summary(low.total.batched)### 去除总counts数batch outlier的细胞
sce[,!low.total.batched]
5. 查看特征质量矩阵perFeatureQCMetrics
# Pretending that the first 10 cells are empty wells, for demonstration.
# 通过subset设定空wells，结果分别统计空/非空wells中的基因平均值等
per.feat <- perFeatureQCMetrics(sce, subsets=list(Empty=1:10))
summary(per.feat$mean)
summary(per.feat$detected)
summary(per.feat$subsets_Empty_ratio)assay(sce)[,1:10]## 每一个基因的在细胞中的平均count数
#Calculate the average count for each feature after normalizing observations
#using per-cell size factors.
ave <- calculateAverage(sce)
summary(ave)6. 归一化normalization
#####normalization# example_sce <- mockSCE()
# sizeFactors(example_sce)
# example_sce <- computeLibraryFactors(example_sce)
# sizeFactors(example_sce) <- librarySizeFactors(example_sce)
assay(example_sce)
assay(altExp(example_sce))library(scRNAseq)
sce2 <- ZeiselBrainData()
colnames(colData(sce2))library(scuttle)
sce2 <- quickPerCellQC(sce2, subsets=list(Mito=grep("mt-", rownames(sce))),sub.fields=c("subsets_Mito_percent", "altexps_ERCC_percent"))
# 转化为dgCMatrix类型，快速运算
counts(sce) <- as(counts(sce), "dgCMatrix")#Define per-cell size factors from the library sizes (i.e., total sum of counts per cell).
summary(librarySizeFactors(sce))
# Define per-cell size factors from the geometric mean of counts per cell.
#summary(geometricSizeFactors(sce))## 对所有细胞归一化。添加sizeFactor到colData(sce)
sizeFactors(sce) <- librarySizeFactors(sce)
sce <- computeLibraryFactors(sce)
colnames(colData(sce))## 跟据clusters分别计算sizeFactors
library(scran)
clusters <- quickCluster(sce)# Scaling normalization of single-cell RNA-seq data by
# deconvolving size factors from cell pools.
sce <- computePooledFactors(sce, clusters=clusters)
summary(sizeFactors(sce))sce2 <- computeSpikeFactors(sce, "ERCC")
summary(sizeFactors(sce2))# assay(altExp(sce2,"ERCC"))[1:3,1:5]7. 表达谱数据转化
sce <- logNormCounts(sce)
assayNames(sce)
#counts(sce)[1:3,1:2]
#logcounts(sce)[1:3,1:2]assay(sce, "cpm") <- calculateCPM(sce)# calculateFPKM 需要features的lengths8. 其他函数
## 按细胞特征分组统计
agg.sce <- aggregateAcrossCells(sce, ids=sce$level1class)
#applySCE(sce, aggregateAcrossCells,ids=sce$level1class)
head(assay(agg.sce))
colData(agg.sce)[,c("ids", "ncells")]agg.sce <- aggregateAcrossCells(sce, ids=colData(sce)[,c("level1class", "tissue")])
head(assay(agg.sce))
colData(agg.sce)[,c("level1class", "tissue", "ncells")]## 按基因特征分组统计
agg.feat <- sumCountsAcrossFeatures(sce,ids=list(GeneSet1=1:10, GeneSet2=11:50, GeneSet3=1:100),average=TRUE, exprs_values="logcounts")
agg.feat[,1:10]# summarizeAssayByGroup
# aggregateAcrossFeatures参考
http://www.bioconductor.org/packages/release/bioc/html/scuttle.html
scuttle包对单细胞数据质控相关推荐
- SingleR包注释单细胞数据
SingleR包通过利用纯细胞类型的参考转录组数据集独立推断每个单细胞的起源细胞,从单细胞RNA测序数据执行无偏细胞类型识别. 1. 载入单细胞数据 # if (!requireNamespace(& ...
- splatter包生成单细胞RNA测序数据
Splatter是一个模拟单细胞RNA测序计数数据的软件包.它提供了一个简单的界面,用于创建可复制且文档充分的复杂模拟.可以从真实数据估计参数,并提供用于比较真实数据集和模拟数据集的函数. # if ...
- 这个R包自动注释单细胞数据的平均准确率为83％，使用后我的结果出现了点问题|附全代码...
估计大家都有一个这样的感觉就是单细胞数据具有一定的数据依赖性,好多的marker在相同的组织中,别人的数据就表达的十分明显,在你的数据中就是不太显著,比如NK细胞的KLRF1.于是,细胞自动注释也就应 ...
- R语言分析单细胞数据Day1——下载Seurat包并进行预处理（一）
Task.1 安装Seurat,准备处理single cell data 安装Seurat时,只能安装3.2.3以下的版本,太高就不兼容! install.packages('remotes') %安 ...
- 代码分析 | 单细胞转录组质控详解
前言 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 (重磅综述:三万字 ...
- 推荐几个单细胞数据分享和展示平台 | 短视频演示
Broad的单细胞数据分享和展示平台 可选择子类展示 映射单个基因的颜色到t-SNE/UMAP图 分屏展示Cluster着色图和单基因着色图 多基因热图.Dotplot.Boxplot.Violinp ...
- 一个R包完成单细胞基因集富集分析 （全代码）
singleseqgset | 单细胞RNA-Seq基因集富集分析 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (Ch ...
- 让你的单细胞数据动起来！|iCellR(二)
前情回顾:让你的单细胞数据动起来!|iCellR(一) 在上文我们已经基本完成了聚类分析和细胞在不同条件下的比例分析,下面继续开始啦,让你看看iCellR的强大! 1. 每个cluster的平均表达量 ...
- 【Bioinfo Blog 006】【R Code 005】——GEO表达谱数据质控
目录 一.数据下载及质控 1.1 GEO数据下载 1.1.1 GEOquery包安装 1.1.2 .cel数据下载 1.2 读取.cel文件 1.2.1 Affy包安装 1.2.2 利用ReadAff ...
最新文章
- java datahandler_Java Web Services：使用DataHandler类发送文件
- TortoiseGit入门（图文教程） Git，Github，puttygen，SSH
- Android图片缓存框架Glide
- 觉得UtraWebGrid老不稳定
- aix用户登录次数受限问题（3004-300 输入了无效的登录名或password）
- Lucene实战(第2版)》
- JAVA开发工具整理
- BZOJ4912 SDOI2017天才黑客（最短路+虚树）
- 商用台式电脑配置_装机不求人，10分钟电脑配置挑选速成攻略
- python异常处理与导入模块与导入包
- 在WCF中使用消息队列MSMQ
- Linux Vue环境搭建
- 在macOS下制作黑苹果镜像
- 2020年中级数据库系统工程师考试时间表与考试大纲
- csrss32.exe
- 短视频app开发，随机生成中文名字
- 微信公共号开发简单入门
- 计算机学渣和你说说从毕业到工作
- linux mate主题目录,七大顶级Linux桌面：Cinnamon和MATE_服务器_服务器产业-中关村在线...
- 国内外几个主流的CMS系统推荐