文章/答案/技术大牛

发布

社区首页 >专栏 >拟时序分析的State表达矩阵和差异基因

拟时序分析的State表达矩阵和差异基因

用户11414625

发布于 2024-12-20 08:43:26

17400

代码可运行

文章被收录于专栏：生信星球520生信星球520

运行总次数：0

代码可运行

需求

我们的线上直播课每月会开一次直播答疑，既往报名的老学员都可以参与，不加钱不限时，主打一个关爱，如果你不知道的话现在知道了哦。

两个需求：

第一个是做拟时序的不同发育阶段(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), scRNAState)：split函数用于根据scRNAState中的每个细胞的状态将列名（样本/细胞）分组。scRNA
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函数的结果，即一个向量，其中包含了每个状态组的基因平均表达量。