这里是佳奥，我们继续补充相关的背景知识。
1 细胞周期推断
也是通过一些基因的表达来进行分类的。
rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input.Rdata')
a[1:4,1:4]
head(df)
##载入第0步准备好的表达矩阵，及细胞的一些属性（hclust分群，plate批次，检测到的基因数量）
##注意 变量a是原始的counts矩阵，变量 dat是logCPM后的表达量矩阵。
group_list=df$g
plate=df$plate
table(plate)
a[1:4,1:4]
library(scran)
# https://mp.weixin.qq.com/s/nFSa5hXuKHrGu_othopbWQ
sce <- SingleCellExperiment(list(counts=dat))
#list() 创建列表
##这个是小鼠的
library(org.Mm.eg.db)
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds",
package="scran"))
ensembl <- mapIds(org.Mm.eg.db, keys=rownames(sce),
keytype="SYMBOL", column="ENSEMBL")
#取探针名创建一个向量
#rownames(sce) 取行名（即实验检测到的基因）
if(F){
assigned <- cyclone(sce, pairs=mm.pairs, gene.names=ensembl)
save(assigned,file = 'cell_cycle_assigned.Rdata')
}
load(file = 'cell_cycle_assigned.Rdata')
head(assigned$scores)
table(assigned$phases)
draw=cbind(assigned$score,assigned$phases) #合并assigned$score列和assigned$phases列
colnames(draw)
attach(draw)
library(scatterplot3d)
scatterplot3d(G1, S, G2M, angle=20,
color = rainbow(3)[as.numeric(as.factor(assigned$phases))],
grid=TRUE, box=FALSE)
detach(draw)
library(pheatmap)
cg=names(tail(sort(apply(dat,1,sd)),100))
n=t(scale(t(dat[cg,])))
# pheatmap(n,show_colnames =F,show_rownames = F)
library(pheatmap)
df$cellcycle=assigned$phases
ac=df
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,
filename = 'all_cells_top_100_sd_all_infor.png')
dev.off()
head(ac)
table(ac[,c(1,5)])
> table(ac[,c(1,5)])
cellcycle
g G1 G2M S
1 298 10 4
2 272 22 6
3 121 0 0
4 32 2 1
QQ截图20220901102109.png
QQ截图20220901102213.png
2 rpkm概念与三种计算方法
这一部分讲解rpkm值如何在算法上与counts值统一。
即如何从5 Mb大小的counts矩阵得到23 Mb大小的rpkm矩阵：
rm(list = ls())
options(stringsAsFactors = F)
load(file = '../input.Rdata')
a[1:4,1:4]
head(df)
##载入第0步准备好的表达矩阵，及细胞的一些属性（hclust分群，plate批次，检测到的基因数量）
##注意 变量a是原始的counts矩阵，变量 dat是log2CPM后的表达量矩阵。
library("TxDb.Mmusculus.UCSC.mm10.knownGene")
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
?TxDb
##part1 下面是定义基因长度为 非冗余exon长度之和
if(F){
exon_txdb=exons(txdb)
genes_txdb=genes(txdb)
genes_txdb
?GRanges
o = findOverlaps(exon_txdb,genes_txdb)
o
## exon - 1 : chr1 4807893-4807982
## 1 6523
# genes_txdb[6523] # chr1 4807893-4846735 , 18777
t1=exon_txdb[queryHits(o)]
t2=genes_txdb[subjectHits(o)]
t1=as.data.frame(t1)
t1$geneid=mcols(t2)[,1]
# 如果觉得速度不够，就参考R语言实现并行计算
# http://www.bio-info-trainee.com/956.html
# lapply : 遍历列表向量内的每个元素，并且使用指定函数来对其元素进行处理。返回列表向量。
# 函数split()可以按照分组因子，把向量，矩阵和数据框进行适当的分组;
# 它的返回值是一个列表，代表分组变量每个水平的观测。
g_l = lapply(split(t1,t1$geneid),function(x){
# x=split(t1,t1$geneid)[[1]]
head(x)
tmp=apply(x,1,function(y){
y[2]:y[3]
})
length(unique(unlist(tmp)))
# sum(x[,4])
})
head(g_l)
g_l=data.frame(gene_id=names(g_l),length=as.numeric(g_l))
save(g_l,file = 'step7-g_l.Rdata')
}
load(file = 'step7-g_l.Rdata')
##part2 下面是定义基因长度为 最长转录本长度
if(F){
t_l=transcriptLengths(txdb)
head(t_l)
t_l=na.omit(t_l)
head(t_l)
t_l=t_l[order(t_l$gene_id,t_l$tx_len,decreasing = T),]
head(t_l)
str(t_l)
t_l=t_l[!duplicated(t_l$gene_id),]
head(t_l)
g_l=t_l[,c(3,5)]
}
head(g_l)
library(org.Mm.eg.db)
s2g=toTable(org.Mm.egSYMBOL)
head(s2g)
g_l=merge(g_l,s2g,by='gene_id') #把g_l,s2g两个数据框以'gene_id'为连接进行拼接
#merge函数可以实现对两个数据表进行匹配和拼接的功能。
# 参考counts2rpkm，定义基因长度为非冗余CDS之和
# http://www.bio-info-trainee.com/3298.html
a[1:4,1:4]
ng=intersect(rownames(a),g_l$symbol) #取a数据框的行名与g_l数据框的symbol列的交集
#intersect()取交集
# 有了counts矩阵和对应的基因长度信息，就很容易进行各种计算了：
exprSet=a[ng,]
lengths=g_l[match(ng,g_l$symbol),2]
head(lengths)
head(rownames(exprSet))
# http://www.biotrainee.com/thread-1791-1-1.html
exprSet[1:4,1:4]
total_count<- colSums(exprSet)
head(total_count)
head(lengths)
total_count[4]
lengths[1]
1*10^9/(1122*121297)
rpkm <- t(do.call( rbind,
lapply(1:length(total_count),
function(i){
10^9*exprSet[,i]/lengths/total_count[i]
}) ))
rpkm[1:4,1:4]
# 下面可以比较一下 自己根据counts值算出来的RPKM和作者提供的RPKM区别。
a=read.table('../GSE111229_Mammary_Tumor_fibroblasts_768samples_rpkmNormalized.txt.gz',header = T ,sep = '\t')
# 每次都要检测数据
a[1:4,1:4]
rpkm_paper=a[ng,]
rpkm_paper[1:4,1:4]
rpkm[1:4,1:4]
有关背景讲解就到这里。
下一篇我们进入单细胞转录组数据初探部分。
我们下一篇再见！
