前言
线上课程教学
课题设计、定制生信分析
云服务器租赁
加微信备注99领取使用
跑代码时卡顿、电脑不给力让人抓狂!找果叔试用稳定高速的服务器,让分析顺畅无比!
代码学不会?bug 频繁出现,束手无策?实操生信分析课程赶快学起来!滴滴果叔领取体验课程哦~
线上课程教学
课题设计、定制生信分析
云服务器租赁
加微信备注99领取使用
pip install -r requirements.txt
pip install .
conda env create --name spatialcells --file=conda.yaml
pip install .
conda.yaml 指定如下:
name: spatial-cells-env
channels:
- conda-forge
- defaults
dependencies:
- python>=3.7, <3.11
- ipykernel
- matplotlib>=3.7.0
- pandas>=2.0.3
- seaborn>=0.12.2
- shapely>=2.0
- tqdm
- pip
- pip:
- anndata==0.9.2
- scanpy==1.9.4
加载相关的python包
import pandas as pd
import matplotlib.pyplot as plt
import anndata as ad
import spatialcells as spc
adata = ad.read_h5ad("../../data/MEL1_adata.h5ad")
spc.prep.setGate(adata, "KERATIN_cellRingMask", 6.4, debug=True)
spc.prep.setGate(adata, "SOX10_cellRingMask", 7.9, debug=True)
spc.prep.setGate(adata, "CD3D_cellRingMask", 7, debug=True)
marker = ["SOX10_cellRingMask_positive"]
communitycolumn = "COI_community"
ret = spc.spatial.getCommunities(adata, marker, eps=60, newcolumn=communitycolumn)
fig, ax = plt.subplots(figsize=(10, 8))
spc.plt.plotCommunities(
adata, ret, communitycolumn, plot_first_n_clusters=10, s=2, fontsize=10, ax=ax
)
ax.invert_yaxis()
plt.show()
communityIndexList = [6, 3, 14, 51, 29, 47, 39, 44, 22]
boundary = spc.spatial.getBoundary(
adata, communitycolumn, communityIndexList, alpha=130
)
boundary = spc.spa.pruneSmallComponents(boundary, min_edges=50, holes_min_edges=500)
roi_boundary = spc.spa.getExtendedBoundary(boundary, offset=2000)
markersize = 1
fig, ax = plt.subplots(figsize=(10, 7))
## all points
ax.scatter(
*zip(*adata.obs[["X_centroid", "Y_centroid"]].to_numpy()),
s=markersize,
color="grey",
alpha=0.2
)
# Points in selected commnities
xy = adata.obs[adata.obs[communitycolumn].isin(communityIndexList)][
["X_centroid", "Y_centroid"]
].to_numpy()
ax.scatter(xy[:, 0], xy[:, 1], s=markersize, color="r")
# Bounds of points in selected commnities
spc.plt.plotBoundary(boundary, ax=ax, label="Boundary", color="b")
spc.plt.plotBoundary(roi_boundary, ax=ax, label="ROI boundary", color="g")
ax.invert_yaxis()
plt.show()
spc.spatial.assignPointsToRegions(
adata,
[ ],
[ ],
assigncolumn="region",
default="BG",
)
point_size = 1
fig, ax = plt.subplots(figsize=(10, 7))
for region in sorted(set(adata.obs["region"])):
tmp = adata.obs[adata.obs.region == region]
ax.scatter(
*zip(*tmp[["X_centroid", "Y_centroid"]].to_numpy()),
s=point_size,
alpha=0.7,
label=region
)
spc.plt.plotBoundary(boundary, ax=ax, label="Boundary", color="purple")
spc.plt.plotBoundary(roi_boundary, ax=ax, label="ROI boundary", color="r")
plt.legend(loc="upper right")
ax.invert_yaxis()
plt.show()
def merge_pheno(row):
if row["phenotype_large_cohort"] in [
"T cells",
"Cytotoxic T cells",
"Exhausted T cells",
]:
return "T cells"
elif row["phenotype_large_cohort"] in ["Melanocytes"]:
return "Tumor cells"
else:
return "Other cells"
def cell_type(row):
if row["SOX10_cellRingMask_positive"]:
return "SOX10+"
elif row["CD3D_cellRingMask_positive"]:
return "CD3D+"
else:
return "Other cells"
# Applying the function to create the new columns
adata.obs["pheno1"] = pd.Categorical(adata.obs.apply(merge_pheno, axis=1))
adata.obs["Cell Types"] = pd.Categorical(adata.obs.apply(cell_type, axis=1))
spc.msmt.getRegionComposition(adata, "pheno1")
melano = adata[
(adata.obs.SOX10_cellRingMask_positive) & (adata.obs.region.isin(["Tumor_ROI", "Tumor"]))
]
tcells = adata[
(adata.obs.CD3D_cellRingMask_positive)
& (adata.obs.region.isin(["Tumor_ROI", "Tumor"]))
]
fig, ax = plt.subplots(figsize=(10, 7))
ax.invert_yaxis()
ax.set_aspect("equal")
plt.scatter(
tcells.obs["X_centroid"],
tcells.obs["Y_centroid"],
s=0.5,
label="T cells",
color="green",
alpha=0.5,
)
spc.plt.plotBoundary(roi_boundary, ax=ax, label="ROI boundary", color="r")
plt.legend(loc="upper right", markerscale=5)
plt.show()
tumor = adata[adata.obs.region.isin(["Tumor_ROI", "Tumor"])]
communitycolumn = "CD3D_cellRingMask_positive"
communityIndexList = [True]
immune_boundary = spc.spatial.getBoundary(
tumor, communitycolumn, communityIndexList, alpha=130
)
immune_boundary = spc.spa.pruneSmallComponents(
immune_boundary, min_edges=25, holes_min_edges=30, min_area=30000
)
markersize = 0.1
fig, ax = plt.subplots(figsize=(10, 7))
## all points
ax.scatter(
*zip(*adata.obs[["X_centroid", "Y_centroid"]].to_numpy()),
s=markersize,
color="grey",
alpha=0.2
)
# Points in selected commnities
xy = tumor.obs[tumor.obs[communitycolumn].isin(communityIndexList)][
["X_centroid", "Y_centroid"]
].to_numpy()
ax.scatter(xy[:, 0], xy[:, 1], s=markersize, color="green", alpha=1, label="T cells")
# Bounds of points in selected commnities
spc.plt.plotBoundary(
immune_boundary, ax=ax, label="Immune Cell Region Boundary", color="k", linewidth=1
)
spc.plt.plotBoundary(roi_boundary, ax=ax, label="ROI boundary", color="r")
# ax.set_xlim(0, 20000)
# ax.set_ylim(0, 13000)
ax.invert_yaxis()
ax.set_aspect("equal")
ax.set_axis_off()
# plt.legend(loc="upper right", markerscale=5, fontsize=13.5)
# plt.savefig("immune_cell_region1.png", dpi=400)
plt.show()
spc.spatial.assignPointsToRegions(
melano, [immune_boundary], ["T"], assigncolumn="tumor_isolated_region", default="F"
)
point_size = 0.5
fig, ax = plt.subplots(figsize=(10, 7))
ax.scatter(
*zip(*adata.obs[["X_centroid", "Y_centroid"]].to_numpy()),
s=markersize,
color="grey",
alpha=0.2
)
colors = ["red", "orange"]
labels = ["Immune-isolated Tumor Cells", "Immune-rich Tumor Cells"]
for i, region in enumerate(sorted(set(melano.obs["tumor_isolated_region"]))):
tmp = melano.obs[melano.obs.tumor_isolated_region == region]
ax.scatter(
*zip(*tmp[["X_centroid", "Y_centroid"]].to_numpy()),
s=point_size,
alpha=0.5,
color=colors[i],
label=labels[i]
)
# Bounds of points in selected commnities
spc.plt.plotBoundary(
immune_boundary, ax=ax, label="Immune Cell Region Boundary", color="k", linewidth=1
)
spc.plt.plotBoundary(roi_boundary, ax=ax, label="ROI boundary", color="b")
ax.invert_yaxis()
ax.set_axis_off()
# plt.savefig("roi_region1.png", dpi=400)
plt.show()
print("Percentage of tumor cells in immune-isolated regions: ")
melano.obs["tumor_isolated_region"].value_counts() / len(melano.obs)
roi_area = spc.msmt.getRegionArea(roi_boundary)
tumor_area = spc.msmt.getRegionArea(boundary)
immune_area = spc.msmt.getRegionArea(immune_boundary)
tumor_immune_overlap = boundary.intersection(immune_boundary)
overlap_area = spc.msmt.getRegionArea(tumor_immune_overlap)
print(f"Area of ROI: {roi_area:.2f}")
print(f"Area of main tumor cell region: {tumor_area:.2f}")
print(f"Area of immune cell region: {immune_area}")
print(f"Area of overlap between tumor and immune cell regions: {overlap_area:.2f}")
print(
f"Percentage of tumor region that has overlap with "
f"immune cell region: {overlap_area / tumor_area:.3f}"
)
识别与免疫细胞相邻的肿瘤细胞。
dists = spc.msmt.getMinCellTypesDistance(melano, tcells)
adata.obs.loc[
(adata.obs.SOX10_cellRingMask_positive) & (adata.obs.region.isin(["Tumor_ROI", "Tumor"])), "dist"
] = dists
threshold = 20
adata.obs["dist_binned"] = adata.obs["dist"] <= threshold
infiltrated = adata.obs[
(adata.obs.SOX10_cellRingMask_positive)
& (adata.obs.region == "Tumor")
& (adata.obs.dist_binned == True)
]
non_infiltrated = adata.obs[
(adata.obs.SOX10_cellRingMask_positive)
& (adata.obs.region == "Tumor")
& (adata.obs.dist_binned == False)
]
fig, ax = plt.subplots(figsize=(10, 6))
ax.invert_yaxis()
region = adata.obs[(adata.obs.region == "Tumor")]
plt.scatter(region["X_centroid"], region["Y_centroid"], s=1, alpha=0.2, color="grey")
plt.scatter(
infiltrated["X_centroid"],
infiltrated["Y_centroid"],
s=1,
alpha=0.5,
color="green",
label=f"Infiltrated (distance to t-cells <= {threshold}um)",
)
plt.legend(markerscale=10)
plt.show()
不会分析还想用生信工具助力发文咋办?有这顾虑的朋友,想一步到位就带着想法来,不论是代码实操还是在线文章结果复现,果叔照样能提供,还有大家都想要的服务器,找果叔获取就对了!
定制生信分析
服务器租赁
扫码咨询果叔
往期回顾
01 |
02 羡慕了!不做实验照样高分 “开挂”!中南大学雷光华团队玩转MR,3表2图成就1区7.6分佳绩!纯生信发文妙招你悟了吗?! |
03 |
04 |