scVelo于2020年發(fā)表在Nature Biotechnology榔至,對最初的RNA velocity研究及其附帶的軟件velocyto進行了多項改進蛙婴。
本教程是介紹如何使用python包scVelo對單細胞RNA-seq數(shù)據(jù)(scRNA-seq)進行RNA速率分析。
參考資料:
scvelo github網(wǎng)站:https://github.com/theislab/scvelo
scvelo官方文檔:https://scvelo.readthedocs.io/index.html
Seurat to RNA-Velocity教程:https://github.com/basilkhuder/Seurat-to-RNA-Velocity#multiple-sample-integration
因為scvelo為python包枝哄,seurat數(shù)據(jù)在使用scvelo進行分析時,需要獲取seurat對象中的meta-data數(shù)據(jù)阻荒,上面的教程可以指導如何生成scvelo輸入文件挠锥。
scvelo實戰(zhàn)教程:https://smorabit.github.io/tutorials/8_velocyto/
怎么理解RNA速率?
La Manno等人2018年在Nature雜志上發(fā)表《RNA velocity of single cells》,提出了RNA速率分析方法财松。
RNA豐度是分析單個細胞狀態(tài)的有力指標瘪贱。單細胞RNA測序通過定量的方式以高準確度、高靈敏度和高通量來揭示RNA豐度辆毡。然而菜秦,這種方法僅捕獲細胞群某個時間點的靜態(tài)快照,這對分析細胞群在時間維度的現(xiàn)象研究(例如胚胎發(fā)生或組織再生)構成了挑戰(zhàn)舶掖。
作者提出了RNA速率球昨,即基因表達狀態(tài)的時間導數(shù),在常見的單細胞RNA測序方案中眨攘,可以通過新生(未剪接)和成熟(剪接)mRNA的相對豐度來估計基因剪接和降解的速率主慰,區(qū)分未剪接和剪接的mRNA,來直接估計細胞基因表達的動態(tài)分化鲫售。剛轉錄出的mRNA包含外顯子和內含子共螺,經(jīng)過splicing切除內含子后,得到用于編碼蛋白的spliced mRNA情竹。spliced mRNA的豐度由未成熟mRNA的splicing速度和降解速率共同決定藐不。通過計算未剪接轉錄本和剪接轉錄本之間的比率來推斷細胞命運的狀態(tài)(過渡與穩(wěn)定)和方向性(軌跡)。
簡而言之,RNA速率分析使我們能夠使用轉錄動力學的數(shù)學模型推斷在scRNA-seq實驗中未直接觀察到的轉錄動力學雏蛮。我們可以使用RNA速率來確定在給定的目標細胞群中是否誘導或抑制了目標基因涎嚼。此外,我們可以通過偽時間軌跡推斷這些信息來預測細胞命運決定挑秉。
問題1:通過splicing和spliced mRNA信息如何推斷分化方向法梯?
在時間軸上,mRNA是未剪接到剪接的過程犀概。scRNA-seq實驗只是捕獲細胞群某個時間點的靜態(tài)快照立哑,也就是細胞群一瞬間的轉錄組信息,時間點S1阱冶。假設分析A刁憋,B,C三種細胞類型木蹬,我們可以計算到在S1時間點至耻,每個細胞類型的mRNA剪接率,ABC剪接率由低到高镊叁,怎么能得到分化方向是A->B->C這個結論?只是一個時間點尘颓,而不是一個時間軸,這塊理解還是很斷線...
分析步驟:
本教材將演示如何將處理過/標準化的Seurat對象與RNA速率分析結合使用晦譬。Seurat是基于R語言的斩郎,但所有RNA Velocity軟件/包都是Python弹澎,因此我們將seurat單細胞數(shù)據(jù)遷移至python環(huán)境。
- 第一步:velocyto生成單個樣本的loom文件,每個樣本單獨生成匣沼;這步前提要cellranger分析夭问;
- 第二步:獲取seurat分析中的meta-data數(shù)據(jù)倡蝙,一般只對關注的某些細胞類型做RNA速率分析予颤,獲取umap坐標信息,cluster和celltype信息生棍;這步前提要seurat聚類分析颤霎;
- 第三步:運行python程序,調用scvelo涂滴,整合loom文件和meta-data數(shù)據(jù)友酱,運行RNA速率,繪圖柔纵;
如果有多個樣本缔杉,需要整合多個樣本的loom文件,生成combined.loom搁料;
安裝程序:
- scVelo(用于RNA Velocity)
- Velocyto 或 Kallisto Bustools(用于生成我們的初始RNA Velocity對象)
- Anndata(用于操作我們的RNA Velocity對象)
- Seurat
- Samtools -- 可選(Velocyto運行Samtools排序未分類.bam)
一或详、velocyto生成loom文件
首先进苍,我們需要將Seurat分析中使用的每個單細胞樣本生成loom文件(一種為基因組數(shù)據(jù)集設計的文件格式,例如單細胞)鸭叙。loom文件與Seurat 中使用的文件格式不同。必須根據(jù)樣品的原始FASTQ或BAM文件生成loom文件拣宏。生成Loom文件有兩種方法是:Velocyto命令行和Kallisto Bustools(KB.)
Kallisto Bustools沒有嘗試過沈贝,我們用Velocyto實操。
velocyto github網(wǎng)站:https://github.com/velocyto-team/velocyto.py
velocyto文檔:https://velocyto.org/velocyto.py
1.1 安裝velocyto包的依賴
conda install numpy scipy cython numba matplotlib scikit-learn h5py click
1.2 安裝velocyto
pip install velocyto
velocyto --help
1.3 示例腳本
velocyto run10x是針對10X的樣本測序數(shù)據(jù)的命令勋乾,運行前要進行cellranger分析宋下,cellranger分析完后生成bam文件。
velocyto將樣本的bam文件生成loom文件辑莫,此步驟耗時学歧,建議在HPC計算節(jié)點執(zhí)行。
velocyto.py的命令行參數(shù):參考
velocyto run10x -m repeat_msk.gtf mypath/sample01 somepath/refdata-cellranger-mm10-1.2.0/genes/genes.gtf
-s, --metadatatable FILE 樣本的注釋文件各吨,csv格式枝笨,行是樣本,列是注釋信息揭蜒。(跟seurat横浑,monocle等R包中的metadata格式是一樣的)
-m, --mask FILE 去掉一些重復表達的基因, .gtf文件屉更,可以從UCSE gene browser上下載
-l, --logic TEXT 用于篩選的邏輯值(這個參數(shù)還沒有研究過)
-M, --multimap 考慮非唯一映射(作者不建議這樣)
-@, --samtools-threads INTEGER 用samtools對bam文件進行sort時用到的進程數(shù)
--samtools-memory INTEGER 每個進程分配的內存(單位Mb)
-t, --dtype TEXT loom文件層的dtype徙融。如果每個細胞的每個基因超過6000個molecules/reads,這個參數(shù)建議設置為32瑰谜。(默認是16)
mypath/sample01:
sample01的cellranger結果目錄欺冀,該目錄包含outs、outs/analys和outs/filtered_feature_bc_matrix子文件夾萨脑,包含bam文件隐轩。
genes.gtf:
cellranger分析流程的參考基因組的注釋文件
path/cellranger/reference/refdata-gex-mm10-2020-A/genes/genes.gtf
path/cellranger/reference/refdata-gex-GRCh38-2020-A/genes/genes.gtf
repeat_msk.gtf:
重復表達的基因.gtf文件,可以從UCSE gene browser上下載
human_allTracks.gtf:小鼠的鏈接
mm10_allTracks.gtf:小鼠的鏈接
百度云盤下載地址:
path/to/RNA_velocity/human_allTracks.gtf
path/to/RNA_velocity/mm10_allTracks.gtf
鏈接:https://pan.baidu.com/s/1-dk1CQUa-39UHIxJ0wWQpg 提取碼:1234
1.4 運行示例項目
# mouse
velocyto run10x -m path/to/RNA_velocity/mm10_allTracks.gtf /path/to/scrna_project/Cellranger_result/sample1/ path/to/cellranger/reference/refdata-gex-mm10-2020-A/genes/genes.gtf -@ 10 --samtools-memory 24000
# human
velocyto run10x -m path/to/RNA_velocity/human_allTracks.gtf /path/to/scrna_project/Cellranger_result/sample1/ path/to/cellranger/reference/refdata-gex-GRCh38-2020-A/genes/genes.gtf -@ 10 --samtools-memory 24000
運行上面的程序老是出現(xiàn)報錯:
no such file or directory: 'samtools'
出現(xiàn)的問題說明:沒有安裝samtools和libcrypto.so.1.0.0
安裝samtools
conda install samtools=1.5
samtools出現(xiàn)的報錯
samtools: error while loading shared libraries: libcrypto.so.1.0.0
嘗試過的方法:
conda install -c bioconda samtools openssl=1.0
conda install samtools openssl=1.0
最終的解決方法:
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/main
conda config --add channels https://mirrors.bfsu.edu.cn/anaconda/pkgs/free
conda install samtools openssl
samtools
重新運行砚哗,成功龙助,寫成python腳本去集群投遞任務。
二蛛芥、從seurat對象中提取meta-data
(R語言運行)已經(jīng)在R中使用Seurat執(zhí)行了主要的數(shù)據(jù)處理(過濾提鸟、歸一化、聚類仅淑、批次處理称勋、降維,細胞類型注釋)涯竟,接下來赡鲜,提取meta-data的主要信息:
- Filtered Cell Ids
- UMAP or TSNE coordinates
- Clusters/Celltype (Optional)
- Clusters/Celltype Colors (Optional)
makedirs <- function(dir_list){
for (i_dir in dir_list){
if (!dir.exists(i_dir) ){
dir.create(i_dir, recursive = TRUE)
}
}
}
#--------------------------------------------------
makedirs(file.path(project_work_dir, "rna_velocity"))
setwd(file.path(project_work_dir, "rna_velocity"))
table(seurat_obj$celltype)
# 獲得每個細胞的UMAP或TSNE坐標空厌,使用 Embeddings函數(shù)
write.csv(Embeddings(seurat_obj, reduction = "umap"), file = "cell_embeddings.csv")
# 獲取每個細胞的barcode
write.csv(Cells(seurat_obj), file = "cellID_obs.csv", row.names = FALSE)
# 提取每個細胞的cluster信息
write.csv(seurat_obj@meta.data[, 'cluster', drop = FALSE], file = "cell_clusters.csv")
# 提取每個細胞的celltype信息
write.csv(seurat_obj@meta.data[, 'celltype', drop = FALSE], file = "cell_celltype.csv")
# 獲取celltype的顏色信息
hue_pal()(length(levels(seurat_obj$celltype)))
# 獲取cluster的顏色信息
hue_pal()(length(levels(seurat_obj$cluster)))
# 繪制umap圖,與RNA速率圖對比看
Idents(seurat_obj) <- "celltype"
alpha.use <- 0.8
p1 <- DimPlot(object = seurat_obj, reduction = "umap", label = TRUE, label.size = 5, pt.size=0.4, raster=FALSE)
p1$layers[[1]]$mapping$alpha <- alpha.use
save_plot("UMAP_clustering.png", p1, base_height = 8, base_aspect_ratio = 1.3, base_width = NULL, dpi=600)
save_plot("UMAP_clustering.pdf", p1, base_height = 8, base_aspect_ratio = 1.3, base_width = NULL)
seurat對象中的barcode是有樣本名稱前綴的银酬,sample1_, sample2_...
三嘲更、整合Loom文件和meta-data
(python運行)整合loom文件和meta-data數(shù)據(jù)
代碼中有對loom文件barcode的處理,與seurat的barcode一致揩瞪。
# 整合多個loom文件
import loompy
files=['path/to/sample_one','path/to/sample_two']
output_filename='/save/path/to/combined.loom'
loompy.combine(files, output_filename, key="Accession")
# 整合loom文件和meta-data數(shù)據(jù)
import scvelo as scv
import pandas as pd
import numpy as np
import os
loom_data = scv.read('/save/path/to/combined.loom', cache=False)
loom_data.obs
# barcode名字去重后綴x赋朦,與seurat導出的barcode名稱一致
loom_data.obs = loom_data.obs.rename(index = lambda x: x.replace(':', '-').replace('x', ''))
loom_data.obs.head()
# 讀取seurat中的meta信息
meta_path = "path/to/rna_velocity"
sample_obs = pd.read_csv(os.path.join(meta_path, "cellID_obs.csv"))
cell_umap= pd.read_csv(os.path.join(meta_path, "cell_embeddings.csv"), header=0, names=["Cell ID", "UMAP_1", "UMAP_2"])
cell_clusters = pd.read_csv(os.path.join(meta_path, "cell_clusters.csv"), header=0, names=["Cell ID", "cluster"])
cell_celltype = pd.read_csv(os.path.join(meta_path, "cell_celltype.csv"), header=0, names=["Cell ID", "celltype"])
# 對細胞文件和RNA剪切速率文件取交集,保留關注的細胞類型
sample_one = loom_data[np.isin(loom_data.obs.index, sample_obs)]
sample_one.obs.head()
# 構建umap坐標, cluster, celltype信息數(shù)據(jù)框
sample_one_index = pd.DataFrame(sample_one.obs.index)
sample_one_index = sample_one_index.rename(columns = {0:'Cell ID'})
umap_ordered = sample_one_index.merge(cell_umap, on = "Cell ID")
umap_ordered.head()
celltype_ordered = sample_one_index.merge(cell_celltype, on = "Cell ID")
celltype_ordered.head()
clusters_ordered = sample_one_index.merge(cell_clusters, on = "Cell ID")
clusters_ordered.head()
# 將umap坐標與cluster信息加入sample_one
umap_ordered = umap_ordered.iloc[:,1:]
clusters_ordered = clusters_ordered.iloc[:,1:]
celltype_ordered = celltype_ordered.iloc[:,1:]
sample_one.obsm['X_umap'] = umap_ordered.values
sample_one.uns['clusters'] = clusters_ordered.values
sample_one.obs['celltype'] = celltype_ordered.values
adata = sample_one
# some gene labels are duplicated (Ensembl IDs are still unique!!)
adata.var_names_make_unique()
# save model to file
adata.write('Allcelltype_dynamicModel.h5ad', compression = 'gzip')
# read h5ad file
adata= scv.read('Allcelltype_dynamicModel.h5ad')
四李破、運行RNA速率
(python運行)可視化繪圖
# Running RNA Velocity
scv.pp.filter_and_normalize(adata,min_shared_counts=30, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.tl.velocity(adata, mode = "stochastic")
scv.tl.velocity_graph(adata)
scv.pl.velocity_embedding(adata, basis='X_umap', arrow_size=5)
ident_colours = ["#F8766D", "#A3A500", "#00BF7D", "#00B0F6", "#E76BF3"]
scv.pl.velocity_embedding_stream(adata, basis='X_umap',color = "celltype", palette = ident_colours)
scv.pl.velocity_embedding_stream(adata, basis='X_umap',color = "celltype", palette = ident_colours, size = 20,alpha =0.8)
scv.pl.velocity_embedding_stream(adata, basis='X_umap',color = "celltype", palette = ident_colours, size = 20,alpha =0.8, save="UMAP_stream.png", figsize=(7,5), legend_fontsize = 9, show=False, title='')
scv.pl.velocity_embedding_stream(adata, basis='X_umap',color = "celltype", palette = ident_colours, size = 20,alpha =0.8, save="UMAP_stream.pdf", figsize=(7,5), legend_fontsize = 9, show=False, title='')
以上是我常用的分析流程宠哄,下面學習下smorabit的教材,看哪些內容可以借鑒嗤攻。
Morabito教程: RNA velocity analysis with scVelo
Step-1:數(shù)據(jù)從Seurat 轉換為Python/anndata
將Seurat對象轉換為可在scVelo中可使用的數(shù)據(jù)文件毛嫉,共導出4個文件,分別是:1)seurat對象的metadata文件妇菱;2)細胞的基因表達矩陣counts.mtx承粤,文件很大;3)pca降維的pca坐標文件pca.csv恶耽,有多列密任,選擇了多少PCA,就有多少列偷俭。4)基因名稱列表gene_names.csv浪讳;
makedirs(file.path(project_work_dir, "rna_velocity"))
setwd(file.path(project_work_dir, "rna_velocity"))
table(seurat_obj$celltype)
# save metadata table:
seurat_obj$barcode <- colnames(seurat_obj)
seurat_obj$UMAP_1 <- seurat_obj@reductions$umap@cell.embeddings[,1]
seurat_obj$UMAP_2 <- seurat_obj@reductions$umap@cell.embeddings[,2]
write.csv(seurat_obj@meta.data, file='metadata.csv', quote=F, row.names=F)
# write expression counts matrix
library(Matrix)
counts_matrix <- GetAssayData(seurat_obj, assay='RNA', slot='counts')
writeMM(counts_matrix, file='counts.mtx')
# write dimesnionality reduction matrix, in this example case pca matrix
write.csv(seurat_obj@reductions$pca@cell.embeddings, file='pca.csv', quote=F, row.names=F)
# write gene names
write.table(
data.frame('gene'=rownames(counts_matrix)),file='gene_names.csv',
quote=F,row.names=F,col.names=F
)
(python運行)seurat數(shù)據(jù)對象轉adata數(shù)據(jù)結構數(shù)據(jù)
import scanpy as sc
import anndata
from scipy import io
from scipy.sparse import coo_matrix, csr_matrix
import numpy as np
import os
import pandas as pd
# load sparse matrix:
X = io.mmread("counts.mtx")
# create anndata object
adata = anndata.AnnData(
X=X.transpose().tocsr()
)
# load cell metadata:
cell_meta = pd.read_csv("metadata.csv")
# load gene names:
with open("gene_names.csv", 'r') as f:
gene_names = f.read().splitlines()
# set anndata observations and index obs by barcodes, var by gene names
adata.obs = cell_meta
adata.obs.index = adata.obs['barcode']
adata.var.index = gene_names
# load dimensional reduction:
pca = pd.read_csv("pca.csv")
pca.index = adata.obs.index
# set pca and umap
adata.obsm['X_pca'] = pca.to_numpy()
adata.obsm['X_umap'] = np.vstack((adata.obs['UMAP_1'].to_numpy(), adata.obs['UMAP_2'].to_numpy())).T
# plot a UMAP colored by sampleID to test:
sc.pl.umap(adata, color=['celltype'], frameon=False, save=True)
# save dataset as anndata format
adata.write('my_data.h5ad')
# reload dataset
adata = sc.read_h5ad('my_data.h5ad')
Step 2: 構建spliced和unspliced計數(shù)矩陣
與我們在Seurat中使用的相同的基于UMI的基因計數(shù)矩陣不同,我們需要一個splicing和spliced轉錄本的矩陣涌萤。 我們可以使用velocy命令行工具或使用Kallisto-Bustools構建這個矩陣淹遵。在這里我使用了 velocyto command line tool,,僅僅是因為我個人從來沒有嘗試過Kallisto-Bustools负溪。
velocy命令行工具直接從cellranger輸出目錄開始運行透揣,需要提供.bam文件,它也可以在任何單細胞平臺上使用川抡。我們還必須為你的物種提供參考.gtf注釋(此處使用mm10)辐真,并且你可以選擇提供.gtf以掩蓋重復區(qū)域(velocy推薦這種方式)。
repeats="/path/to/repeats/mm10_rmsk.gtf"
transcriptome="/path/to/annoation/file/gencode.vM25.annotation.gtf"
cellranger_output="/path/to/cellranger/output/"
velocyto run10x -m $repeats \
$cellranger_output \
$transcriptome
Step 3: 加載數(shù)據(jù)
前面我們已經(jīng)生成了scvelo的輸入數(shù)據(jù)崖堤,可以將其加載到python中侍咱。Velocy為每個樣本創(chuàng)建了一個單獨的拼接和未拼接矩陣,因此我們首先必須將不同的樣本合并為一個對象密幔。此外楔脯,需要處理loom文件的細胞barcode名稱,以使anndata對象與基因表達矩陣數(shù)據(jù)相匹配胯甩。
import scvelo as scv
import scanpy as sc
import cellrank as cr
import numpy as np
import pandas as pd
import anndata as ad
scv.settings.verbosity = 3
scv.settings.set_figure_params('scvelo', facecolor='white', dpi=100, frameon=False)
cr.settings.verbosity = 2
adata = sc.read_h5ad('my_data.h5ad')
# load loom files for spliced/unspliced matrices for each sample:
ldata1 = scv.read('Sample1.loom', cache=True)
ldata2 = scv.read('Sample2.loom', cache=True)
ldata3 = scv.read('Sample3.loom', cache=True)
Variable names are not unique. To make them unique, call `.var_names_make_unique`.
Variable names are not unique. To make them unique, call `.var_names_make_unique`.
Variable names are not unique. To make them unique, call `.var_names_make_unique`.
# rename barcodes in order to merge:
barcodes = [bc.split(':')[1] for bc in ldata1.obs.index.tolist()]
# 查看seurat的barcode命名方式昧廷,也可以是'Sample1_'
barcodes = [bc[0:len(bc)-1] + '_10' for bc in barcodes]
ldata1.obs.index = barcodes
barcodes = [bc.split(':')[1] for bc in ldata2.obs.index.tolist()]
barcodes = [bc[0:len(bc)-1] + '_11' for bc in barcodes]
ldata2.obs.index = barcodes
barcodes = [bc.split(':')[1] for bc in ldata3.obs.index.tolist()]
barcodes = [bc[0:len(bc)-1] + '_9' for bc in barcodes]
ldata3.obs.index = barcodes
# make variable names unique
ldata1.var_names_make_unique()
ldata2.var_names_make_unique()
ldata3.var_names_make_unique()
# concatenate the three loom
ldata = ldata1.concatenate([ldata2, ldata3])
# merge matrices into the original adata object
adata = scv.utils.merge(adata, ldata)
# plot umap to check
sc.pl.umap(adata, color='celltype', frameon=False, legend_loc='on data', title='', save='_celltypes.pdf')
Step 4: 使用scVelo計算RNA速率
我們使用scVelo來計算RNA速率堪嫂。scVelo允許用戶使用2018年原始出版物中的穩(wěn)態(tài)模型(steady-state model)以及 2020年出版物中更新的動態(tài)模型( dynamical model)。
首先木柬,我們檢查每個細胞類型 中spliced和unspliced的占比皆串。接下來我們執(zhí)行一些預處理,然后我們使用穩(wěn)態(tài)模型(隨機選項)計算RNA速率眉枕。
scv.pl.proportions(adata, groupby='celltype_full')
# pre-process
scv.pp.filter_and_normalize(adata)
scv.pp.moments(adata)
WARNING: Did not normalize X as it looks processed already. To enforce normalization, set `enforce=True`.
Normalized count data: spliced, unspliced.
WARNING: Did not modify X as it looks preprocessed already.
computing neighbors
finished (0:00:14) --> added
'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)
computing moments based on connectivities
finished (0:00:17) --> added
'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)
# compute velocity
scv.tl.velocity(adata, mode='stochastic')
scv.tl.velocity_graph(adata)
computing velocities
finished (0:01:15) --> added
'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph
finished (0:15:00) --> added
'velocity_graph', sparse matrix with cosine correlations (adata.uns)
Step 5: 可視化RNA速率場
我們可以通過可視化二維降維頂部的矢量場愚战,對所有基因和所有細胞的RNA速率進行廣泛的可視化。
scv.pl.velocity_embedding(adata, basis='umap', frameon=False, save='embedding.pdf')
computing velocity embedding
finished (0:00:06) --> added
'velocity_umap', embedded velocity vectors (adata.obsm)
saving figure to file ./figures/scvelo_embedding.pdf
scv.pl.velocity_embedding_grid(adata, basis='umap', color='celltype', save='embedding_grid.pdf', title='', scale=0.25)
scv.pl.velocity_embedding_stream(adata, basis='umap', color=['celltype', 'condition'], save='embedding_stream.pdf', title='')
我們還可以可視化我們感興趣的基因的動態(tài)變化齐遵。在這里,我們展示了Ptgds基因的未剪接轉錄本與剪接轉錄本的比率塔插,以及在feature plots圖中顯示基因的速率和基因表達值梗摇。
# plot velocity of a selected gene
scv.pl.velocity(adata, var_names=['Ptgds'], color='celltype')
Morabito教程是把seurat轉成python環(huán)境的Anndata對象,再和loom數(shù)據(jù)進行merge想许;
而我們常見的做法是將seurat的meta-data和loom數(shù)據(jù)進行merge伶授;
后面的下游分析和對RNA速率算法理解留到下一篇文章。未完待續(xù)...