169_适配体SELEX分析¶
一句话概述¶
SELEX-seq(系统进化配体指数富集测序)通过高通量测序追踪随机DNA/RNA文库经过多轮筛选后的序列富集动态，利用生物信息学分析序列富集度、结构预测和亲和力排名，加速高亲和力适配体(Aptamer)的发现。
核心知识点总览¶
| 知识点 | 说明 |
|---|---|
| SELEX | 体外筛选高亲和力核酸适配体的方法 |
| 高通量SELEX(HT-SELEX) | 结合NGS追踪每轮筛选的序列变化 |
| 序列富集分析 | 比较不同轮次中序列频率变化 |
| 基序(Motif)发现 | 找到富集序列中的保守结构元件 |
| 二级结构预测 | 预测适配体的折叠结构(茎环/G-四链体等) |
| 亲和力排名 | 基于富集度和结构特征预测结合亲和力 |
| 聚类分析 | 将相似序列聚类为家族 |
| 应用场景 | 诊断探针、药物递送、生物传感器 |
各步骤详解¶
第一步：SELEX-seq实验和数据概述¶
白话解释： SELEX就像"钓鱼比赛"：从10^14种随机DNA序列的"大池塘"里，经过多轮"钓鱼"(与靶蛋白结合→洗掉不结合的→留下结合的→扩增→下一轮)，最终能紧紧咬住靶标的序列会越来越多。HT-SELEX在每一轮都做测序，看哪些序列在变多(被富集)——这些就是潜在的高亲和力适配体。
数据结构：
Round0 (初始文库): 10M unique sequences, 均匀分布
Round3: 部分序列开始富集
Round6: 少量序列高度富集
Round10: 若干序列主导(>1% each)
第二步：数据预处理和序列提取¶
代码示例：
# === 质控和接头去除 ===
# SELEX文库结构: 5'固定区-随机区(N20-N60)-3'固定区
fastp -i round3_R1.fq.gz -I round3_R2.fq.gz \
-o clean_R1.fq.gz -O clean_R2.fq.gz \
--adapter_sequence AGATCGGAAGAG
# === 提取随机区 ===
# 已知5'引物(20nt) + 随机区(40nt) + 3'引物(20nt)
cutadapt -g FIXED_5PRIME_SEQUENCE \
-a FIXED_3PRIME_SEQUENCE \
-o random_region.fq.gz \
--discard-untrimmed \
clean_R1.fq.gz
# === 序列计数和去重 ===
zcat random_region.fq.gz | awk 'NR%4==2' | sort | uniq -c | \
sort -rn > round3_counts.txt
# 输出: count sequence
# === Python处理 ===
from collections import Counter
from Bio import SeqIO
import pandas as pd
def process_selex_round(fastq_file, primer5="ATCGATCG", primer3="GCTAGCTA"):
"""提取随机区并计数"""
sequences = Counter()
for record in SeqIO.parse(fastq_file, "fastq"):
seq = str(record.seq)
# 找到引物位置并提取随机区
start = seq.find(primer5)
end = seq.find(primer3)
if start >= 0 and end > start:
random_region = seq[start+len(primer5):end]
if 35 <= len(random_region) <= 45: # 期望40nt
sequences[random_region] += 1
return sequences
# 处理各轮次
rounds = {}
for r in range(0, 11):
rounds[f"Round{r}"] = process_selex_round(f"round{r}_clean.fq")
print(f"Round{r}: {len(rounds[f'Round{r}'])} unique, {sum(rounds[f'Round{r}'].values())} total")
第三步：富集分析¶
代码示例：
import numpy as np
import pandas as pd
# === 计算各序列的富集度 ===
def calculate_enrichment(counts_later, counts_earlier, pseudocount=1):
"""计算两轮间的富集倍数"""
# 归一化为频率
total_later = sum(counts_later.values())
total_earlier = sum(counts_earlier.values())
enrichment = {}
all_seqs = set(counts_later.keys()) | set(counts_earlier.keys())
for seq in all_seqs:
freq_later = (counts_later.get(seq, 0) + pseudocount) / (total_later + pseudocount)
freq_earlier = (counts_earlier.get(seq, 0) + pseudocount) / (total_earlier + pseudocount)
enrichment[seq] = freq_later / freq_earlier
return enrichment
# 计算Round0→Round10的富集
enrichment_0_10 = calculate_enrichment(rounds["Round10"], rounds["Round0"])
# Top富集序列
top_sequences = sorted(enrichment_0_10.items(), key=lambda x: x[1], reverse=True)[:100]
# === 富集动力学追踪 ===
# 追踪top序列在各轮次的频率变化
top_seqs = [s[0] for s in top_sequences[:20]]
dynamics = pd.DataFrame(index=top_seqs)
for round_name, counts in rounds.items():
total = sum(counts.values())
dynamics[round_name] = [counts.get(seq, 0)/total for seq in top_seqs]
# 可视化富集曲线
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(10, 6))
for i, seq in enumerate(top_seqs[:5]):
ax.plot(range(11), dynamics.loc[seq].values, label=f"Seq_{i+1}")
ax.set_xlabel("SELEX Round")
ax.set_ylabel("Frequency")
ax.set_title("Aptamer Enrichment Dynamics")
ax.legend()
plt.savefig("enrichment_dynamics.png", dpi=150)
第四步：结构预测和亲和力排名¶
代码示例：
import RNA # ViennaRNA
# === 二级结构预测 ===
def predict_structure(sequence):
"""预测RNA/DNA适配体二级结构"""
structure, mfe = RNA.fold(sequence)
return structure, mfe
# 对top适配体预测结构
structural_features = []
for seq in top_seqs[:50]:
struct, mfe = predict_structure(seq)
# G-四链体检测
g4_potential = detect_g_quadruplex(seq)
structural_features.append({
'sequence': seq,
'structure': struct,
'MFE': mfe,
'stem_count': struct.count('('), # 配对碱基数
'loop_fraction': struct.count('.')/len(struct),
'G4_potential': g4_potential,
'GC_content': (seq.count('G')+seq.count('C'))/len(seq)
})
features_df = pd.DataFrame(structural_features)
# === G-四链体检测 ===
import re
def detect_g_quadruplex(seq):
"""检测G-四链体基序"""
pattern = r'(G{3,}).{1,7}(G{3,}).{1,7}(G{3,}).{1,7}(G{3,})'
match = re.search(pattern, seq)
return match is not None
# === 综合亲和力排名 ===
# 结合富集度和结构特征
features_df['enrichment'] = [enrichment_0_10[s] for s in features_df['sequence']]
features_df['affinity_score'] = (
0.5 * np.log2(features_df['enrichment']) +
0.3 * (-features_df['MFE']) + # 更稳定结构可能更好
0.2 * features_df['stem_count'] / features_df.sequence.str.len()
)
# 排名
features_df_sorted = features_df.sort_values('affinity_score', ascending=False)
print("Top 10 aptamer candidates:")
print(features_df_sorted[['sequence', 'enrichment', 'MFE', 'affinity_score']].head(10))
第五步：序列聚类和基序发现¶
代码示例：
from sklearn.cluster import DBSCAN
from Levenshtein import distance as lev_distance
import numpy as np
# === 基于编辑距离的聚类 ===
# 计算距离矩阵
n = len(top_seqs[:100])
dist_matrix = np.zeros((n, n))
for i in range(n):
for j in range(i+1, n):
d = lev_distance(top_seqs[i], top_seqs[j])
dist_matrix[i, j] = d
dist_matrix[j, i] = d
# DBSCAN聚类
clustering = DBSCAN(eps=5, min_samples=3, metric='precomputed')
clusters = clustering.fit_predict(dist_matrix)
print(f"发现 {len(set(clusters))-1} 个序列家族 (+噪声)")
# === 基序发现(MEME) ===
# 将top序列写入FASTA
with open("top_aptamers.fasta", "w") as f:
for i, seq in enumerate(top_seqs[:50]):
f.write(f">aptamer_{i}\n{seq}\n")
实战命令¶
# === SELEX-seq完整分析流程 ===
# 1. 质控
for r in 0 3 6 10; do
fastp -i round${r}.fq.gz -o round${r}_clean.fq.gz
done
# 2. 序列提取和计数
python extract_random_regions.py --rounds 0,3,6,10 --output counts/
# 3. 富集分析
python enrichment_analysis.py --input counts/ --output enrichment/
# 4. 结构预测
python structure_predict.py --top 100 --output structures/
# 5. 基序发现
meme top_aptamers.fa -dna -oc meme_out/ -nmotifs 5
面试常问点¶
Q1：HT-SELEX与传统SELEX的优势？¶
A： 传统SELEX需要8-15轮筛选后克隆测序少量序列。HT-SELEX在早期轮次(3-6轮)就可以通过NGS追踪百万级序列的富集动态：(1) 减少筛选轮数(避免过度富集导致多样性丢失)；(2) 发现中等亲和力但有独特结构的适配体；(3) 提供富集动力学信息帮助理解选择压力；(4) 可同时识别多个独立的适配体家族。
Q2：如何判断一个序列是否是真正的高亲和力适配体？¶
A： 多维度评估：(1) 富集度高(相对Round0频率增加>100倍)；(2) 富集动力学合理(逐轮递增而非突然出现)；(3) 预测结构稳定(负MFE)且有明确的结合口袋结构；(4) 与同家族序列有保守基序(功能性序列元件)；(5) 最终需要实验验证：SPR、ITC或滤膜结合实验测定Kd。
Q3：SELEX数据中PCR偏差如何影响结果？¶
A： PCR会引入两种偏差：(1) GC含量偏差——GC极端的序列扩增效率不同；(2) 结构偏差——强二级结构的序列可能扩增困难。影响：某些序列的富集可能是PCR偏好而非真实结合选择。应对：(1) 使用emulsion PCR或digital PCR减少偏差；(2) 分析中校正GC含量和结构复杂度；(3) 对照实验(无靶标SELEX)排除PCR偏好的序列。
Q4：G-四链体适配体有什么特殊性？¶
A： G-四链体(G4)是4条G-tract通过Hoogsteen氢键形成的特殊结构，在适配体中常见(如凝血酶适配体TBA)。特点：(1) 结构非常稳定(K+稳定)；(2) 序列模式：(G3+N1-7)4；(3) 靶标结合主要通过loop区和G-quartet平面；(4) 分析时需要使用专门的G4结构预测工具(QGRS Mapper)而非普通RNA折叠。
易错点¶
1. 未正确去除固定区序列¶
错误： 分析中包含引物序列，导致虚假保守性。 正确做法： 严格提取随机区，去除5'和3'固定序列。
2. 过度筛选导致多样性丢失¶
错误： 分析Round15的数据，只剩下1-2条高频序列。 正确做法： HT-SELEX的优势是在早期轮次(Round3-6)就可以看到富集趋势。过度筛选会丢失有价值的中等亲和力适配体。
3. 将PCR duplicates当作真实富集¶
错误： 不区分PCR复制和真实选择导致的频率增加。 正确做法： 使用UMI标记每个原始分子，或分析多个独立PCR重复的一致性。
4. 忽略负选择控制¶
错误： 不做negative control SELEX。 正确做法： 平行进行无靶标/非特异靶标SELEX，排除非特异富集的序列(如强PCR偏好、强自身结构但无结合的序列)。
补充知识¶
分析工具¶
- AptaSUITE: 综合SELEX-seq分析Java工具
- APTANI: 适配体鉴定R包
- FASTAptamer: 适配体序列分析Perl工具
- Mfold/ViennaRNA: 结构预测
- QGRS Mapper: G-四链体预测
适配体数据库¶
- Aptamer Base: 已知适配体及其靶标数据库
- SELEX DB: SELEX实验数据汇总