概要
本文主要討論Seurat對(duì)象導(dǎo)入到Monocle中直接進(jìn)行分析的可行性栓袖，分兩種情況：
①經(jīng)過數(shù)據(jù)清洗叽赊、標(biāo)準(zhǔn)化和聚類的Seurat對(duì)象導(dǎo)入
②未經(jīng)過任何處理的Seurat對(duì)象導(dǎo)入
以下先進(jìn)行Monocle包的簡單介紹，再分這兩種情況進(jìn)行嘗試塔橡。
為什么要分這兩種情況進(jìn)行嘗試？
- Seurat包中也有將數(shù)據(jù)標(biāo)準(zhǔn)化的步驟癞谒，作者的建議是在Monocle中要再次進(jìn)行標(biāo)準(zhǔn)化弹砚，但是他自己也沒有嘗試過桌吃，所以不確定會(huì)怎么樣逗物。
- Seurat包中有個(gè)
ScaleData的命令翎卓，目的是去除測序產(chǎn)生的批次效應(yīng)和技術(shù)噪音，但對(duì)于我們的數(shù)據(jù)（按不同時(shí)間缺血處理的脾臟锐帜，根據(jù)錐蟲感染小鼠的時(shí)間進(jìn)行測序）畜号，我們要觀察的就是這些不同時(shí)間批次之間的差別缴阎，有可能這個(gè)命令會(huì)將這個(gè)差別掩蓋了。因此如果直接輸入已經(jīng)聚類好的Seurat對(duì)象简软，也許會(huì)出現(xiàn)問題蛮拔。
關(guān)于Monocle
http://cole-trapnell-lab.github.io/monocle-release/
【Introduction】
Monocle介紹了使用RNA-Seq進(jìn)行單細(xì)胞軌跡分析的策略，能夠?qū)⒓?xì)胞按照模擬的時(shí)間順序進(jìn)行排列痹升，顯示它們的發(fā)展軌跡如細(xì)胞分化等生物學(xué)過程。Monocle從數(shù)據(jù)中用無監(jiān)督或半監(jiān)督學(xué)習(xí)獲得這個(gè)軌跡疼蛾。
無監(jiān)督：利用Monocle的自己一套工具或Seurat生成一個(gè)基因列表
半監(jiān)督：通過自身的知識(shí)積累人為輸入一些認(rèn)為重要的基因
Monocle不是通過實(shí)驗(yàn)將細(xì)胞純化為離散狀態(tài)肛跌，而是使用算法來學(xué)習(xí)每個(gè)細(xì)胞必須經(jīng)歷的基因表達(dá)變化的序列，作為動(dòng)態(tài)生物過程的一部分察郁。一旦它了解了基因表達(dá)變化的整體“軌跡”衍慎，Monocle就可以將每個(gè)細(xì)胞放置在軌跡中的適當(dāng)位置。然后皮钠，可以使用Monocle的差異分析工具包來查找在軌跡過程中受到調(diào)控的基因稳捆。如果該過程有多個(gè)結(jié)果，Monocle將重建“分支”軌跡麦轰。這些分支對(duì)應(yīng)于細(xì)胞“決策”，Monocle提供了強(qiáng)大的工具來識(shí)別受其影響的基因并參與制作它們侧纯。網(wǎng)站中也提供了分析分支的方法鞠评。Monocle依靠Reversed Graph Embedding的機(jī)器學(xué)習(xí)技術(shù)來構(gòu)建單細(xì)胞軌跡茂蚓。
除了構(gòu)建單細(xì)胞軌跡之外壕鹉，它還能夠做差異表達(dá)分析和聚類來揭示重要的基因和細(xì)胞剃幌。這與Seurat的功能相似。
【W(wǎng)orkflow以及與Seurat的異同】
①創(chuàng)建CellDataSet對(duì)象（下簡稱CDS對(duì)象）
若要?jiǎng)?chuàng)建新的CDS對(duì)象晾浴，則需要整理出3個(gè)輸入文件（基因-細(xì)胞表達(dá)矩陣负乡、細(xì)胞-細(xì)胞特征注釋矩陣、基因-基因特征注釋矩陣）脊凰，但方便的是抖棘，Monocle支持從Seurat中直接導(dǎo)入對(duì)象，通過importCDS命令實(shí)現(xiàn)狸涌。
在創(chuàng)建之后切省，也會(huì)進(jìn)行數(shù)據(jù)過濾和標(biāo)準(zhǔn)化，不同的是Seurat是基于作圖的方法去篩選掉異常的細(xì)胞帕胆，而Monocle的過濾方法則是根據(jù)表達(dá)數(shù)據(jù)的數(shù)學(xué)關(guān)系來實(shí)現(xiàn)朝捆。
上限：
下限：
其中為基因表達(dá)總數(shù), n 為細(xì)胞數(shù)，
為表達(dá)水平方差
大概就是以一個(gè)大致的細(xì)胞表達(dá)水平為基準(zhǔn)懒豹，表達(dá)量太高的跟太低的去除掉芙盘。②通過已知的Marker基因分類細(xì)胞
方法一：通過一些現(xiàn)有的生物/醫(yī)學(xué)知識(shí)手動(dòng)賦予它們細(xì)胞名，將這些細(xì)胞分為幾大類脸秽，然后再聚類細(xì)胞儒老。
方法二：與Seurat包的分類細(xì)胞方法類似，篩選出差異表達(dá)基因用于聚類记餐，然后再根據(jù)每一個(gè)cluster的marker賦予細(xì)胞名驮樊。③聚類細(xì)胞
采用的也是t-SNE算法。-
④將細(xì)胞按照偽時(shí)間的順序排列在軌跡上
-
Step1：選擇輸入基因用于機(jī)器學(xué)習(xí)
這個(gè)過程稱為feature selection（特征選擇）片酝，這些基因?qū)壽E的形狀有著最重要的影響巩剖。我們應(yīng)該要選擇的是最能反映細(xì)胞狀態(tài)的基因。
如果直接通過Seurat輸出的一些重要的基因（比如每個(gè)cluster中的高FC值基因）作為輸入對(duì)象的話就能夠?qū)崿F(xiàn)一個(gè)“無監(jiān)督”分析钠怯〖涯В或者也可以利用生物學(xué)知識(shí)手動(dòng)選擇一些重要的基因進(jìn)行“半監(jiān)督”分析。 -
Step2：數(shù)據(jù)降維
利用Reversed graph embedding算法將數(shù)據(jù)降維晦炊。
相對(duì)于PCA來說這個(gè)算法更能夠反應(yīng)數(shù)據(jù)的內(nèi)部結(jié)構(gòu)（據(jù)monocle網(wǎng)站說是這樣） -
Step3：將細(xì)胞按照偽時(shí)間排序
這個(gè)過程稱為manifold learning（流形學(xué)習(xí)）鞠鲜。Monocle利用軌跡來描述細(xì)胞如何從一個(gè)狀態(tài)轉(zhuǎn)換到另一個(gè)狀態(tài)宁脊。得到的是一個(gè)樹狀圖，樹的根部或樹干表示的是細(xì)胞最初的狀態(tài)（如果有的話）贤姆，隨著細(xì)胞的變化就沿著分叉樹枝發(fā)展榆苞。一個(gè)細(xì)胞的偽時(shí)間值（pseudotime value）為它的位置沿著樹枝到根部的距離。
-
Step1：選擇輸入基因用于機(jī)器學(xué)習(xí)
⑤差異表達(dá)分析
還沒細(xì)看
情況①：經(jīng)過清洗霞捡、標(biāo)準(zhǔn)化和聚類的Seurat對(duì)象導(dǎo)入
spleen
An object of class seurat in project 10X_spleen
15655 genes across 1940 samples.
clustered_spleen_monocle <- importCDS(spleen, import_all = TRUE)
head(pData(clustered_spleen_monocle))
nGene nUMI orig.ident percent.mito res.0.6 Size_Factor
AAACCTGAGGTGATAT 1072 3999 10X_spleen 0.01150863 3 NA
AAACCTGAGTCATCCA 1562 4640 10X_spleen 0.03709295 5 NA
AAACCTGCAGTGAGTG 995 3489 10X_spleen 0.03955288 3 NA
AAACGGGAGACTGGGT 1704 7397 10X_spleen 0.02852508 0 NA
AAACGGGAGACTTGAA 976 3102 10X_spleen 0.04932302 2 NA
AAACGGGAGATAGGAG 1236 4548 10X_spleen 0.02682498 1 NA
res.0.6為cluster的編號(hào)坐漏，將此列名更改為“Cluster”，方便后面使用碧信。
names(pData(clustered_spleen_monocle))[names(pData(clustered_spleen_monocle))=="res.0.6"]="Cluster"
range(pData(clustered_spleen_monocle)$Cluster)
[1] "0" "8"
確認(rèn)此列的范圍為0到8赊琳，也就是共9個(gè)cluster。
計(jì)算SizeFactor和Dispersions
計(jì)算用于方便之后的分析使用砰碴。
clustered_spleen_monocle <- estimateSizeFactors(clustered_spleen_monocle)
clustered_spleen_monocle <- estimateDispersions(clustered_spleen_monocle)
跳過-數(shù)據(jù)清洗
由于此數(shù)據(jù)已經(jīng)經(jīng)過數(shù)據(jù)清洗躏筏，因此沒有必要在Monocle中再一次處理。
數(shù)據(jù)標(biāo)準(zhǔn)化
根據(jù)作者的建議呈枉，即使在Seurat包中已經(jīng)標(biāo)準(zhǔn)化處理過的數(shù)據(jù)趁尼，在轉(zhuǎn)化到Monocle中時(shí)仍然需要再一次進(jìn)行標(biāo)準(zhǔn)化。
#將表達(dá)矩陣中所有值進(jìn)行l(wèi)og標(biāo)準(zhǔn)化
L <- log(exprs(clustered_spleen_monocle[expressed_genes,]))
#將每個(gè)基因都標(biāo)準(zhǔn)化猖辫，melt方便作圖
library(reshape)
melted_dens_df <- melt(Matrix::t(scale(Matrix::t(L))))
#作圖酥泞，查看標(biāo)準(zhǔn)化的基因表達(dá)值的分布
qplot(value, geom = "density", data = melted_dens_df) +
+ stat_function(fun = dnorm, size = 0.5, color = 'red') +
+ xlab("Standardized log(FPKM)") +
+ ylab("Density")
利用細(xì)胞Marker分類細(xì)胞
首先，因?yàn)橐环N細(xì)胞下又可以細(xì)分成更加小的類別啃憎，因此在用marker給細(xì)胞類時(shí)芝囤，就要考慮到它們的對(duì)應(yīng)關(guān)系。如CD4基因所對(duì)應(yīng)的細(xì)胞為CD4+ T細(xì)胞荧飞，而CD4+T細(xì)胞屬于T細(xì)胞中的一種凡人，我們就要告訴Monocle CD4+T細(xì)胞屬于T細(xì)胞的子集，讓它不要再分類的過程中將它們分成兩類叹阔。
Monocle提供了一個(gè)將細(xì)胞分級(jí)的函數(shù)newCellTypeHierarchy.
cth <- newCellTypeHierarchy()
將marker和細(xì)胞對(duì)應(yīng)起來挠轴，以及排列好它們的從屬關(guān)系。
cth <- addCellType(cth, "T cell",
classify_func=function(x) {x["CD3D",] > 0})
cth <- addCellType(cth, "CD4+ T cell",
classify_func=function(x) {x["CD4",] > 0},
parent_cell_type_name = "T cell")
cth <- addCellType(cth, "CD8+ T cell",
classify_func=function(x) {
x["CD8A",] > 0 | x["CD8B",] > 0
},
parent_cell_type_name = "T cell")
cth <- addCellType(cth, "B cell",
classify_func=function(x) {x["MS4A1",] > 0})
cth <- addCellType(cth, "Monocyte",
classify_func=function(x) {x["CD14",] > 0})
cth <- addCellType(cth, "Neutrophil",
classify_func=function(x) {x["S100A9",] > 0})
下面這一步將細(xì)胞進(jìn)行分類耳幢。
clustered_spleen_monocle <- classifyCells(clustered_spleen_monocle, cth, 0.1)
查看細(xì)胞分類的情況岸晦。
table(pData(clustered_spleen_monocle)$CellType)
偽時(shí)間分析
Step1: 選擇定義細(xì)胞發(fā)展的基因
diff_test_res <- differentialGeneTest(clustered_spleen_monocle,fullModelFormulaStr = "~Cluster")
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
Step2: 數(shù)據(jù)集降維
clustered_spleen_monocle <- setOrderingFilter(clustered_spleen_monocle, ordering_genes)
plot_ordering_genes(clustered_spleen_monocle)
clustered_spleen_monocle <- reduceDimension(clustered_spleen_monocle, max_components = 2, reduction_method = "DDRTree")
Step3: 將細(xì)胞按照偽時(shí)間排序
clustered_spleen_monocle <- orderCells(clustered_spleen_monocle)
查看能用于上色區(qū)分的變量：
colnames(pData(clustered_spleen_monocle))
[1] "nGene" "nUMI" "orig.ident" "percent.mito" "Cluster"
[6] "Size_Factor" "CellType" "Pseudotime" "State"
其實(shí)這一步按照正常情況下來說，是應(yīng)該最好有一個(gè)時(shí)間的變量（如Hour或者Time）睛藻，用來區(qū)別不同時(shí)間批次處理產(chǎn)生的數(shù)據(jù)启上，子亮部分的數(shù)據(jù)就可以根據(jù)不同的寄生時(shí)間來給細(xì)胞上色，從而觀察隨著寄生時(shí)間的變化細(xì)胞狀態(tài)（發(fā)育/分化）的變化店印。而像脾臟數(shù)據(jù)這一種冈在，雖然按照了4個(gè)時(shí)間點(diǎn)進(jìn)行了處理，但是數(shù)據(jù)并沒有按照不同時(shí)間點(diǎn)區(qū)分出來按摘，因此只能夠根據(jù)細(xì)胞的分化的過程來確定哪些是原始狀態(tài)包券。
plot_cell_trajectory(clustered_spleen_monocle, color_by = "Cluster")
plot_cell_trajectory(clustered_spleen_monocle, color_by = "CellType")
plot_cell_trajectory(clustered_spleen_monocle, color_by = "State")
這是一種樹狀圖纫谅，有三條細(xì)胞軌跡，表示細(xì)胞狀態(tài)主要分為3個(gè)階段溅固，中間的數(shù)字1表示一個(gè)分叉付秕。
上圖的細(xì)胞依據(jù)不同的Cluster進(jìn)行上色，根據(jù)之前的Seurat聚類分析侍郭，Cluster5（淺藍(lán)色）對(duì)應(yīng)的是中性粒細(xì)胞询吴，在此圖位于上面那一分支的頂端；Cluster0（紅色）對(duì)應(yīng)的是B細(xì)胞亮元，主要位于右邊一支的最頂端猛计；左下角頂端的藍(lán)色有可能是NK細(xì)胞，但不確定苹粟。這樣看來比較適合做初始狀態(tài)的是右邊那一支有滑。與下圖對(duì)比得到的結(jié)果也差不多跃闹。
plot_cell_trajectory(clustered_spleen_monocle, color_by = "CellType") + facet_wrap(~CellType, nrow = 1)
上圖將每個(gè)細(xì)胞的分布軌跡分開顯示嵌削，比較明顯地看到B細(xì)胞集中的位置在右支頂端，然后集中為T細(xì)胞望艺，中間摻雜一些中性粒細(xì)胞（也有可能沒分好）苛秕。但是大部分的細(xì)胞還是沒有分出來，這個(gè)結(jié)果還需要再處理一下找默。
由于Monocle不能分辨哪一條軌跡才是“樹根”艇劫，也就是不知道哪一個(gè)細(xì)胞狀態(tài)才是更初始的，可設(shè)定root_state參數(shù)來設(shè)定哪條軌跡才是初始狀態(tài)惩激。然后賦予每一個(gè)細(xì)胞偽時(shí)間值店煞，就可以觀察基因在偽時(shí)間里的表達(dá)變化。待處理好細(xì)胞分類之后风钻，就可以接著做這一步顷蟀。
head(pData(clustered_spleen_monocle))
nGene nUMI orig.ident percent.mito Cluster Size_Factor CellType
AAACCTGAGGTGATAT 1072 3999 10X_spleen 0.01150863 3 0.8327676 T cell
AAACCTGAGTCATCCA 1562 4640 10X_spleen 0.03709295 5 0.9661104 Ambiguous
AAACCTGCAGTGAGTG 995 3489 10X_spleen 0.03955288 3 0.7269267 Monocyte
AAACGGGAGACTGGGT 1704 7397 10X_spleen 0.02852508 0 1.5411513 Ambiguous
AAACGGGAGACTTGAA 976 3102 10X_spleen 0.04932302 2 0.6462960 B cell
AAACGGGAGATAGGAG 1236 4548 10X_spleen 0.02682498 1 0.9475674 Ambiguous
情況②：未經(jīng)過標(biāo)準(zhǔn)化處理的Seurat對(duì)象導(dǎo)入
創(chuàng)建CellDataSet對(duì)象（下簡稱CDS對(duì)象）
創(chuàng)建Seurat對(duì)象spleen_monocle，先去除一些測序質(zhì)量差的細(xì)胞：
留下所有在>=3個(gè)細(xì)胞中表達(dá)的基因min.cells = 3骡技；
留下所有檢測到>=200個(gè)基因的細(xì)胞min.genes = 200鸣个。
spleen_monocle <- CreateSeuratObject(raw.data = spleen.data, min.cells = 3, min.genes = 200, project = "10X_spleen")
從Monocle中導(dǎo)入Seurat對(duì)象
library(monocle)
spleen_monocle <- importCDS(spleen_monocle, import_all = TRUE)
查看數(shù)據(jù)：
spleen_monocle
CellDataSet (storageMode: environment)
assayData: 15655 features, 1959 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGAGGTGATAT AAACCTGAGTCATCCA ... TTTGTCAGTCGTCTTC (1959 total)
varLabels: nGene nUMI orig.ident Size_Factor
varMetadata: labelDescription
featureData
featureNames: RP11-34P13.7 FO538757.2 ... AC240274.1 (15655 total)
fvarLabels: gene_short_name
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
head(fData(spleen_monocle))
gene_short_name
RP11-34P13.7 RP11-34P13.7
FO538757.2 FO538757.2
AP006222.2 AP006222.2
RP4-669L17.10 RP4-669L17.10
RP11-206L10.9 RP11-206L10.9
LINC00115 LINC00115
head(pData(spleen_monocle))
nGene nUMI orig.ident Size_Factor
AAACCTGAGGTGATAT 1072 3999 10X_spleen NA
AAACCTGAGTCATCCA 1562 4640 10X_spleen NA
AAACCTGCAGTGAGTG 995 3489 10X_spleen NA
AAACGGGAGACTGGGT 1704 7397 10X_spleen NA
AAACGGGAGACTTGAA 976 3102 10X_spleen NA
AAACGGGAGATAGGAG 1236 4548 10X_spleen NA
15655個(gè)基因，1959個(gè)細(xì)胞布朦，與之前創(chuàng)建的的Seurat對(duì)象相符囤萤。
計(jì)算SizeFactor和Dispersions
計(jì)算用于方便之后的分析使用。
spleen_monocle <- estimateSizeFactors(spleen_monocle)
spleen_monocle <- estimateDispersions(spleen_monocle)
數(shù)據(jù)清洗
根據(jù)前文所說的表達(dá)式是趴，可以利用nUMI值進(jìn)行過濾
head(pData(spleen_monocle))
nGene nUMI orig.ident Size_Factor num_genes_expressed
AAACCTGAGGTGATAT 1072 3999 10X_spleen 0.8207242 1070
AAACCTGAGTCATCCA 1562 4640 10X_spleen 0.9521386 1560
AAACCTGCAGTGAGTG 995 3489 10X_spleen 0.7164140 995
AAACGGGAGACTGGGT 1704 7397 10X_spleen 1.5188634 1704
AAACGGGAGACTTGAA 976 3102 10X_spleen 0.6369493 976
AAACGGGAGATAGGAG 1236 4548 10X_spleen 0.9338638 1236
觀察到這里多了一個(gè)Size_Factor的列
#設(shè)置上下限：
upper_bound_spleen <- 10^(mean(log10(pData(spleen_monocle)$nUMI))
+ 2*sd(log10(pData(spleen_monocle)$nUMI)))
lower_bound_spleen <- 10^(mean(log10(pData(spleen_monocle)$nUMI))
- 2*sd(log10(pData(spleen_monocle)$nUMI)))
#可視化
qplot(nUMI,data = pData(spleen_monocle), geom = "density") + geom_vline(xintercept = lower_bound_spleen) + geom_vline(xintercept = upper_bound_spleen)
留下兩條豎線中間的部分：
spleen_monocle <- spleen_monocle[,pData(spleen_monocle)$nUMI > lower_bound_spleen & pData(spleen_monocle)$nUMI < upper_bound_spleen]
spleen_monocle
CellDataSet (storageMode: environment)
assayData: 15655 features, 1864 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGAGGTGATAT AAACCTGAGTCATCCA ... TTTGTCAGTCGTCTTC (1864
total)
varLabels: nGene nUMI orig.ident Size_Factor
varMetadata: labelDescription
featureData
featureNames: RP11-34P13.7 FO538757.2 ... AC240274.1 (15655 total)
fvarLabels: gene_short_name
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
過濾后剩下1864個(gè)細(xì)胞涛舍，15655個(gè)基因
數(shù)據(jù)標(biāo)準(zhǔn)化
# 將表達(dá)矩陣中所有值進(jìn)行l(wèi)ognormalize.
L <- log(exprs(spleen_monocle[expressed_genes,]))
# Standardize each gene, so that they are all on the same scale,
# Then melt the data with plyr so we can plot it easily
library("reshape")
melted_dens_df <- melt(Matrix::t(scale(Matrix::t(L))))
# Plot the distribution of the standardized gene expression values.
qplot(value, geom = "density", data = melted_dens_df) +
stat_function(fun = dnorm, size = 0.5, color = 'red') +
xlab("Standardized log(FPKM)") +
ylab("Density")
分類細(xì)胞
disp_table <- dispersionTable(spleen_monocle)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
spleen_monocle <- setOrderingFilter(spleen_monocle, unsup_clustering_genes$gene_id)
plot_ordering_genes(spleen_monocle)
setOrderingFilter將一些用于之后聚類的基因標(biāo)記起來；
plot_ordering_genes根據(jù)這些基因的平均表達(dá)水平展示了基因的表達(dá)差異程度唆途，紅線表示Monocle基于這種關(guān)系對(duì)分散的期望富雅。我們用于聚類的基因顯示為黑點(diǎn)缤削，而其他基因顯示為灰點(diǎn)。（這里有點(diǎn)不太明白縱坐標(biāo)的dispersion經(jīng)驗(yàn)值是什么意思）
聚類細(xì)胞
作碎石圖：
plot_pc_variance_explained(spleen_monocle, return_all = F, max_components = 25)
選擇前8個(gè)成分進(jìn)行聚類
spleen_monocle <- reduceDimension(spleen_monocle, max_components = 2, num_dim = 8, reduction_method = 'tSNE', verbose = T)
spleen_monocle <- clusterCells(spleen_monocle)
如果說每個(gè)時(shí)間點(diǎn)的細(xì)胞都聚在一起的話吹榴，那么大致上看此圖正好分成4個(gè)模塊亭敢。
將細(xì)胞按照偽時(shí)間的順序排列在軌跡上
Step1:選擇定義細(xì)胞發(fā)展的基因
這里可以用Monocle包中的dpFeature命令進(jìn)行選擇基因。
#這個(gè)命令運(yùn)行時(shí)間特別長
diff_test_res <- differentialGeneTest(spleen_monocle,fullModelFormulaStr = "~Cluster")
ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
另一種方式是基于生物知識(shí)人為選擇基因：
HSPA1A*基因是之前用Seurat包發(fā)現(xiàn)的一個(gè)比較有意思的基因图筹，幾乎在所有的cluster中都有不同程度的表達(dá)（下圖）帅刀，并且查閱文獻(xiàn)發(fā)現(xiàn)此基因表達(dá)熱休克蛋白，在機(jī)體發(fā)生缺血/缺氧時(shí)此蛋白的表達(dá)是一種保護(hù)機(jī)制远剩，可以作為心臟驟涂勰纾患者的預(yù)后標(biāo)志物。還有學(xué)者研究基于腦缺血的持續(xù)時(shí)間跟HSPA1A表達(dá)量的關(guān)系
瓜晤，盡管該論文做了兩組（缺血時(shí)間30分鐘和60分鐘）的最后的結(jié)論是這兩組的表達(dá)量差異沒有顯著性锥余，但是脾臟缺血是經(jīng)過了12h, 24h 和72h的處理，時(shí)間跨度大很多痢掠，因此還是具有比較好的研究價(jià)值驱犹。
因此我認(rèn)為選擇HSPA1A作為反映細(xì)胞狀態(tài)的標(biāo)志基因也是可以的。
12月1號(hào)后——
Step2: 數(shù)據(jù)集降維
spleen_monocle <- setOrderingFilter(spleen_monocle, ordering_genes)
plot_ordering_genes(spleen_monocle)
spleen_monocle <- reduceDimension(spleen_monocle, max_components = 2, reduction_method = "DDRTree")
Step3: 將細(xì)胞按照偽時(shí)間排序
spleen_monocle <- orderCells(spleen_monocle)
colnames(pData(spleen_monocle))
[1] "nGene"
[2] "nUMI"
[3] "orig.ident"
[4] "Size_Factor"
[5] "Cluster"
[6] "peaks"
[7] "halo"
[8] "delta"
[9] "rho"
[10] "nearest_higher_density_neighbor"
plot_cell_trajectory(spleen_monocle, color_by = "Cluster")
plot_cell_trajectory(spleen_monocle, color_by = "State")
GM_state <- function(cds){
+ if (length(unique(pData(cds)$State)) > 1){
+ T0_counts <- table(pData(cds)$State, pData(cds)$Cluster)[,"0"]
+ return(as.numeric(names(T0_counts)[which
+ (T0_counts == max(T0_counts))]))
+ } else {
+ return (1)
+ }
+ }
plot_cell_trajectory(spleen_monocle, color_by = "Pseudotime")
plot_cell_trajectory(spleen_monocle, color_by = "State") +
+ facet_wrap(~State, nrow = 1)
作圖
> blast_genes <- row.names(subset(fData(spleen_monocle), gene_short_name %in% c("CCL5","TRAC","GNLY")))
> plot_genes_jitter(spleen_monocle[blast_genes,],
+ grouping = "State",
+ min_expr = 0.1)
> blast_genes <- row.names(subset(fData(spleen_monocle), gene_short_name %in% c("HSPA1A","MS4A1","S100A8")))
> plot_genes_jitter(spleen_monocle[blast_genes,],
+ grouping = "State",
+ min_expr = 0.1)
下面是沒有跑過的Ｗ慊Ｐ劬浴！
sp_expressed_genes <- row.names(subset(fData(spleen_monocle),
num_cells_expressed >= 10))
spleen_monocle_filtered <- spleen_monocle[HSMM_expressed_genes,]
my_genes <- row.names(subset(fData(spleen_monocle_filtered),
gene_short_name %in% c("CDK1", "MEF2C", "MYH3")))
cds_subset <- spleen_monocle_filtered[my_genes,]
plot_genes_in_pseudotime(cds_subset, color_by = "Cluster")
Reference
[1]Jenei, Z.M., Széplaki, G., Merkely, B. et al. Cell Stress and Chaperones (2013) 18: 447. https://doi.org/10.1007/s12192-012-0399-2
[2]Choi JI, Kim SD, Kim SH, Lim DJ, Ha SK. Semi-quantitative analyses of hippocampal heat shock protein-70 expression based on the duration of ischemia and the volume of cerebral infarction in mice. J Korean Neurosurg Soc. 2014;55(6):307-12.
