如何从AnnData数据集中筛选3'UTR含特定RNA序列的基因

下面的流程可以直接对接你用scanpy处理好的AnnData对象，不需要手动处理基因组注释、坐标转换这些生信操作，按步骤改对应参数就能跑。

1. 安装依赖库 #

先在终端运行命令装需要的包：

pip install scanpy biopython pyensembl pandas

第一次运行前需要加载对应物种的基因组注释，以人类hg38、Ensembl 109版本为例，小鼠就把species改成"mouse"、release选对应mm10的版本就行：

import pyensembl
# 替换成你自己研究物种对应的注释参数
ensembl_data = pyensembl.EnsemblRelease(release=109, species="human")
ensembl_data.index()

这步第一次跑会自动下载对应版本的注释文件，等进度条跑完就行，不用手动找资源下载。

2. 批量提取基因3'UTR序列 #

先读入你自己的AnnData文件，遍历里面的所有基因提取最长转录本的3'UTR序列（选最长转录本是为了覆盖所有可能的结合位点，避免不同转录本剪接带来的漏检）：

import scanpy as sc
import pandas as pd

# 替换成你自己的h5ad文件存储路径
adata = sc.read_h5ad("./your_scrna_result.h5ad")
all_genes = adata.var_names.tolist()

# 替换成你要找的目标RNA序列，比如基序、miRNA结合位点序列，直接写带U的RNA序列就行
target_rna_seq = "YOUR_TARGET_SEQUENCE_HERE"
# 转成DNA形式的序列（U替换成T）和基因组提取的序列做匹配
target_dna_seq = target_rna_seq.replace("U", "T").upper()

gene_utr_map = {}
for gene in all_genes:
    try:
        # 拿到该基因对应的所有转录本
        tx_list = ensembl_data.transcripts_by_name(gene)
        all_utr = []
        for tx in tx_list:
            utr_seq = tx.three_prime_utr_sequence
            if utr_seq:
                all_utr.append(str(utr_seq).upper())
        # 保留最长的3'UTR序列
        if all_utr:
            gene_utr_map[gene] = max(all_utr, key=len)
    except:
        # 注释库中匹配不到的基因直接跳过，不影响整体流程
        continue

3. 筛选目标基因并对接后续分析 #

直接做精确序列匹配，把结果直接写回AnnData的var元数据里，和你之前的scanpy分析流程完全兼容：

# 筛选3'UTR包含目标序列的基因
matched_gene_list = [g for g, seq in gene_utr_map.items() if target_dna_seq in seq]

# 把匹配结果写到adata的属性里，后续分群、差异分析直接调用
adata.var["is_target_gene"] = adata.var_names.isin(matched_gene_list)

# 可选：把匹配到的基因列表导出成表格
pd.Series(matched_gene_list, name="target_genes").to_csv("./matched_3utr_genes.csv", index=False)

# 后续如果要取这些基因的子集做分析，直接用下面的代码就行，和你之前的scanpy操作没有区别
# adata_target = adata[:, adata.var["is_target_gene"]]

常见问题处理 #

• 如果你adata里存的基因名是Ensembl ID（比如ENSG开头的ID）而不是基因名，把上面代码里的transcripts_by_name换成transcripts_by_id即可
• 如果你要找的是带简并碱基的基序（比如存在可变碱基的序列模式），把字符串精确匹配换成正则表达式匹配就可以，不需要额外改其他逻辑
• 如果自动下载注释速度慢，可以手动把对应版本的GTF压缩包放到用户目录下的.pyensembl缓存文件夹，程序会自动识别，不需要重新下载
• 如果你研究的是其他非模式物种，只要有对应物种的Ensembl注释版本，把species和release参数改成对应值就能跑，流程是通用的

内容的提问来源于stack exchange，提问作者Greenline