公众号里的文章大多数需要编程基础，如果因为代码看不懂，而跟不上正文的节奏，可以来找我学习，相当于给自己一个新手保护期。我的课程都是循环开课，点进去咨询微信↓ 生信分析直播课程(每月初开一期，春节休一个月) 生信新手保护学习小组（每月两期） 单细胞陪伴学习小组（每月两期)
需求
我们的线上直播课每月会开一次直播答疑，既往报名的老学员都可以参与，不加钱不限时，主打一个关爱，如果你不知道的话现在知道了哦。
两个需求：
第一个是做拟时序的不同发育阶段(state)之间的差异分析，并且希望加上每个基因在当前state和其他state各自的基因表达量平均值。
第二个是想要一个平均表达量矩阵，每列是一个State，每行是一个基因。
说明一下：图是公司给的结果，里面的数值是未log的tpm，我们采用的是logNormlize之后的数据，更好。
学生说，跟公司要代码，人家不给。哈哈，那能给你才怪。所以来找我咯。
输入数据
首先是前面的monocle单样本代码，p2+p1/p3之前的行原封不动的运行，得到sc_cds，和原来的输入数据scRNA(seurat对象)一起存为Rdata。
拟时序分析完成后，CellDataSet对象里面已经有了state信息。
rm(list = ls())
library(Seurat)
## Loading required package: SeuratObject
## Loading required package: sp
##
## Attaching package: 'SeuratObject'
## The following objects are masked from 'package:base':
##
## intersect, t
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
load("sc_exp.Rdata")
head(sc_cds@phenoData@data) #state信息
## orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5
## AAACCCAAGTAGGGTC a 10768 3213 7.030089 7
## AAAGGTAGTCTTGAGT a 5565 1936 5.804133 7
## AACAAAGTCTTCCGTG a 8503 3063 7.397389 7
## AACCACATCCTTCACG a 12678 3846 6.830730 7
## AACGGGATCATTTCCA a 8032 2501 2.079183 1
## AACGTCATCTCGGCTT a 10596 3221 5.841827 1
## seurat_clusters celltype Size_Factor Pseudotime State
## AAACCCAAGTAGGGTC 7 FCGR3A+ Mono 1.5413473 0.6942694 1
## AAAGGTAGTCTTGAGT 7 FCGR3A+ Mono 0.7965823 0.6812000 1
## AACAAAGTCTTCCGTG 7 FCGR3A+ Mono 1.2171319 0.2854422 1
## AACCACATCCTTCACG 7 FCGR3A+ Mono 1.8147476 6.7191558 3
## AACGGGATCATTTCCA 1 CD14+ Mono 1.1497123 17.0878087 5
## AACGTCATCTCGGCTT 1 CD14+ Mono 1.5167270 17.6253461 4
检查sc_cds和scRNA中的细胞顺序是否一致，确认是一致的才可以直接在scRNA里添加sc_cds里面的列。
head(scRNA@meta.data)
## orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5
## AAACCCAAGTAGGGTC a 10768 3213 7.030089 7
## AAAGGTAGTCTTGAGT a 5565 1936 5.804133 7
## AACAAAGTCTTCCGTG a 8503 3063 7.397389 7
## AACCACATCCTTCACG a 12678 3846 6.830730 7
## AACGGGATCATTTCCA a 8032 2501 2.079183 1
## AACGTCATCTCGGCTT a 10596 3221 5.841827 1
## seurat_clusters celltype
## AAACCCAAGTAGGGTC 7 FCGR3A+ Mono
## AAAGGTAGTCTTGAGT 7 FCGR3A+ Mono
## AACAAAGTCTTCCGTG 7 FCGR3A+ Mono
## AACCACATCCTTCACG 7 FCGR3A+ Mono
## AACGGGATCATTTCCA 1 CD14+ Mono
## AACGTCATCTCGGCTT 1 CD14+ Mono
identical(rownames(scRNA@meta.data),rownames(sc_cds@phenoData@data))
## [1] TRUE
需求1：State间差异分析
#确认一致，可以添加啦
scRNA$State = sc_cds@phenoData@data$State
Idents(scRNA) = scRNA$State #设为scRNA的Idents
设置Idents为State
Idents是Seurat对象的重要内容之一，是细胞的分类。在常规的Seurat分析流程里，这个Idents是一直在变化的。
数据读取进来时，这个Idents是样本名称，如果是单样本数据那就只有一个分类，名字可以自己设置(CreateSeuratObject的project参数)
完成FindCluster步骤后，Idents会变成亚群编号（0，1，2，3…）
完成细胞类型注释后，Idents会被设置成细胞名称。
它的Idents是什么，FindAllmarkers就会找什么分类之间的差别。
我们现在要找State之间的差异基因，就设置为State啦。
计算差异基因
用monocle自带的差异分析函数得出的结果信息不太齐全，还是seurat好啊。
scRNA.markers <- FindAllMarkers(scRNA,
only.pos = TRUE,
min.pct = 0.25)
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
scRNA.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
## # A tibble: 10 × 7
## # Groups: cluster [5]
## p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
## 1 2.02e-10 3.60 0.396 0.048 5.25e- 6 1 ARID5B
## 2 9.16e- 9 3.95 0.302 0.027 2.38e- 4 1 P2RY10
## 3 7.48e-16 8.33 0.333 0 1.94e-11 2 ASPM
## 4 7.48e-16 8.33 0.333 0 1.94e-11 2 SYT9
## 5 1.68e- 8 2.44 0.549 0.168 4.37e- 4 3 PRDM1
## 6 2.91e- 5 2.79 0.353 0.101 7.57e- 1 3 C1QA
## 7 3.14e- 6 2.43 0.32 0.067 8.16e- 2 4 LINC01094
## 8 1.36e- 3 2.52 0.26 0.093 1 e+ 0 4 CD9
## 9 4.91e- 6 4.20 0.302 0.064 1.28e- 1 5 NAF1
## 10 2.03e- 5 4.03 0.395 0.134 5.27e- 1 5 HSPA6
因为将State设置为了Idents，这里的”cluster”列其实是State。把列名改掉，避免歧义。
colnames(scRNA.markers)[6] = "State"
计算平均表达量
即计算每个基因在当前state和其他state的平均表达量，学生说为了方便直接查看和比较表达量。
exp = scRNA@assays$RNA$data
a = apply(scRNA.markers,1,function(x){
c(mean(exp[x[7],scRNA$State==x[6]]),
mean(exp[x[7],scRNA$State!=x[6]]))
})
a = as.data.frame(t(a))
colnames(a) = c("Target_State","Other_State")
head(a)
## Target_State Other_State
## VMO1 2.468474 0.3722993
## TNFSF13B 1.709487 0.2883641
## CDKN1C 1.612884 0.2331031
## JPT1 2.585694 1.0113633
## CYTIP 2.436601 1.0251027
## PAPSS2 1.812296 0.5325872
scRNA.markers = cbind(scRNA.markers,a)
head(scRNA.markers)
## p_val avg_log2FC pct.1 pct.2 p_val_adj State gene
## VMO1 1.285754e-27 3.447251 1.000 0.245 3.341417e-23 1 VMO1
## TNFSF13B 8.299190e-23 3.268468 0.925 0.245 2.156793e-18 1 TNFSF13B
## CDKN1C 1.761194e-21 3.300786 0.868 0.163 4.576990e-17 1 CDKN1C
## JPT1 2.570100e-21 2.154897 1.000 0.646 6.679175e-17 1 JPT1
## CYTIP 1.086084e-18 1.783399 1.000 0.626 2.822514e-14 1 CYTIP
## PAPSS2 1.802407e-18 2.258809 0.962 0.408 4.684095e-14 1 PAPSS2
## Target_State Other_State
## VMO1 2.468474 0.3722993
## TNFSF13B 1.709487 0.2883641
## CDKN1C 1.612884 0.2331031
## JPT1 2.585694 1.0113633
## CYTIP 2.436601 1.0251027
## PAPSS2 1.812296 0.5325872
搞定啦。
需求2：state平均表达量矩阵
我本以为应该要很多行代码才能搞定。结果只要两行！chatgpt写了十几行还报错。还是自己写的好使嗯。
我写的哦👇。如有雷同是别人抄我的。
计算！
dat = split(colnames(exp),scRNA$State)
re = sapply(dat,function(x){rowMeans(exp[,x])})
head(re)
## 1 2 3 4 5
## TSPAN6 0.00000000 0.0000000 0.0000000 0.0000000 0.0000000
## DPM1 0.20975283 0.5425002 0.2493284 0.3039274 0.2149655
## SCYL3 0.03053863 0.0000000 0.0831260 0.0000000 0.0000000
## C1orf112 0.01484739 0.0000000 0.0000000 0.0000000 0.0000000
## FGR 0.69361333 0.3610408 0.9225549 1.1443649 0.9587395
## CFH 0.00000000 0.0000000 0.0000000 0.0000000 0.0000000
让chatgpt讲讲这两句代码的意思，超多料。让它写他不一定行，让它解释肯定是行的！
第一行代码：dat = split(colnames(exp), scRNA$State)
colnames(exp)：这部分获取表达矩阵exp的列名，通常列名代表不同的样本或细胞。split(colnames(exp), scRNA$State)：split函数用于根据scRNA$State中的每个细胞的状态将列名（样本/细胞）分组。scRNA$State是一个向量，其中包含了与exp矩阵中每一列相对应的状态信息（例如，不同的细胞类型或条件）。dat：这个变量将存储split函数的结果，即一个列表，其中每个元素包含了对应于相同状态的所有列名。
第二行代码：re = sapply(dat, function(x) {rowMeans(exp[,x])})
sapply(dat, function(x) {rowMeans(exp[,x])})：sapply函数用于对dat列表中的每个元素（即每个状态组的列名）应用一个函数。这里应用的函数是rowMeans，它计算exp矩阵中每一行（基因）在指定列（细胞/样本）上的平均值。function(x) {rowMeans(exp[,x])}：这是一个匿名函数，它接受一个向量x（代表一组列名），并计算exp矩阵中这些列的平均表达量。re：这个变量将存储sapply函数的结果，即一个向量，其中包含了每个状态组的基因平均表达量。
总结
这两行代码的目的是将表达矩阵exp根据细胞状态分组，并计算每个组中所有基因的平均表达量。这样，你可以比较不同状态（如不同细胞类型或不同疾病状态）下基因表达的差异。这种方法在单细胞RNA测序数据分析中非常常见，用于识别状态特异性的基因表达模式。
检查!
验证一个，看上面的代码是否正确啊。
mean(exp["DPM1",scRNA$State==1])
## [1] 0.2097528
re["DPM1",1]
## [1] 0.2097528
搞定啦。
