使用PHATE復(fù)現(xiàn)Science Immunology文章的結(jié)果

在上篇文章中仪搔，我們初步探索了PHATE的使用方法憎账，發(fā)現(xiàn)它在揭示一些連續(xù)分化過程中不同細(xì)胞狀態(tài)之間的微小局部差異具有很好的效果领铐，同時也能保留細(xì)胞全局的整體結(jié)構(gòu)毛雇。在本節(jié)教程中泼菌，我們將復(fù)現(xiàn)演示近期發(fā)表在Science Immunology期刊上的一篇文章的結(jié)果括丁，進(jìn)一步學(xué)習(xí)PHATE的相關(guān)使用方法陨亡。

image.png

原文鏈接：A reservoir of stem-like CD8+ T cells in the tumor-draining lymph node preserves the ongoing anti-tumor immune response. 2021, Science Immunology

文章結(jié)果圖：


image.png

image.png

本文復(fù)現(xiàn)圖：


image.png
image.png

文章數(shù)據(jù)下載

文章處理后的基因表達(dá)矩陣文件存放在GEO數(shù)據(jù)庫中犀填，檢索號為 GSE182509舵稠，這里只提供了pkl格式的表達(dá)矩陣超升，需要用python的pickle包進(jìn)行讀取，我已將其轉(zhuǎn)換為TSV文件存放在我的百度云盤中哺徊，有需要的可以下載使用室琢。

鏈接：https://pan.baidu.com/s/1IoSIYoEfTzZarLWXWvgvzg
提取碼：gkd9

image.png

加載示例數(shù)據(jù)

# 導(dǎo)入所需的python包
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import phate
import scprep
import magic

# 加載示例數(shù)據(jù)
chronic = scprep.io.load_tsv("GSM5530565_Chronic_LCMV_processed.matrix.txt")
acute = scprep.io.load_tsv("GSM5530564_Acute_LCMV_processed.matrix.txt")
tumor8 = scprep.io.load_tsv("GSM5530560_Lung_week8_processed.matrix.txt")
dln8 = scprep.io.load_tsv("GSM5530561_LN_week8_processed.matrix.txt")
tumor17 = scprep.io.load_tsv("GSM5530562_Lung_week17_processed.matrix.txt")
dln17 = scprep.io.load_tsv("GSM5530563_LN_week17_processed.matrix.txt")

# 查看示例數(shù)據(jù)
chronic.head()
#                    Mrpl15 (ENSMUSG00000033845)  ...  CAAA01147332.1 (ENSMUSG00000095742)
#0                                                ...
#AAACCTGAGAGTAAGG-1                          0.0  ...                             0.000000
#AAACCTGGTGTGAAAT-1                          0.0  ...                             0.000000
#AAACGGGAGGCTCATT-1                          0.0  ...                             0.000000
#AAACGGGGTAGCTTGT-1                          0.0  ...                             0.000000
#AAACGGGGTCGAATCT-1                          0.0  ...                             0.552073
#[5 rows x 11595 columns]
chronic.shape
#(1185, 11595)
acute.shape
#(10768, 12960)
tumor8.shape
#(806, 11749)
dln8.shape
#(1742, 12116)
tumor17.shape
#(731, 11150)
dln17.shape
#(876, 11595)

數(shù)據(jù)質(zhì)控預(yù)處理和數(shù)據(jù)歸一化

由于下載的數(shù)據(jù)是預(yù)先已經(jīng)進(jìn)行質(zhì)控過濾和歸一化處理的，這里將不再進(jìn)行處理落追。詳細(xì)用法見上期 使用PHATE進(jìn)行單細(xì)胞高維數(shù)據(jù)的可視化

使用PHATE進(jìn)行低維嵌入降維可視化

### analysis for chronic sample ###
#Embedding Data Using PHATE
# Instantiating the PHATE estimator
phate_operator = phate.PHATE(n_jobs=-2)

Y_phate = phate_operator.fit_transform(chronic)
Y_phate
#array([[ -5.34563999,   6.11670333],
#       [ 29.8771147 ,   8.33029219],
#       [  9.44218856, -25.94568946],
#       ...,
#       [-10.86402299,  -1.15554   ],
#       [-14.65615727,  -1.03794057],
#       [-16.53150258,   4.50993521]])

細(xì)胞聚類分群

# cell clustering
clusters = phate.cluster.kmeans(phate_operator, n_clusters=8)
clusters
#array([0, 2, 1, ..., 4, 4, 4], dtype=int32)

# 聚類結(jié)果可視化
scprep.plot.scatter2d(Y_phate, c=clusters, figsize=(8,7), cmap="Spectral",
                      ticks=False, label_prefix="PHATE", title= "Chronic LCMV")
plt.savefig("plot_phate_chronic_2d_by_cluster.png")
image.png

使用MAGIC進(jìn)行數(shù)據(jù)填充

#Creating the MAGIC operator
magic_op = magic.MAGIC()

# Running MAGIC for all genes
chronic_magic = magic_op.fit_transform(chronic, genes='all_genes')
chronic_magic.head()
#                    Mrpl15 (ENSMUSG00000033845)  ...  CAAA01147332.1 (ENSMUSG00000095742)
#0                                                ...
#AAACCTGAGAGTAAGG-1                     0.150924  ...                             0.076723
#AAACCTGGTGTGAAAT-1                     0.153466  ...                             0.067584
#AAACGGGAGGCTCATT-1                     0.146601  ...                             0.064838

# rename columns names
chronic_magic.columns = [i.split(" ")[0] for i in chronic_magic.columns.tolist()]

# marker基因可視化
markers = ["Sell", "Ccr7", "Tcf7", "Slamf6", "Xcl1", "Il7r", "Eomes", "Tbx21", "Gzmb", "Prf1", "Pdcd1", "Havcr2", "Cd101", "Cx3cr1", "Cxcr6"]

for marker in markers:
        # 2d plot
        scprep.plot.scatter2d(Y_phate, c=chronic_magic[marker], figsize=(8,7), cmap="Reds",
                        ticks=False, label_prefix="PHATE", title=marker + " magic expression")
        plt.savefig("plot_chronic_magic_marker_2d_" + marker + ".png")
image.png

根據(jù)這些marker基因的表達(dá)情況盈滴，我們將不同的cluster進(jìn)行細(xì)胞類型的注釋，得到以下的細(xì)胞注釋結(jié)果轿钠。

Tnaive：Sell, Ccr7
Tsl: Tcf7, Slamf6, Xcl1
Ttrans: Cx3cr1, Cxcr6
Tex: Pdcd1, Havcr2, Cd101

image.png

多樣本合并分析

### analysis for combined five sample ###
# Merge all datasets and create a vector representing the time point of each sample
alldata = [chronic,tumor8,tumor17,dln8,dln17]

EBT_counts, sample_labels = scprep.utils.combine_batches(
    alldata,
    ["Chronic","Early Tumor","Late Tumor","Early LN","Late LN"],
    append_to_cell_names=True
)
del alldata # removes objects from memory

EBT_counts.head()
#                            Mrpl15 (ENSMUSG00000033845)  ...  CAAA01147332.1 (ENSMUSG00000095742)
#AAACCTGAGAGTAAGG-1_Chronic                          0.0  ...                             0.000000
#AAACCTGGTGTGAAAT-1_Chronic                          0.0  ...                             0.000000
#AAACGGGAGGCTCATT-1_Chronic                          0.0  ...                             0.000000
#AAACGGGGTAGCTTGT-1_Chronic                          0.0  ...                             0.000000
#AAACGGGGTCGAATCT-1_Chronic                          0.0  ...                             0.552073
#[5 rows x 10246 columns]

EBT_counts.shape
#(5340, 10246)

sample_labels
#AAACCTGAGAGTAAGG-1_Chronic    Chronic
#AAACCTGGTGTGAAAT-1_Chronic    Chronic
#AAACGGGAGGCTCATT-1_Chronic    Chronic
#AAACGGGGTAGCTTGT-1_Chronic    Chronic
#AAACGGGGTCGAATCT-1_Chronic    Chronic
#                               ...
#Name: sample_labels, Length: 5340, dtype: object

PHATE降維可視化

#Embedding Data Using PHATE
# Instantiating the PHATE estimator
phate_operator = phate.PHATE(n_jobs=-2)

Y_phate = phate_operator.fit_transform(EBT_counts)
Y_phate
#array([[ 79.05651647,  15.42592929],
#       [ 20.72815444,  25.65566379],
#       [ 90.57712893,  -4.77917562],
#       ...,
#       [-52.39011592,  39.20142516],
#       [-28.51731009,  -8.66499775],
#       [-46.00734805, -17.37265621]])

scprep.plot.scatter2d(Y_phate, c=sample_labels, figsize=(10,8), cmap="Spectral",
                      ticks=False, label_prefix="PHATE")
plt.savefig("plot_phate_2d_by_sample.png")
test4.png
# 3D visualization
phate_operator.set_params(n_components=3)
Y_phate_3d = phate_operator.transform()
Y_phate_3d
#array([[ 77.85894712,  13.21732854, -11.29708591],
#       [ 17.23971363,  19.47621167, -27.54015408],
#       [ 87.69788489,   4.0231593 ,  13.16082805],
#       ...,
#       [-46.46850791,  42.68065699,  -2.93836416],
#       [-20.53357918,  -4.9821892 , -22.70265013],
#       [-36.18351133,  -9.56494547, -29.67206442]])

scprep.plot.scatter3d(Y_phate_3d, c=sample_labels, figsize=(8,6), cmap="Spectral",
                      ticks=False, label_prefix="PHATE")
plt.savefig("plot_phate_3d_by_sample.png")
test5.png

細(xì)胞聚類分群

# cell clustering
clusters = phate.cluster.kmeans(phate_operator, n_clusters=12)
clusters
#array([8, 6, 1, ..., 2, 4, 4], dtype=int32)

# save meta data
meta = pd.merge(pd.DataFrame(sample_labels),pd.DataFrame(clusters,index=sample_labels.index,columns=["cluster"]),left_index=True,right_index=True)
meta.head()
#                           sample_labels  cluster
#AAACCTGAGAGTAAGG-1_Chronic       Chronic        8
#AAACCTGGTGTGAAAT-1_Chronic       Chronic        6
#AAACGGGAGGCTCATT-1_Chronic       Chronic        1
#AAACGGGGTAGCTTGT-1_Chronic       Chronic        8
#AAACGGGGTCGAATCT-1_Chronic       Chronic        8

meta.to_csv("metadata.csv")

scprep.plot.scatter3d(Y_phate_3d, c=clusters, figsize=(8,6), cmap="Spectral",
                      ticks=False, label_prefix="PHATE")
plt.savefig("plot_phate_3d_by_cluster.png")
test6.png
scprep.plot.scatter2d(Y_phate, c=clusters, figsize=(8,6), cmap="Spectral",
                      ticks=False, label_prefix="PHATE")
plt.savefig("plot_phate_2d_by_cluster.png")
test7.png
# to save as a gif:
scprep.plot.rotate_scatter3d(Y_phate_3d, c=sample_labels,
                             figsize=(8,6), cmap="Spectral",
                             ticks=False, label_prefix="PHATE", filename="phate.gif")
phate.gif

使用MAGIC進(jìn)行數(shù)據(jù)填充

# rename columns names
EBT_counts.columns = [i.split(" ")[0] for i in EBT_counts.columns.tolist()]

# MAGIC imputation
#Creating the MAGIC operator
magic_op = magic.MAGIC()

#Running MAGIC with gene selection
#bmmsc_magic = magic_op.fit_transform(bmmsc_data, genes=["Mpo", "Klf1", "Ifitm1"])
#bmmsc_magic.head()

#Visualizing gene-gene relationships
#fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16, 6))

#scprep.plot.scatter(x=bmmsc_data['Mpo'], y=bmmsc_data['Klf1'], c=bmmsc_data['Ifitm1'],  ax=ax1,
#                    xlabel='Mpo', ylabel='Klf1', legend_title="Ifitm1", title='Before MAGIC')

#scprep.plot.scatter(x=bmmsc_magic['Mpo'], y=bmmsc_magic['Klf1'], c=bmmsc_magic['Ifitm1'], ax=ax2,
#                    xlabel='Mpo', ylabel='Klf1', legend_title="Ifitm1", title='After MAGIC')
#plt.tight_layout()
#plt.show()

# Running MAGIC for all genes
EBT_counts_magic = magic_op.fit_transform(EBT_counts, genes='all_genes')
EBT_counts_magic.head()

markers = ["Sell", "Ccr7", "Tcf7", "Slamf6", "Xcl1", "Il7r", "Eomes", "Tbx21", "Gzmb", "Prf1", "Pdcd1", "Havcr2", "Cd101", "Cx3cr1", "Cxcr6"]

for marker in markers:
        # 2d plot
        scprep.plot.scatter2d(Y_phate, c=EBT_counts_magic[marker], figsize=(8,6), cmap="Reds",
                        ticks=False, label_prefix="PHATE", title=marker + " magic expression")
        plt.savefig("plot_magic_marker_2d_" + marker + ".png")
        # 3d plot
        scprep.plot.scatter3d(Y_phate_3d, c=EBT_counts_magic[marker], figsize=(8,6), cmap="Reds",
                        ticks=False, label_prefix="PHATE", title=marker + " magic expression")
        plt.savefig("plot_magic_marker_3d_" + marker + ".png")
image.png

基因差異表達(dá)分析

# Differential analysis
# By samples
de_samples = scprep.stats.differential_expression_by_cluster(EBT_counts_magic,clusters=sample_labels,direction="up")
de_samples
de_samples["Chronic"]

for key,value in de_samples.items():
        print(value.head(n=10))

# save DE data
de_samples["Chronic"].to_csv("DEs_Chronic.csv")
de_samples["Early Tumor"].to_csv("DEs_Early_Tumor.csv")
de_samples["Early LN"].to_csv("DEs_Early_LN.csv")

# By clusters
de_clusters = scprep.stats.differential_expression_by_cluster(EBT_counts_magic,clusters=clusters,direction="up")
de_clusters
de_clusters[0]

for key,value in de_clusters.items():
        print(value.head(n=10))
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末巢钓，一起剝皮案震驚了整個濱河市，隨后出現(xiàn)的幾起案子疗垛，更是在濱河造成了極大的恐慌症汹，老刑警劉巖，帶你破解...
    沈念sama閱讀 206,126評論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件继谚，死亡現(xiàn)場離奇詭異烈菌，居然都是意外死亡阵幸，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,254評論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門芽世，熙熙樓的掌柜王于貴愁眉苦臉地迎上來挚赊，“玉大人，你說我怎么就攤上這事济瓢≤睿” “怎么了？”我有些...
    開封第一講書人閱讀 152,445評論 0贊 341
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵旺矾，是天一觀的道長蔑鹦。 經(jīng)常有香客問我，道長箕宙，這世上最難降的妖魔是什么嚎朽？ 我笑而不...
    開封第一講書人閱讀 55,185評論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任，我火速辦了婚禮柬帕，結(jié)果婚禮上哟忍，老公的妹妹穿的比我還像新娘。我一直安慰自己陷寝，他們只是感情好锅很，可當(dāng)我...
    茶點故事閱讀 64,178評論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布。 她就那樣靜靜地躺著凤跑，像睡著了一般爆安。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上仔引，一...
    開封第一講書人閱讀 48,970評論 1贊 284
  • 城市分裂傳說
    那天扔仓，我揣著相機(jī)與錄音，去河邊找鬼咖耘。 笑死当辐，一個胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的鲤看。 我是一名探鬼主播，決...
    沈念sama閱讀 38,276評論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼耍群，長吁一口氣：“原來是場噩夢啊……” “哼义桂！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起蹈垢，我...
    開封第一講書人閱讀 36,927評論 0贊 259
  • 萬榮殺人案實錄
    序言：老撾萬榮一對情侶失蹤慷吊，失蹤者是張志新（化名）和其女友劉穎，沒想到半個月后曹抬，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體溉瓶，經(jīng)...
    沈念sama閱讀 43,400評論 1贊 300
  • ?護(hù)林員之死
    正文 獨居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 35,883評論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時候發(fā)現(xiàn)自己被綠了堰酿。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片疾宏。...
    茶點故事閱讀 37,997評論 1贊 333
  • 活死人
    序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖触创，靈堂內(nèi)的尸體忽然破棺而出坎藐，到底是詐尸還是另有隱情，我是刑警寧澤哼绑，帶...
    沈念sama閱讀 33,646評論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布岩馍，位于F島的核電站，受9級特大地震影響抖韩，放射性物質(zhì)發(fā)生泄漏蛀恩。R本人自食惡果不足惜，卻給世界環(huán)境...
    茶點故事閱讀 39,213評論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一茂浮、第九天 我趴在偏房一處隱蔽的房頂上張望双谆。 院中可真熱鬧，春花似錦励稳、人聲如沸佃乘。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,204評論 0贊 19
  • 一樁弒父案驹尼，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽趣避。三九已至，卻和暖如春新翎，著一層夾襖步出監(jiān)牢的瞬間程帕，已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 31,423評論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工地啰， 沒想到剛下飛機(jī)就差點兒被人妖公主榨干…… 1. 我叫王不留愁拭，地道東北人。 一個月前我還...
    沈念sama閱讀 45,423評論 2贊 352
  • 代替公主和親
    正文 我出身青樓亏吝，卻偏偏與公主長得像岭埠，于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子蔚鸥，可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 42,722評論 2贊 345