前言
在前两期的推文中：bigSCale：大规模scRNA-seq分析工具和 bigSCale代码实操（一）：对scRNA-seq聚类、拟时序分析，分别介绍了bigScale的分析框架和原理以及如何使用bigScale2对单细胞数据进行聚类分型、表型分型以及拟时序分析。这也是bigSCale2的核心工具，那么，在本期推文中，我们来结合代码实操介绍一下bigSCale2的后两个工具--bigSCale2 GRN（Gene Regulatory Networks）和bigSCale2 iCells。
bigSCale2 GRN 是一种基于单细胞数据推断基因调控网络并量化基因中心性的一个工具，本次推文将展示如何从健康胰腺和2型糖尿病(T2D)胰腺推断 GRN 并进行比较。bigSCale2 iCells可以将任何大小的大型数据集减少为质量明显更高的小型数据集，本次推文会介绍.mtx和.h5两种输入格式的大型数据实操。
代码流程
第一步：构建健康胰腺的GRN
data("pancreas")
results.ctl=compute.network(expr.data = expr.ctl,gene.names = gene.names)
#查看结果
results.ctl$graph
IGRAPH 8974736 UN-- 2268 9490 --
+ attr: name (v/c), bigSCale.Degree (v/n), bigSCale.PAGErank (v/n), bigSCale.Closeness (v/n), bigSCale.Betweenness (v/n)
+ edges from 8974736 (vertex names):
[1] TGFBR3--RGS16 TGFBR3--SIX2 TGFBR3--ENTPD3 TGFBR3--GLP1R TGFBR3--ASB9 PINK1 --RNF220 PINK1 --APH1A PINK1 --PCBP1
[9] PINK1 --TEX261 PINK1 --TADA3 PINK1 --NELFA PINK1 --HNRNPA0 PINK1 --PNRC1 PINK1 --C7orf26 PINK1 --BRI3 PINK1 --CIZ1
[17] PINK1 --ENDOG PINK1 --TMEM203 PINK1 --LDB1 PINK1 --SRSF9 PINK1 --ERH PINK1
+ ... omitted several edges
我们可以看出健康胰腺网络有2268个节点(基因)和9490条边。BigScale2计算了基因(节点)的中心性并将它们存储在输出列表中,中心性是衡量基因在调控网络中的重要性的一个指标,让我们来看看健康网络中最核心的基因。
DT::datatable(results.ctl$centrality)
image.png
此外，BigScale2提供了将网络结果导出为cytoscape可直接读取的形式。
toCytoscape(G = results.ctl$graph,file.name = 'Ph_net.json')
image.png
第二步：构建健康胰腺和T2D胰腺的GRN
构建模型并提取两个网络的结果：
model=compute.network.model(expr.data = cbind(expr.ctl,expr.t2d))
results.ctl=compute.network(expr.data = expr.ctl,gene.names = gene.names,model = model)
results.t2d=compute.network(expr.data = expr.t2d,gene.names = gene.names,model = model)
健康的胰腺有2301个节点和9513个边，而2型糖尿病有5387个节点和48583个边。
在这里，作者建议进一步操纵网络，以获得两者之间相似数量的边，从而使两个网络的比对结果更有意义。
output=homogenize.networks(list(result.ctl,result,t2d))
result.ctl=output[[1]]
result.t2d=output[[2]]
处理后这两个网络有非常相似的边数量，分别为23745和23456条边。
最后，我们可以比较两个基因调控网络中节点中心性之间的差异：
comparison=compare.centrality(list(results.ctl$centrality,results.t2d$centrality),c('Control','TD2'))
DT::datatable(comparison$Degree)# 查看结果
image.png
第三步：iCells 大数据降维
2百万（2M）个细胞降维为1万5千（15K）个细胞（.mtx格式）
数据下载：https://oncoscape.v3.sttrcancer.org/atlas.gs.washington.edu.mouse.rna/downloads
out=iCells(file.dir = "gene_counts.txt",target.cells = 15000)
27万（270K）个细胞降维为8千（8K）个细胞（.h5格式）
rhdf5::h5ls("ica_cord_blood_h5.h5")
group name otype dclass dim
0 / GRCh38 H5I_GROUP
1 /GRCh38 barcodes H5I_DATASET STRING 384000
2 /GRCh38 data H5I_DATASET INTEGER 260473471
3 /GRCh38 gene_names H5I_DATASET STRING 33694
4 /GRCh38 genes H5I_DATASET STRING 33694
5 /GRCh38 indices H5I_DATASET INTEGER 260473471
6 /GRCh38 indptr H5I_DATASET INTEGER 384001
7 /GRCh38 shape H5I_DATASET INTEGER 2
8 / PYTABLES_FORMAT_VERSION H5I_GROUP
9 / VERSION H5I_GROUP
10 / titldfdffde H5I_GROUP
11 / title H5I_GROUP
12 / version H5I_GROUP
#将数据转换为 .mtx 格式
out=bigscale.convert.h5(input.file = "ica_cord_blood_h5.h5",output.file= "ica_cord_blood.mtx.tar.gz",counts.field = "GRCh38",filter.cells = 400)
barcodes=rhdf5::h5read(file = "ica_cord_blood_h5.h5",name = "GRCh38/barcodes")
barcodes[1]
[1] "MantonCB1_HiSeq_1-AAACCTGAGGAGTTGC-1"
IDs=as.factor(apply(X = barcodes,MARGIN = 1,FUN = substr,start=7,stop=9))
IDs=IDs[out$filtered.cells]
out=iCells(file.dir = "ica_cord_blood.mtx.tar.gz",target.cells = 8000,sample.conditions = IDs)
image.png
小结
截止到目前，有关bigSCale2的三个工具都已经介绍完了，可以发现bigSCale2是一个较为全面的单细胞数据分析工具。值得一提的是bigScale2 Core没有使用降维的方法，因此具有较为灵敏的标记识别能力。此外，bigSCale2还提供了推断单细胞数据集的基因调控网络及网络之间的比较和可视化的方法。过大的单细胞数据常常困扰着科研人员，而bigSCale2提供了将任何大小的大型数据集压缩成更高质量的较小数据集的工具，极大的方便了许多单细胞下游分析。
好啦，本期分享到这里就结束了，我们下期再会~
网友评论