10X_final

Plan for 10X genomics data

数据处理

策略	进展	备注
平衡各时间点测序深度，删除ESC中的feeder细胞	完成	22012 genes across 67127 cells
对每个时间点采样500个细胞，在PC上尝试新算法		测试数据

数据可视化

OSK 重编程数据

方案1(终止)

7W细胞的diffusion map无法在超算正常跑完，考虑采用更快速的censor参数运行，censor参数考虑mRNA的丢失情况，censor参数需要的是read的读长，因此输入的数据须为基因表达值。

(Option1) 取D7/8的细胞找高变异基因；(Option2) 取所有细胞找高变异基因。
加入censor参数，计算diffusion map。
对前几个DCs，构建k-NN graph，利用Louvain clustering聚类。
FR展开实现可视化。

方案2 (最保守但成功率最大)

大量细胞的diffusion map难以计算，考虑先用少量细胞计算出流形，在此基础上预测其他细胞在该流形上的定位。

计算1W个细胞的diffusion components
计算其他细胞在已有DC上的定位（destiny -> dm_predict）

方案3（终止）

考虑用UMAP算法代替diffusion map。UMAP算法为最近发表的流形学习算法，其不仅能像tSNE一样保留数据的局部特征（看到稳定的细胞状态），同时能保留一定的全局特征（如时间点之间的差异），更重要的是计算速度非常快（比tSNE快两个数量级）。

对数据做一次或二次AP聚类（全部细胞HVG或D7/D8 HVG）。
计算每个基因cluster的特征表达（为适应10X数据的drop out问题，考虑尝试Zero-inflated factor analysis代替PCA，但计算起来会慢），删除干扰分析的基因cluster（与细胞周期相关的）。
利用UMAP把数据降到~10维。
构建10维数据的k-NN网络，利用FR展开。

方案4（终止）

与方案3类似，采用另一套流形学习方法PAHTE，这套方法定义一套全新的核函数来度量细胞的相似性，并通过随机游走消除数据噪音。计算速度尚未测试过。

对数据做一次或二次AP聚类（全部细胞HVG或D7/D8 HVG）。
计算每个基因cluster的特征表达（为适应10X数据的drop out问题，考虑尝试Zero-inflated factor analysis代替PCA，但计算起来会慢），删除干扰分析的基因cluster（与细胞周期相关的）。
利用PHATE把数据降到2维。

PS

为了保持10x数据与C1数据的一致性，这次分析暂不考虑URD算法（Science斑马鱼发育路径研究），D7/8的细胞里可能有多支NR细胞，包括少量神经细胞，角质化细胞样细胞等。

IDEA

1. 真正意义的Single-cell Orientation Tracing

做过众多测试，我们认为先对细胞聚类，找marker基因，用marker基因的特征来计算Diffusion，一方面可以降低噪音，另一方面可以拿到较为稳定的分支结构。

该策略有几个重要优势：

先聚类并找每类的marker基因，可以把分群能力强的基因富集出来。同时一个基因被允许分到不同的cluster中。基因富集基因属于机器里 特征工程 的特征选择。数据决定分析的上限，算法决定分析的下限。
利用PC1表示每个cluster marker基因可以稳定群基因的特征，减少冗余信息和避免不同cluster因基因数差异导致的不平衡，确保最终的分支能展示出来。

存在的问题：

细胞分类容易受到细胞周期的影响，尤其是OSK的数据，周期影响尤为显著。

分支呈现可能不充分。例如OSK的多能性分支中不能将Best cell分出来。我相信多能性分支还有继续细分的可能。

解决办法：

在高变异基因中去除与细胞周期高相关的基因/计算出每个细胞所处的周期，回归矫正
对于每类细胞，都做二次聚类。二次聚类所用的基因为该类的marker基因，如果有新的细胞类型产生，检验它们的maker基因，计算所有细分的marker基因在所有细胞中的PC1.

计划 2018-5-30

重要紧急事件

1. 完成7W细胞OSK重编程分支图

方案1-3并行进行

方案1

用D3-8的5000个细胞找HVG，聚类，找DEG，能够做出分支图。将DEG放到所有D3-8中，看是否能出来相似的分支，如果能，就与MEF, D0, D1, ESC细胞进行拼接，再展示结果。

方案2

用所有D3-8细胞找HVG，聚类，找DEG（完成）。把DEG挑选出来做SOT（1h）。确认每个SOT的gene group保留数据的主要特征（1h）。然后再做DM+FR展示看是否出现分群（3h）。

方案3

用所有D3-8细胞找HVG，聚类，找DEG（完成）。直接计算DEG的PC1 pattern。确认每群marker基因的生物学功能，保留数据的主要特征（1h）。然后再做DM+FR展示看是否出现分群（3h）。

方案4

目前的拼图存在几个问题：1. NR与MEF细胞混杂在一起；2. RP与ES的联系不明显；3. ESC存在异质性。说明全局的kNN构建存在问题。问题1，3说明时间点的信号不够强，2可以用更多的分析结果辅助证明RP以及其中best cell的存在。解决1，3的问题，可以用每个时间点特异高表达的基因来构建全局kNN。

方案5

SOT在鉴定D3-8的分支失效，主要问题在基因聚类，一方面容易错分，另一方面两次AP可能抹掉大量信息。考虑仅用AP结果作为删除周期特征的手段，删去细胞周期基因后，对基因表达或第一次AP的结果做PCA，取前几个PC做diffusion map和FR。

紧急事件

1. CiPS预处理

基本流程

细胞质控，标准化，缩尾处理，找HVG，聚类，找DEG（8h）

2. CiPS分支图

方案1

SOT对差异表达基因聚类，观察每个gene group代表的意义（与细胞周期和换液的关系），挑选有意义的gene group。采样40000个细胞做diffusion+FR，如果效果好则预测所有细胞的DC pattern。（1day）

方案2

计算DEG的PC1 pattern，去除不重要的pattern，采样40000细胞做diffusion+FR，如果效果好则预测所有细胞的DC pattern。（1day）

重要事件

基于OSK的网络图聚类，平均每个cluster细胞的表达值，用一个点表示一个cluster，简化分支模型，并可在此结构上计算基因的伪时间表达（6h）。
计算7F每个分支上基因的伪时间表达，找出特征基因（4h）。
绘制OSK 所有细胞的UMAP。
寻找Best cells。1. 从D3-8的图中鉴定出RP分支，单独取出来做Louvain聚类/SC3聚类.
Maker基因的上游调控转录因子及Motif分析。
计算10X NR和C1 NR的相似性。

问题

OSK数据存在的问题

多能性分支富集出Fos家族的Motif，目标基因与管家基因有关，须想办法去掉Fos的TF再分析。
Wnt信号通路可能会抑制角质细胞的生成。

CiPS数据潜在问题

换液会导致细胞在D22天前后出现断层，导致细胞不连续。最科学的办法是用计算的方法矫正换液带来的异质性。另一种方案是做两次AP，把换液影响的基因去除。
可能需要考虑细胞周期的影响。
9万多个细胞在计算Diffusion map时存在困难。考虑矫正/去除换液影响后，采样4W个细胞跑diffusion，再预测所有细胞的DCs。

进程

SOT: 5/30 - 14:20 - 计算D3-8 top100 marker基因的欧式距离，作为AP聚类的度量。利用相关系数的度量无法将神经细胞的特征独立出来。欧式距离计算时间很长。Rerun: 5/30 - 16:09

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/osk/10x/2_2_cls_D38/1_4_manifold/2_SOT

umap: 5/30 - 15:27 - 执行所有细胞的PCA，作为UMAP的输入。

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/osk/10x/3_1_manifold/umap/1_Extract_log_data.R

cips: 5/30 - 15:01 - cips数据细胞QC。Finished： 5/30 - 15:50

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/ciPS/All_cells/preprocess/1_Cell_QC.R

7F: 5/30 - 15:45 - 开始计算7F每个分支的基因表达趋势

 
xxxxxxxxxx
F:/7f/Explore/Tracing/SCENIC2/Psudo_exp.R # br1_sce

Trajactory1: Initial state -> Failed end
Trajactory2: Initial state -> Reprogramming potential -> Non-iPS state
Trajactory3: Initial state -> Reprogramming potential -> Pluripotent state

cips: 5/30 - 19:57 - cips数据标准化，聚类，Find marker genes. 5/30 - 20：30终止

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/ciPS/All_cells/preprocess

umap: 5/30 - 20:45 - 完成所有细胞的UMAP（包括D2）

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/osk/10x/3_1_manifold/umap

umap: 5/30 - 23:39 - 完成所有细胞的UMAP（去除D2）

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/osk/10x/3_1_manifold/umap

SOT: 5/31 - 2:02 - 完成Round1 AP cluster 的PCA（9 pcs）和diffusion map

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/osk/10x/2_2_cls_D38/1_4_manifold/2_SOT/DM_on_top100.R

cips: 5/31 - 2:09 - 重启cips数据标准化，聚类，Find marker genes.

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/ciPS/All_cells/preprocess/2_Trim_Values.R

7F: 5/31 - 2:12 - 检验7因子三条分支随伪时间发生差异表达的基因。

 
xxxxxxxxxx
F:/7f/Explore/Tracing/SCENIC2/Psudo_exp.R

7F: 5/31 - 2:42 - 完成平滑基因表达。

 
xxxxxxxxxx
F:/7f/Explore/Tracing/SCENIC2/Psudo_exp.R

SOT: 5/31 - 10:45 - 完成DM on 9pcs的diffusion map，分支非常明显。 - 下一步在其上plot marker genes

参数（top 100 marker genes in each cell cluster, ）

 
xxxxxxxxxx
F:/SingleCell/CellPress/Review/10X/ALL_CELL/D38/2_1_4_Cls_without_cc/SOT_TOP100

SOT: 5/31 - 10:51 - 用宽松的阈值鉴定基因group，去除细胞周期后，用round 1 AP cluster做PCA再做diffusion map.

参数（DEG - padj < 1e-4, ~1000 genes，Pearson correlation，gr2 - 细胞周期，top 7 pcs）

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/osk/10x/2_2_cls_D38/1_4_manifold/2_SOT
F:\SingleCell\CellPress\Review\10X\ALL_CELL\D38\2_1_4_Cls_without_cc\SOT_025

SOT: 5/31 - 15：35 - 直接用去除周期后的round1 AP cluster做diffusion map。

参数（DEG - padj < 1e-4, ~1000 genes，Pearson correlation，gr2 - 细胞周期）

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/osk/10x/2_2_cls_D38/1_4_manifold/2_SOT/DM_on_lfc025_pearson_on_ap1.R

SOT: 5/31 - 15：37 - 直接用去除周期后的round1 AP cluster做diffusion map。

参数（DEG - top 100 in each cluster, ~500 genes，Pearson correlation，gr2 - 细胞周期）

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/osk/10x/2_2_cls_D38/1_4_manifold/2_SOT/DM_on_top100_on_ap1.R

Fig1: 6/1 - 9:33 - 抽提D3-8基因表达数据，制作水螅图。包括：

2D/3D时间点图
Louvain聚类展示
找末端3/4个分支的差异表达基因

聚类：

clusters	注释	颜色
15	多能性	666666
9	多能性根部	547ABA
12	神经	D3405A
10	代谢	934C9C
2	NR1-Cd34	A1BBA8
4	NR1-根部	DEBDB3
6	起点	FBDF8A

参数：all_cell_dcs.RData; usedDim = 5 dcs, k = 200, grid = "grid", seed = 3.

Color:

 
xxxxxxxxxx
Day = gsub("(mef|d0|d1|d2|d3|d5|d6|d7|d8|esc).*", "\\1", rownames(umap2d))
Day = factor(Day, levels = c("mef", "d0", "d1", "d3", "d5", "d6", "d7", "d8", "esc"))
color_scale = c("#93574F", "#AAC8EE", "#BCBE3F", "#D17EB2", "#EB803A", "#D72D2E", "#2A9B49", "#227EAD", "#9469B0")
names(color_scale) = c("mef", "d0", "d1", "d3", "d5", "d6", "d7", "d8", "esc")

 
xxxxxxxxxx
F:\SingleCell\CellPress\Review\10X\ALL_CELL\Conclusion_v1\Figure1_OSK\Branches

CiPS: 6/2 - 6:50 - 用15个clusters的marker基因的PC1计算3600个细胞的diffusion map

 
x
/public/home/lhlin/project/single_cell/ciPS/All_cells/HVG_Based/TOP100/2_Diffusion.R

CiPS: 6/3 - 15:49 - 用3600个细胞的diffusion map结果预测剩余细胞的在流形空间中的位置

 
xxxxxxxxxx
/public/home/lhlin/project/single_cell/ciPS/All_cells/HVG_Based/TOP100/2_Diffusion.R