單個樣本分析接上一篇：Python圖文復現2022||02-數據介紹與下載，結合視頻觀看效果更佳~
視頻相關代碼如下：
- 視頻如下：https://www.youtube.com/watch?v=uvyG9yLuNSE
- 代碼：https://github.com/mousepixels/sanbomics_scripts
- 主代碼：https://github.com/mousepixels/sanbomics_scripts/blob/main/single_cell_analysis_complete_class.ipynb
讀取數據
先定義一個函數录肯，批量運行多個樣本：這里一定要注意縮進問題
def pp(csv_path):
adata = sc.read_csv(csv_path).T
sc.pp.filter_genes(adata, min_cells = 10)
sc.pp.highly_variable_genes(adata, n_top_genes = 2000, subset = True, flavor = 'seurat_v3')
scvi.model.SCVI.setup_anndata(adata)
vae = scvi.model.SCVI(adata)
vae.train()
solo = scvi.external.SOLO.from_scvi_model(vae)
solo.train()
df = solo.predict()
df['prediction'] = solo.predict(soft = False)
df.index = df.index.map(lambda x: x[:-2])
df['dif'] = df.doublet - df.singlet
doublets = df[(df.prediction == 'doublet') & (df.dif > 1)]
adata = sc.read_csv(csv_path).T
adata.obs['Sample'] = csv_path.split('_')[2] #'raw_counts/GSM5226574_C51ctr_raw_counts.csv'
adata.obs['doublet'] = adata.obs.index.isin(doublets.index)
adata = adata[~adata.obs.doublet]
sc.pp.filter_cells(adata, min_genes=200) #get rid of cells with fewer than 200 genes
#sc.pp.filter_genes(adata, min_cells=3) #get rid of genes that are found in fewer than 3 cells
adata.var['mt'] = adata.var_names.str.startswith('mt-') # annotate the group of mitochondrial genes as 'mt'
adata.var['ribo'] = adata.var_names.isin(ribo_genes[0].values)
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt', 'ribo'], percent_top=None, log1p=False, inplace=True)
upper_lim = np.quantile(adata.obs.n_genes_by_counts.values, .98)
adata = adata[adata.obs.n_genes_by_counts < upper_lim]
adata = adata[adata.obs.pct_counts_mt < 20]
adata = adata[adata.obs.pct_counts_ribo < 2]
return adata
這個過程比較久穆律，會依次讀取GSE171524數據集中的26例樣本，并進行上面的函數里面定義的分析。
這個過程中就可以跑去聽聽視頻了岭参。
import os
# 注意dir改成自己的路徑
dir = '/path/data/GSE171524/'
out = []
for file in os.listdir(dir+'raw_counts/'):
out.append(pp(dir + 'raw_counts/' + file))
將多個樣本連接合并在一起结笨，合并后共有105264個細胞：
adata = sc.concat(out)
adata
AnnData object with n_obs × n_vars = 105264 × 29236
obs: 'Sample', 'doublet', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'pct_counts_ribo'
var: 'n_cells'
過濾并保存：
sc.pp.filter_genes(adata, min_cells = 10)
from scipy.sparse import csr_matrix
adata.X = csr_matrix(adata.X)
adata.X
adata.write_h5ad(dir+'combined.h5ad')
預處理
先讀取上次保存的數據：105264的細胞
import scanpy as sc
import scvi
import seaborn as sns
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# 注意dir改成自己的路徑
dir = '/path/data/GSE171524/'
adata = sc.read_h5ad(dir+'combined.h5ad')
adata
AnnData object with n_obs × n_vars = 105264 × 29236
obs: 'Sample', 'doublet', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'pct_counts_ribo'
var: 'n_cells'
每個樣本中的各種指標統(tǒng)計：
adata.obs.groupby('Sample').count()
共26個樣本：
image-20221001224829602.png
低表達過濾以及數據標準化
sc.pp.filter_genes(adata, min_cells = 100)
adata.layers['counts'] = adata.X.copy()
sc.pp.normalize_total(adata, target_sum = 1e4)
sc.pp.log1p(adata)
adata.raw = adata
# 查看每個細胞的指標
adata.obs.head()
結果圖：
image-20221001225018337.png
去Sample間的批次以及聚類
這個過程運行時間也會有丟丟長
scvi.model.SCVI.setup_anndata(adata, layer = "counts", categorical_covariate_keys=["Sample"], continuous_covariate_keys=['pct_counts_mt', 'total_counts', 'pct_counts_ribo'])
model = scvi.model.SCVI(adata)
#may take a while without GPU
model.train()
adata.obsm['X_scVI'] = model.get_latent_representation()
adata.layers['scvi_normalized'] = model.get_normalized_expression(library_size = 1e4)
sc.pp.neighbors(adata, use_rep = 'X_scVI')
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution = 0.5)
sc.pl.umap(adata, color = ['leiden', 'Sample'], frameon = False)
plt.savefig(dir+"01-intergration_umap.png")
結果圖：
image-20221002000706224.png
保存數據：
adata.write_h5ad(dir + 'integrated.h5ad')
差異表達分析
使用沒有矯正的數據做差異表達分析
# 更改聚類數
sc.tl.leiden(adata, resolution = 1)
sc.tl.rank_genes_groups(adata, 'leiden')
# 可視化
sc.pl.rank_genes_groups(adata, n_genes=20, sharey=False)
plt.savefig(dir+"02-intergration_rank_genes_groups.png", dpi=300)
部分cluster的top基因：
image-20221002204205291.png
markers = sc.get.rank_genes_groups_df(adata, None)
markers = markers[(markers.pvals_adj < 0.05) & (markers.logfoldchanges > .5)]
markers
# 保存
markers.to_csv(dir+'markers.csv')
差異結果：
image-20221002002640069.png
使用模型標準化后的值做差異表達分析：
markers_scvi = model.differential_expression(groupby = 'leiden')
# 設定閾值
markers_scvi = markers_scvi[(markers_scvi['is_de_fdr_0.05']) & (markers_scvi.lfc_mean > .5)]
markers_scvi
# 保存
markers_scvi.to_csv(dir+'scvi_markers.csv')
結果：
image-20221002003124276.png
重新聚類后的結果可視化：總共得到34個cluster
sc.pl.umap(adata, color = ['leiden'], frameon = False, legend_loc = "on data")
plt.savefig(dir+"02-intergration_umap_1.png")
## 有多少cluster
adata.obs['leiden']
#Name: leiden, Length: 105264, dtype: category
#Categories (34, object): ['0', '1', '2', '3', ..., '30', '31', '32', '33']
結果圖：
image-20221002003157201.png
細胞類型注釋
文章中進行了三次注釋包晰，第一次注釋大類，主要為9個類：
image-20221002202142318.png
這幾個類的基因文章中沒有提供炕吸，就用我們自己收集的基因來進行注釋好了伐憾。
在視頻資源中，視頻speaker老師給了一張圖：是目前免疫細胞分類很詳細的一張圖了：
image-20221002215228373.png
基因可視化：
markers = ['EPCAM','KRT8','KRT18','KRT19', # epithelial cells
'CD68', 'CTSS', 'FCN1','CD163', # myeloid cells
'ACTA2', 'DCN', 'ACTB' , # fibroblasts
'PECAM1','CD34','VWF', # endothelial cells
'PTPRC','CD3D','CD3E','CD3G','CD8A', 'CD4', # T and natural killer (NK) lymphocytes
'CD79A', 'MS4A1','CD19', # B lymphocytes and plasma cells
'SNAP25', 'SYT1', # neuronal cells
'CPA3', # mast cells
'CST3', 'LAMP3', 'HLA-DQB2', 'HLA-DPB1', 'BIRC3', # antigen-presenting cells (APCs; primarily dendritic cells)
'MKI67','TOP2A', # Cycling
'HBB','HBD' # Erythroid
]
# 看某一個基因的表達情況
markers[markers.names=="HBB"]
markers[markers.group=="30"]
#, layer = 'scvi_normalized'
# , vmax = 5
sc.pl.umap(adata, color = ["HBB"], frameon = False, layer = 'scvi_normalized', vmax =4)
plt.savefig(dir+"03-intergration_umap_Erythroid_1.png")
一marker為基礎赫模，繪制如下類似圖树肃，vmax參數可以進行調整讓結果看起來更明顯，注視不出來的可以看看每個cluster高表達的基因瀑罗，查一查功能就立馬可以推斷出來了：
image-20221002003232078.png
結合marker表達以及cluster特異性高表達基因：
image-20221003010745511.png
細胞注釋如下：
for x in range(0,35):
print(f'"{x}":"", ')
cell_type = {"0":"T Cell",
"1":"Epithelial Cell",
"2":"Myeloid Cell",
"3":"Fibroblast",
"4":"Myeloid Cell",
"5":"T Cell",
"6":"Myeloid Cell",
"7":"Epithelial Cell",
"8":"Endothelial",
"9":"Fibroblast",
"10":"Myeloid Cell",
"11":"pDC",
"12":"Epithelial Cell",
"13":"Fibroblast",
"14":"Fibroblast",
"15":"Epithelial Cell",
"16":"Epithelial Cell",
"17":"Myeloid Cell",
"18":"Epithelial Cell",
"19":"Epithelial Cell",
"20":"Neuronal Cell",
"21":"Epithelial Cell",
"22":"B Cell",
"23":"Mast Cell",
"24":"T Cell",
"25":"pDC",
"26":"Myeloid Cell",
"27":"Endothelial",
"28":"Fibroblast",
"29":"Myeloid Cell",
"30":"Erythroid",
"31":"B Cell",
"32":"Neuronal Cell",
"33":"Myeloid Cell"
}
adata.obs['cell type'] = adata.obs.leiden.map(cell_type)
sc.pl.umap(adata, color = ['cell type'], frameon = False, legend_loc = "on data")
plt.tight_layout()
plt.savefig(dir+"04-intergration_umap_anno.png", dpi=300)
adata
adata.uns['scvi_markers'] = markers_scvi
adata.uns['markers'] = markers
# 保存
adata.write_h5ad(dir + 'integrated_anno.h5ad')
model.save(dir + 'model.model')
細胞注釋結果如下：
image-20221003011549203.png
這里注釋出來了一串文獻中沒有的紅細胞胸嘴。
下次進行詳細注釋~~~
