作者，Evil Genius

生活總是無情，牛馬也很無奈，除了接受現(xiàn)實(shí)，別無其他辦法。

公司大調(diào)整的日子，也進(jìn)入倒計(jì)時(shí)。

不過還是推薦大家走好每一步，生物信息這個(gè)行業(yè)不算什么暴利行業(yè)，但是有個(gè)相對(duì)安穩(wěn)的生活還是能做到的，將來大家畢業(yè)參加工作，比如進(jìn)了諾禾、華大等，生活還是不成問題的。

但是隨著現(xiàn)在越來越卷，公司的進(jìn)入門檻也在提高，18年我參加工作基本上是個(gè)碩士就行了，不要求有什么代碼語言基礎(chǔ)，來了公司會(huì)培訓(xùn)，25年不行了，不是科班出身，沒有代碼基礎(chǔ)，沒有相關(guān)的一些基礎(chǔ)，很難進(jìn)去。

以后的路怎么走，誰也不知道，最好自己再說吧。

visium分析思路

今天我們更新腳本

這個(gè)圖片在課上講過，強(qiáng)調(diào)的是在不同條件下，配受體共定位分布差異，那么今天我們需要更新的是，借助單細(xì)胞 + 空間的多組學(xué)數(shù)據(jù)，分析不同條件下的共定位配受體，并且分析條件之間差異的共定位配受體，23年培訓(xùn)的共定位繪圖大家有空復(fù)習(xí)一下。

所有的分析都建立在一定的基礎(chǔ)之上，linux 、R、python生信必備，雖然讓大家好好學(xué)習(xí)這個(gè)，11月打算開番外，2025番外--linux、R、python培訓(xùn), 不過實(shí)際情況也能想到，沒多少人把基礎(chǔ)當(dāng)回事，用的時(shí)候才知道著急了。

分兩步進(jìn)行

第一步，單細(xì)胞分析不同組別的差異配受體對(duì)。

第二步，首先驗(yàn)證同組間配受體的空間共定位情況，去掉單細(xì)胞分析得到共定位較差的配受體對(duì)，然后分析不同組間共定位差異的配受體對(duì)，進(jìn)行空間分析與展示，這個(gè)時(shí)候基本上就找到了靶點(diǎn)。

首先實(shí)現(xiàn)第一步，實(shí)現(xiàn)單細(xì)胞不同條件下的差異配受體對(duì)，做好細(xì)胞類型，放在adata.obs['celltype']下面，不同組別放在adata.obs['batch']下面。

import scanpy as sc
import squidpy as sq
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import mannwhitneyu
from statsmodels.stats.multitest import multipletests
import warnings
import sys
import os
warnings.filterwarnings('ignore')

# 設(shè)置參數(shù)
input_file = sys.argv[1] if len(sys.argv) > 1 else "data.h5ad"
output_dir = "./lr_analysis_results"
os.makedirs(output_dir, exist_ok=True)

print("=== 單細(xì)胞差異配受體對(duì)分析 ===")

# 1. 加載數(shù)據(jù)
print(f"\n1. 加載數(shù)據(jù): {input_file}")
adata = sc.read(input_file)

print(f"數(shù)據(jù)維度: {adata.n_obs} 個(gè)細(xì)胞, {adata.n_vars} 個(gè)基因")
print(f"批次分布: {adata.obs['batch'].value_counts().to_dict()}")
print(f"細(xì)胞類型分布: {adata.obs['celltype'].value_counts().to_dict()}")

# 2. 使用squidpy的ligrec分析配受體相互作用
print("\n2. 使用squidpy ligrec分析配受體相互作用...")

# 第一步：分析每個(gè)批次內(nèi)的配受體相互作用
batches = adata.obs['batch'].unique()
print(f"分析批次: {batches}")

# 存儲(chǔ)每個(gè)批次的結(jié)果
batch_results = {}

for batch in batches:
    print(f"\n分析批次: {batch}")
    batch_adata = adata[adata.obs['batch'] == batch].copy()
    
    # 使用squidpy的ligrec分析配受體相互作用
    sq.gr.ligrec(
        batch_adata,
        cluster_key="celltype",
        copy=True,
        use_raw=False,
        n_perms=1000
    )
    
    # 獲取結(jié)果
    if 'ligrec' in batch_adata.uns:
        lr_results = batch_adata.uns['ligrec']
        batch_results[batch] = lr_results
        
        # 篩選顯著的結(jié)果 (pvalue < 0.05)
        significant_pairs = lr_results[lr_results['pvalue'] < 0.05]
        print(f"批次 {batch} 中發(fā)現(xiàn) {len(significant_pairs)} 個(gè)顯著配受體對(duì)")
        
        # 保存結(jié)果
        significant_pairs.to_csv(f"{output_dir}/significant_lr_pairs_{batch}.csv", index=False)
        
        # 顯示top 10顯著對(duì)
        if len(significant_pairs) > 0:
            top_pairs = significant_pairs.nsmallest(10, 'pvalue')
            print(f"批次 {batch} Top 10 顯著配受體對(duì):")
            for idx, row in top_pairs.iterrows():
                print(f"  {row['ligand']} - {row['receptor']}: pvalue={row['pvalue']:.4f}")

# 3. 比較不同批次間的差異配受體
print("\n3. 比較不同批次間的差異配受體...")

# 收集所有批次中出現(xiàn)的配受體對(duì)
all_lr_pairs = set()
for batch, results in batch_results.items():
    if len(results) > 0:
        pairs = set(zip(results['ligand'], results['receptor']))
        all_lr_pairs.update(pairs)

print(f"總共發(fā)現(xiàn) {len(all_lr_pairs)} 個(gè)獨(dú)特的配受體對(duì)")

# 分析配受體對(duì)在不同批次中的出現(xiàn)情況
lr_presence = {}
for ligand, receptor in all_lr_pairs:
    pair_name = f"{ligand}_{receptor}"
    presence = []
    for batch in batches:
        if batch in batch_results and len(batch_results[batch]) > 0:
            batch_pairs = set(zip(batch_results[batch]['ligand'], batch_results[batch]['receptor']))
            present = 1 if (ligand, receptor) in batch_pairs else 0
        else:
            present = 0
        presence.append(present)
    lr_presence[pair_name] = presence

# 創(chuàng)建出現(xiàn)矩陣
presence_df = pd.DataFrame(lr_presence, index=batches).T

print("\n配受體對(duì)在不同批次中的出現(xiàn)情況:")
print(presence_df)

# 4. 識(shí)別批次特異的配受體對(duì)
print("\n4. 識(shí)別批次特異的配受體對(duì)...")

batch_specific_pairs = {}
for batch in batches:
    batch_specific = []
    for pair_name in presence_df.index:
        row = presence_df.loc[pair_name]
        # 只在當(dāng)前批次出現(xiàn)，其他批次不出現(xiàn)
        if row[batch] == 1 and row.sum() == 1:
            batch_specific.append(pair_name)
    
    batch_specific_pairs[batch] = batch_specific
    print(f"批次 {batch} 特異的配受體對(duì): {len(batch_specific_pairs[batch])} 個(gè)")

# 5. 識(shí)別共同出現(xiàn)的配受體對(duì)
print("\n5. 識(shí)別共同出現(xiàn)的配受體對(duì)...")
common_pairs = []
for pair_name in presence_df.index:
    row = presence_df.loc[pair_name]
    if row.sum() == len(batches):  # 在所有批次中都出現(xiàn)
        common_pairs.append(pair_name)

print(f"在所有批次中都出現(xiàn)的配受體對(duì): {len(common_pairs)} 個(gè)")

# 6. 可視化結(jié)果
print("\n6. 可視化分析結(jié)果...")

# 創(chuàng)建熱圖顯示配受體對(duì)的出現(xiàn)情況
plt.figure(figsize=(12, 8))
sns.heatmap(presence_df, cmap="YlGnBu", cbar_kws={'label': '出現(xiàn) (1) / 不出現(xiàn) (0)'})
plt.title("配受體對(duì)在不同批次中的出現(xiàn)情況")
plt.tight_layout()
plt.savefig(f"{output_dir}/lr_presence_heatmap.png", dpi=300, bbox_inches='tight')
plt.show()

# 創(chuàng)建批次特異性配受體對(duì)的可視化
plt.figure(figsize=(10, 6))
batch_counts = [len(pairs) for pairs in batch_specific_pairs.values()]
plt.bar(batch_specific_pairs.keys(), batch_counts)
plt.title("各批次特異的配受體對(duì)數(shù)量")
plt.ylabel("特異配受體對(duì)數(shù)量")
plt.xticks(rotation=45)
plt.tight_layout()
plt.savefig(f"{output_dir}/batch_specific_lr_pairs.png", dpi=300, bbox_inches='tight')
plt.show()

# 7. 保存匯總結(jié)果
print("\n7. 保存匯總結(jié)果...")

# 創(chuàng)建匯總表格
summary_data = []
for batch in batches:
    if batch in batch_results and len(batch_results[batch]) > 0:
        total_pairs = len(batch_results[batch])
        significant_pairs = len(batch_results[batch][batch_results[batch]['pvalue'] < 0.05])
        specific_pairs = len(batch_specific_pairs[batch])
    else:
        total_pairs = 0
        significant_pairs = 0
        specific_pairs = 0
    
    summary_data.append({
        'batch': batch,
        'total_lr_pairs': total_pairs,
        'significant_pairs': significant_pairs,
        'batch_specific_pairs': specific_pairs
    })

summary_df = pd.DataFrame(summary_data)
summary_df.to_csv(f"{output_dir}/lr_analysis_summary.csv", index=False)
print("匯總結(jié)果:")
print(summary_df)

# 保存共同出現(xiàn)的配受體對(duì)
if common_pairs:
    common_df = pd.DataFrame({'common_lr_pairs': common_pairs})
    common_df.to_csv(f"{output_dir}/common_lr_pairs.csv", index=False)

# 保存批次特異的配受體對(duì)
for batch, pairs in batch_specific_pairs.items():
    if pairs:
        specific_df = pd.DataFrame({f'{batch}_specific_pairs': pairs})
        specific_df.to_csv(f"{output_dir}/{batch}_specific_lr_pairs.csv", index=False)

這樣我們就得到了單細(xì)胞不同條件下，差異的配受體對(duì)。

接下來第二步，分析差異性配受體對(duì)的空間模式，是否空間存在共定位，空間分布是否有差異（臨近包括一個(gè)spot臨近的6個(gè)spot）