单细胞转录组数据校正批次效应实战（中）

刘小泽写于19.6.24、29
和花花毕业流浪，武汉完事转向成都
继续进行剩下两个数据集

第二个数据--CEL-seq2, GSE85241

Muraro et al. (2016) 利用CEL-seq2技术并结合UMI、ERCC得到的

https://www.ncbi.nlm.nih.gov//geo/query/acc.cgi?acc=GSE85241

下面快速使用代码

读数据，看数据

gse85241.df <- read.table("GSE85241_cellsystems_dataset_4donors_updated.csv.gz", sep='\t', header=TRUE, row.names=1)
> dim(gse85241.df)
[1] 19140  3072

提取meta信息

# 还是先看一下
> head(colnames(gse85241.df))
[1] "D28.1_1" "D28.1_2" "D28.1_3" "D28.1_4" "D28.1_5" "D28.1_6"
# 依然是：点号前面的是donor信息
donor.names <- sub("^(D[0-9]+).*", "\\1", colnames(gse85241.df))
> table(donor.names)
donor.names
D28 D29 D30 D31 
768 768 768 768 
# 然后文章使用了8个96孔板，于是可以将点号和下划线之间的数字提取出来
plate.id <- sub("^D[0-9]+\\.([0-9]+)_.*", "\\1", colnames(gse85241.df)) #这句代码中注意使用了一个转义符\\，在R中需要用两个反斜线来表示转义
> table(plate.id)
plate.id
  1   2   3   4   5   6   7   8 
384 384 384 384 384 384 384 384

提取基因、ERCC信息

gene.symb <- gsub("__chr.*$", "", rownames(gse85241.df))
is.spike <- grepl("^ERCC-", gene.symb)
> table(is.spike)
is.spike
FALSE  TRUE 
19059    81

基因转换

library(org.Hs.eg.db)
gene.ids <- mapIds(org.Hs.eg.db, keys=gene.symb, keytype="SYMBOL", column="ENSEMBL")
gene.ids[is.spike] <- gene.symb[is.spike]

keep <- !is.na(gene.ids) & !duplicated(gene.ids)
gse85241.df <- gse85241.df[keep,]
rownames(gse85241.df) <- gene.ids[keep]
> summary(keep)
   Mode   FALSE    TRUE 
logical    1949   17191 
# 去掉了快2000个重复或无表达的基因

创建单细胞对象

# 存储metadata作为colData、基因信息作为rawData、ERCC作为spike-in
sce.gse85241 <- SingleCellExperiment(list(counts=as.matrix(gse85241.df)),
                                     colData=DataFrame(Donor=donor.names, Plate=plate.id),
                                     rowData=DataFrame(Symbol=gene.symb[keep]))
isSpike(sce.gse85241, "ERCC") <- grepl("^ERCC-", rownames(gse85241.df))

质控和标准化

sce.gse85241 <- calculateQCMetrics(sce.gse85241, compact=TRUE)
QC <- sce.gse85241$scater_qc
low.lib <- isOutlier(QC$all$log10_total_counts, type="lower", nmad=3)
low.genes <- isOutlier(QC$all$log10_total_features_by_counts, type="lower", nmad=3)
high.spike <- isOutlier(QC$feature_control_ERCC$pct_counts, type="higher", nmad=3)
data.frame(LowLib=sum(low.lib), LowNgenes=sum(low.genes), 
           HighSpike=sum(high.spike, na.rm=TRUE))

#  LowLib LowNgenes HighSpike
#    577       669       696

# 然后去掉低质量的细胞
discard <- low.lib | low.genes | high.spike
sce.gse85241 <- sce.gse85241[,!discard]
> summary(discard)
   Mode   FALSE    TRUE 
logical    2346     726

可以看到文库小的有577个，基因表达少的有669个，高spike-in的有696个，但是最后只去掉了726个，这是因为，有的细胞同时存在以上两种或三种低质量情况，因此并不能简单认为总共去除细胞数=577+669+696

聚类

clusters <- quickCluster(sce.gse85241, min.mean=0.1, method="igraph")
> table(clusters)
clusters
  1   2   3   4   5   6 
237 248 285 483 613 480

标准化

sce.gse85241 <- computeSumFactors(sce.gse85241, min.mean=0.1, clusters=clusters)
summary(sizeFactors(sce.gse85241))
sce.gse85241 <- computeSpikeFactors(sce.gse85241, general.use=FALSE)
summary(sizeFactors(sce.gse85241, "ERCC"))
sce.gse85241 <- normalize(sce.gse85241)

鉴定HVGs

block <- paste0(sce.gse85241$Plate, "_", sce.gse85241$Donor)
fit <- trendVar(sce.gse85241, block=block, parametric=TRUE) 
dec <- decomposeVar(sce.gse85241, fit)
plot(dec$mean, dec$total, xlab="Mean log-expression", 
     ylab="Variance of log-expression", pch=16)
is.spike <- isSpike(sce.gse85241)
points(dec$mean[is.spike], dec$total[is.spike], col="red", pch=16)
curve(fit$trend(x), col="dodgerblue", add=TRUE)

这张图中的ERCC表达量就有一些比较高的，但是占比不高，另外总体波动不大

# 选出来这些基因
dec.gse85241 <- dec
dec.gse85241$Symbol <- rowData(sce.gse85241)$Symbol
dec.gse85241 <- dec.gse85241[order(dec.gse85241$bio, decreasing=TRUE),]
> head(dec.gse85241,2)
DataFrame with 2 rows and 7 columns
                            mean            total              bio
                       <numeric>        <numeric>        <numeric>
ENSG00000115263 7.66453729345785 6.66863456231166 6.63983282676052
ENSG00000089199 4.62375793902937 6.46558866721711 6.34422879524839
                              tech   p.value       FDR      Symbol
                         <numeric> <numeric> <numeric> <character>
ENSG00000115263 0.0288017355511366         0         0         GCG
ENSG00000089199   0.12135987196872         0         0        CHGB

第三个数据--Smart-seq2, E-MTAB-5061

Segerstolpe et al. (2016)利用Smart-seq2，添加了ERCC，这个数据和上面两个不同，它存放在ArrayExpress数据库，当然也是用链接规律的：https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-5061/ (这个文件比较大，压缩文件151M，解压后700多M)

读入数据

文件较大，先读入样本，也就是第一行(nrow=1)，看下数量

header <- read.table("pancreas_refseq_rpkms_counts_3514sc.txt", 
    nrow=1, sep="\t", comment.char="")
# 先看下header信息
> header[1,1:4]
        V1            V2            V3            V4
 #samples HP1502401_N13 HP1502401_D14 HP1502401_F14
# 然后将第一个（#samples）去掉
ncells <- ncol(header) - 1L #保存为整数

然后只加载基因名称和表达矩阵

# 这段代码需要再好好理解下
col.types <- vector("list", ncells*2 + 2)
col.types[1:2] <- "character"
col.types[2+ncells + seq_len(ncells)] <- "integer"
e5601.df <- read.table("pancreas_refseq_rpkms_counts_3514sc.txt", 
    sep="\t", colClasses=col.types)
# 最后将基因信息和表达矩阵分离
gene.data <- e5601.df[,1:2]
e5601.df <- e5601.df[,-(1:2)]
colnames(e5601.df) <- as.character(header[1,-1])
dim(e5601.df)
## [1] 26271  3514

判断ERCC

# gene.data[,2]对应测序数据中的基因ID，gene.data[,1]是相应的symbol ID
is.spike <- grepl("^ERCC-", gene.data[,2])
> table(is.spike)
is.spike
FALSE  TRUE 
26179    92

基因转换

library(org.Hs.eg.db)
gene.ids <- mapIds(org.Hs.eg.db, keys=gene.data[,1], keytype="SYMBOL", column="ENSEMBL")
gene.ids[is.spike] <- gene.data[is.spike,2]
# 去掉重复和无表达基因
keep <- !is.na(gene.ids) & !duplicated(gene.ids)
e5601.df <- e5601.df[keep,]
rownames(e5601.df) <- gene.ids[keep]
> summary(keep)
   Mode   FALSE    TRUE 
logical    3367   22904

提取metadata信息

metadata <- read.table("E-MTAB-5061.sdrf.txt", header=TRUE, 
                       sep="\t", check.names=FALSE)
m <- match(colnames(e5601.df), metadata$`Assay Name`)
stopifnot(all(!is.na(m)))
metadata <- metadata[m,]
donor.id <- metadata[["Characteristics[individual]"]]
> table(donor.id)
donor.id
          AZ    HP1502401 HP1504101T2D    HP1504901    HP1506401 
          96          352          383          383          383 
   HP1507101 HP1508501T2D    HP1509101 HP1525301T2D HP1526901T2D 
         383          383          383          384          384

创建单细胞对象

sce.e5601 <- SingleCellExperiment(list(counts=as.matrix(e5601.df)),
    colData=DataFrame(Donor=donor.id),
    rowData=DataFrame(Symbol=gene.data[keep,1]))
isSpike(sce.e5601, "ERCC") <- grepl("^ERCC-", rownames(e5601.df))

后面的操作和之前保持一致了

sce.e5601 <- calculateQCMetrics(sce.e5601, compact=TRUE)
QC <- sce.e5601$scater_qc
low.lib <- isOutlier(QC$all$log10_total_counts, type="lower", nmad=3)
low.genes <- isOutlier(QC$all$log10_total_features_by_counts, type="lower", nmad=3) 
high.spike <- isOutlier(QC$feature_control_ERCC$pct_counts, type="higher", nmad=3)
low.spike <- isOutlier(QC$feature_control_ERCC$log10_total_counts, type="lower", nmad=2)
data.frame(LowLib=sum(low.lib), LowNgenes=sum(low.genes), 
           HighSpike=sum(high.spike, na.rm=TRUE), LowSpike=sum(low.spike))

#  LowLib LowNgenes HighSpike LowSpike
#    162       572       904      359

# 舍弃低质量细胞
discard <- low.lib | low.genes | high.spike | low.spike
sce.e5601 <- sce.e5601[,!discard]
> summary(discard)
   Mode   FALSE    TRUE 
logical    2285    1229 

# 聚类
clusters <- quickCluster(sce.e5601, min.mean=1, method="igraph")
> table(clusters)
clusters
  1   2   3   4   5   6 
305 307 469 272 494 438 

# 标准化
sce.e5601 <- computeSumFactors(sce.e5601, min.mean=1, clusters=clusters)
sce.e5601 <- computeSpikeFactors(sce.e5601, general.use=FALSE)
sce.e5601 <- normalize(sce.e5601)

因为这个数据中donor信息比较多，所以可视化也要特别对待

donors <- sort(unique(sce.e5601$Donor))
> donors
 [1] "AZ"           "HP1502401"    "HP1504101T2D" "HP1504901"   
 [5] "HP1506401"    "HP1507101"    "HP1508501T2D" "HP1509101"   
 [9] "HP1525301T2D" "HP1526901T2D"

一共10个donor，作图可以设置这个参数，调整图片为2列

par(mfrow=c(ceiling(length(donors)/2), 2), 
    mar=c(4.1, 4.1, 2.1, 0.1))

代码作图，注意这段代码和之前的不同

collected <- list() # 第一行可以先不管，目的是创建一个空列表
# 下面进行一个循环，对10个donor进行循环：先取出第一个donor的列信息，然后使用if判断它是不是大于两列（也就是说：这个donor是不是有两个以上的细胞样本），如果只有一列那么就舍去；然后对这个donor的所有列进行标准化，去掉细胞文库差异；接着利用trendVar和decomposeVar鉴定HVGs，然后和之前一样进行可视化；最后将这个donor鉴定出来的HVGs信息放入collected这个列表中，留着以后用
for (x in unique(sce.e5601$Donor)) {
  current <- sce.e5601[,sce.e5601$Donor==x]
  if (ncol(current)<2L) { next }
  current <- normalize(current)
  fit <- trendVar(current, parametric=TRUE) 
  dec <- decomposeVar(current, fit)
  plot(dec$mean, dec$total, xlab="Mean log-expression",
       ylab="Variance of log-expression", pch=16, main=x)
  points(fit$mean, fit$var, col="red", pch=16)
  curve(fit$trend(x), col="dodgerblue", add=TRUE)
  collected[[x]] <- dec
}

因为这个数据中donor信息比较多，因此我们需要将不同donor的HVGs整合成一个数据框(注意是更高级的SV4数据框)

dec.e5601 <- do.call(combineVar, collected)
dec.e5601$Symbol <- rowData(sce.e5601)$Symbol
dec.e5601 <- dec.e5601[order(dec.e5601$bio, decreasing=TRUE),]
> head(dec.e5601,3)
DataFrame with 3 rows and 7 columns
                            mean            total              bio
                       <numeric>        <numeric>        <numeric>
ENSG00000115263 9.79547495804957 24.9059740209558  24.693297105741
ENSG00000118271 10.3601718361198 19.0590510324402 18.9670050741979
ENSG00000089199 8.78499018265489 17.2605488560106 16.9971950283286
                              tech   p.value       FDR      Symbol
                         <numeric> <numeric> <numeric> <character>
ENSG00000115263  0.212676915214769         0         0         GCG
ENSG00000118271 0.0920459582422512         0         0         TTR
ENSG00000089199  0.263353827682004         0         0        CHGB

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 195,783评论 5赞 462
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 82,360评论 2赞 373
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 142,942评论 0赞 325
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 52,507评论 1赞 267
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 61,324评论 5赞 358
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 46,299评论 1赞 273
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 36,685评论 3赞 386
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 35,358评论 0赞 254
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 39,652评论 1赞 293
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 34,704评论 2赞 312
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 36,465评论 1赞 326
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 32,318评论 3赞 313
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 37,711评论 3赞 299
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 28,991评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,265评论 1赞 251
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 41,661评论 2赞 342
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 40,864评论 2赞 335

单细胞转录组数据校正批次效应实战（中）

第二个数据--CEL-seq2, GSE85241

读数据，看数据

提取meta信息

提取基因、ERCC信息

基因转换

创建单细胞对象

质控和标准化

聚类

标准化

鉴定HVGs

第三个数据--Smart-seq2, E-MTAB-5061

读入数据

判断ERCC

基因转换

提取metadata信息

创建单细胞对象

后面的操作和之前保持一致了

推荐阅读更多精彩内容