背景
应该是CHS的项目，需要为二十个品种设计SSRs分子标记，数据是转录组数据，也就是说，如果要设计的话，中间应该是转录组作为输入。
参考文献
基于EST-SSR设计分子标记的参考
这个文章里写的很清楚，使用的是无参组装的转录本，然后他使用了一个工具：
https://github.com/xiaenhua/CandiSSR
老熟人啊我靠
没事，这何尝不是一种NTR
他这个环境依赖有点老了，不过他在他的github仓库中上传了全部依赖，这是很好的习惯，我下载下来后直接上传到服务器上。
这个环境太难配了，放弃，再找一篇文献，然后用更简单的方法配吧。
实操
Trinity组装
使用trinity_assemble_parallel.py这个脚本进行trinity组装
#! /usr/bin/env python3
"""
trinity_assemble_parallel.py
多样本并行Trinity转录组组装脚本
每个样本独立日志，支持单个Trinity多线程和多样本并行
"""
import os
import subprocess
import argparse
from multiprocessing import Pool
def Trinity_assemble_single(R1_file, R2_file, output_dir, threads_per_trinity=20, max_memory="20G"):
"""单样本Trinity组装，日志写入独立文件"""
prefix = os.path.basename(R1_file).split("_")[0]
trinity_out_dir = os.path.join(output_dir, prefix + "_Trinity_assemble")
os.makedirs(trinity_out_dir, exist_ok=True)
log_file = os.path.join(trinity_out_dir, f"{prefix}_trinity.log")
command = [
"Trinity",
"--seqType", "fq",
"--left", R1_file,
"--right", R2_file,
"--CPU", str(threads_per_trinity),
"--max_memory", max_memory,
"--output", trinity_out_dir
]
print(f"[{prefix}] 开始Trinity组装，日志: {log_file}")
with open(log_file, "w") as lf:
subprocess.run(command, stdout=lf, stderr=subprocess.STDOUT, check=True)
assemble_file = os.path.join(trinity_out_dir, "Trinity.fasta")
assemble_rename_file = os.path.join(trinity_out_dir, prefix + "_trinity_assemble.fasta")
os.rename(assemble_file, assemble_rename_file)
print(f"[{prefix}] Trinity完成，输出: {assemble_rename_file}")
return assemble_rename_file
# 顶层函数，multiprocessing 可调用
def run_sample(pair_and_params):
R1_file, R2_file, output_dir, threads_per_trinity, max_memory = pair_and_params
return Trinity_assemble_single(R1_file, R2_file, output_dir, threads_per_trinity, max_memory)
def Trinity_assemble_parallel(sample_file_list, output_dir, threads_per_trinity=20, max_memory="20G", max_parallel=2):
"""
多样本并行Trinity组装
sample_file_list: [(R1_file, R2_file), ...]
threads_per_trinity: 每个Trinity进程使用的线程数
max_memory: 每个Trinity进程使用的内存
max_parallel: 同时运行的Trinity进程数
"""
os.makedirs(output_dir, exist_ok=True)
# 构建 starmap 参数列表，每个元素是一个元组
params_list = [(R1, R2, output_dir, threads_per_trinity, max_memory) for (R1, R2) in sample_file_list]
with Pool(processes=max_parallel) as pool:
results = pool.map(run_sample, params_list)
return results
def read_index_file(index_file, delimiter='\t', col1=0, col2=1):
"""读取索引文件，返回样本(R1, R2)列表"""
sample_pairs = []
with open(index_file, 'r') as f:
for line in f:
if line.startswith("#") or not line.strip():
continue
cols = line.strip().split(delimiter)
if len(cols) <= max(col1, col2):
print(f"跳过行（列数不足）: {line.strip()}")
continue
sample_pairs.append((cols[col1], cols[col2]))
return sample_pairs
def main():
parser = argparse.ArgumentParser(description="多样本Trinity转录组并行组装工具")
parser.add_argument('-i', '--index_txt', required=True, help="索引文件，两列，R1和R2")
parser.add_argument('-o', '--output_dir', required=True, help="Trinity输出文件夹")
parser.add_argument('--threads_trinity', type=int, default=20, help="每个Trinity进程使用的线程数")
parser.add_argument('--max_parallel', type=int, default=2, help="同时运行的Trinity进程数")
parser.add_argument('--max_memory', type=str, default="50G", help="每个Trinity进程使用的内存")
args = parser.parse_args()
sample_pairs = read_index_file(args.index_txt)
print(f"找到 {len(sample_pairs)} 个样本对进行组装")
Trinity_assemble_parallel(
sample_pairs,
args.output_dir,
threads_per_trinity=args.threads_trinity,
max_memory=args.max_memory,
max_parallel=args.max_parallel
)
if __name__ == "__main__":
main()CD-hit去冗余
for f in *Trinity.fasta; do prefix=${f%.fasta}; cd-hit-est -i "$f" -o "${prefix}.cdhit95.fasta" -c 0.95 -n 10 -T 8 -M 0; doneMisa跑SSR
使用以下配置文件进行misa跑SSR
definition(unit_size,min_repeats): 1-20 3-6 3-5 4-5 5-5 6-4
interruptions(max_difference_between_2_SSRs): 100
GFF: false处理SSR输出文件
- 解析 misa 文件。
①-我希望从每个misa文件，得到一个矩阵，一列为SSR的序列，一列为这个SSR所对应的序列id，在母本序列上的起始和终止位置。
由于一个SSR可以反复多次出现，因此一个SSR可以有多个id，多个起始和终止位置。
②-由于每个样本都可以有这样的信息，因此我们可以得到一个新的表格，相当于再添加一列，即样本号，所以我们有样本，SSR序列，SSR序列所对应的id：起始和终止位置，可以吗？
类似于：
Sample SSR sequence_id1:SSR_start,SSR_end； sequence_id2:SSR_start,SSR_end；
也就是说，写一个解析当前路径下所有misa文件的代码，汇总成这些信息，文件名信息如下：
G30201_Trinity_assemble.Trinity.cdhit95.fasta.misa：G30201是样本id
# process_misa_file.py
#!/usr/bin/env python3
import os
import glob
import pandas as pd
input_dir = "."
files = glob.glob(os.path.join(input_dir, "*.misa"))
records = []
for f in files:
sample_id = os.path.basename(f).split("_")[0]
with open(f) as fh:
for line in fh:
line = line.strip()
# 跳过表头
if not line or line.startswith("ID") or line.startswith("SSR nr"):
continue
cols = line.split("\t")
if len(cols) < 7:
continue
seq_id = cols[0]
ssr_seq = cols[3]
start = cols[5]
end = cols[6]
records.append([
sample_id,
ssr_seq,
f"{seq_id}:{start}-{end}"
])
df = pd.DataFrame(records, columns=["Sample", "SSR_sequence", "Location"])
# =========================
# 核心：合并成一行
# =========================
df_merged = (
df.groupby(["Sample", "SSR_sequence"])["Location"]
.apply(lambda x: "; ".join(x))
.reset_index()
)
df_merged.to_csv("SSR_all_samples_merged.tsv", sep="\t", index=False)
print("Done -> SSR_all_samples_merged.tsv")B30108 (AAC)7 TRINITY_DN23_c1_g1_i1_len=2390_path=[0:0-2389]:1635-1655; TRINITY_DN1518_c0_g1_i1_len=2886_path=[0:0-890_1:891-1060_2:1061-1677_3:1678-1686_4:1687-2885]:286-306; TRINITY_DN1518_c0_g1_i3_len=2716_path=[0:0-890_2:891-1507_3:1508-1516_4:1517-2715]:286-306
B30108 (AACA)5 TRINITY_DN744_c0_g3_i1_len=2959_path=[0:0-1863_2:1864-2958]:96-115; TRINITY_DN3180_c0_g1_i6_len=8501_path=[0:0-573_2:574-584_4:585-1366_5:1367-1367_6:1368-1528_7:1529-1634_8:1635-2526_9:2527-2598_10:2599-7308_12:7309-7310_15:7311-8500]:6941-6960; TRINITY_DN3180_c0_g1_i9_len=8526_path=[0:0-573_2:574-584_3:585-585_4:586-1367_6:1368-1528_7:1529-1634_8:1635-2526_9:2527-2598_10:2599-7308_11:7309-7309_13:7310-7334_14:7335-7335_15:7336-8525]:6941-6960
B30108 (AACA)6 TRINITY_DN9750_c0_g2_i1_len=1143_path=[0:0-1142]:1055-1078
B30108 (AACAAA)4 TRINITY_DN4266_c0_g1_i2_len=2116_path=[1:0-624_2:625-938_3:939-1647_5:1648-1685_6:1686-1705_7:1706-2115]:1287-1310
B30108 (AACAGT)5 TRINITY_DN437_c1_g1_i4_len=1030_path=[2:0-353_3:354-453_5:454-482_6:483-1029]:1-30- 计算出样本唯一 SSR 。
好了，现在SSR_all_samples_merged.tsv这个文件中包含了所有的样本的SSR信息，我希望你只根据第一列和第二列，把样本特异存在的SSR，筛选出来，可以做到吗？也就是说，在这个表格基础上做一次过滤，剩下的全部都是样本唯一的SSR
(variety3) root@14dc642d5aaf:/project1/bin# cat process_misa_file2.py
#!/usr/bin/env python3
import pandas as pd
# 读入合并后的文件
df = pd.read_csv("SSR_all_samples_merged.tsv", sep="\t")
# =========================
# 关键步骤：判断 SSR 是否跨样本
# =========================
# 统计每个 SSR_sequence 出现在多少个 Sample 中
ssr_sample_count = df.groupby("SSR_sequence")["Sample"].nunique()
# 找出“只在一个样本中出现”的 SSR
unique_ssr = ssr_sample_count[ssr_sample_count == 1].index
# 过滤原表
df_unique = df[df["SSR_sequence"].isin(unique_ssr)].copy()
# 输出
df_unique.to_csv("SSR_sample_specific.tsv", sep="\t", index=False)
print("Done -> SSR_sample_specific.tsv")- 为样本唯一 SSR 设计两侧引物。
现在我们有了SSR_sample_specific.tsv这个文件后，我们需要设计SSR两侧引物
输入文件的格式是这样的：
(gmata_env) root@14dc642d5aaf:/project1/output/fasta_dir# head SSR_sample_specific.tsv
Sample SSR_sequence Location
B30108 (AAAGG)6 TRINITY_DN3841_c0_g1_i14_len=2724_path=[3:0-154_5:155-183_6:184-216_8:217-937_9:938-966_11:967-1043_14:1044-1111_16:1112-1206_17:1207-1616_18:1617-1825_19:1826-2723]:185-214; TRINITY_DN3841_c0_g1_i18_len=2515_path=[3:0-154_5:155-183_6:184-216_8:217-937_9:938-966_11:967-1043_14:1044-1111_16:1112-1206_17:1207-1616_19:1617-2514]:185-214; TRINITY_DN3841_c0_g1_i5_len=3231_path=[3:0-154_5:155-183_6:184-216_8:217-937_9:938-966_11:967-1043_13:1044-1220_14:1221-1288_15:1289-1618_16:1619-1713_17:1714-2123_18:2124-2332_19:2333-3230]:185-214
其中 TRINITY_DN3841_c0_g1_i14_len=2724_path=[3:0-154_5:155-183_6:184-216_8:217-937_9:938-966_11:967-1043_14:1044-1111_16:1112-1206_17:1207-1616_18:1617-1825_19:1826-2723]:185-214;这样一个部分中:以前代表序列id，:以后代表SSR在序列上的起始和终止位置，请注意，Location可以有多个sequence_id
这里我以这一行作为例子，sample为B30108，你需要找到对应样本B30108的fasta文件，
B30108_Trinity_assemble.Trinity.cdhit95.fasta（所有样本命名格式均相同），
然后找到序列TRINITY_DN3841_c0_g1_i14，根据坐标185-214设计引物。
一个样本的一个SSR序列可能要设计多个引物，比如这一行，你还要找到序列TRINITY_DN3841_c0_g1_i18和序列TRINITY_DN3841_c0_g1_i5。这里你要格外注意序列id的格式
设计引物的代码可以参考这部分：
=========================
FASTA dict cleanup
=========================
for k in list(fasta_dict.keys()):
seq_dict = SeqIO.to_dict(SeqIO.parse(fasta_dict[k], "fasta"))
fasta_dict[k] = {
kk.split(" ")[0]: vv for kk, vv in seq_dict.items()
}
=========================
helper functions
=========================
def gc(seq):
seq = seq.upper()
return (seq.count("G") + seq.count("C")) / len(seq) * 100
def tm(seq):
seq = seq.upper()
return 2*(seq.count("A")+seq.count("T")) + 4*(seq.count("G")+seq.count("C"))
def qc(f, r):
if "NA" in (f, r):
return False
if not (40 <= gc(f) <= 60 and 40 <= gc(r) <= 60):
return False
if abs(tm(f) - tm(r)) > 2:
return False
return True
=========================
primer design
=========================
results = []
for _, row in loc_df.iterrows():
sample = row["Sample"]
ssr = row["SSR"]
seq_id = row["SeqID"]
start = int(row["Start"])
end = int(row["End"])
if sample not in fasta_dict:
continue
if seq_id not in fasta_dict[sample]:
continue
full_seq = str(fasta_dict[sample][seq_id].seq)
left = max(0, start - FLANK - 1)
right = min(len(full_seq), end + FLANK)
flank = full_seq[left:right]
tstart = start - left
tlen = end - start + 1
try:
p = primer3.bindings.designPrimers(
{
'SEQUENCE_TEMPLATE': flank,
'SEQUENCE_TARGET': [tstart, tlen]
},
{
'PRIMER_OPT_SIZE': 20,
'PRIMER_MIN_SIZE': 18,
'PRIMER_MAX_SIZE': 25,
'PRIMER_PRODUCT_SIZE_RANGE': [[100, 300]],
'PRIMER_NUM_RETURN': 1
}
)
f = p.get("PRIMER_LEFT_0_SEQUENCE", "NA")
r = p.get("PRIMER_RIGHT_0_SEQUENCE", "NA")
prod = p.get("PRIMER_PAIR_0_PRODUCT_SIZE", "NA")
pass_qc = qc(f, r)
except:
f, r, prod = "NA", "NA", "NA"
pass_qc = False
results.append([
sample, ssr, seq_id, start, end,
f, r, prod, pass_qc
])
=========================
output
=========================
df = pd.DataFrame(results, columns=[
"Sample", "SSR", "SeqID", "Start", "End",
"Primer_F", "Primer_R", "Product_Size", "PASS"
])
请你给我完整的代码
# process_misa_file3.py
#!/usr/bin/env python3
import os
import re
import pandas as pd
from Bio import SeqIO
import primer3
# ==========================
# 参数设置（可以自己调）
# ==========================
FLANK = 200 # SSR两侧扩展长度（建议100~300）
# ==========================
# 读取SSR表
# ==========================
df = pd.read_csv("SSR_sample_specific.tsv", sep="\t")
# ==========================
# Step1：解析Location字段
# ==========================
records = []
for _, row in df.iterrows():
sample = row["Sample"]
ssr = row["SSR_sequence"]
locs = row["Location"].split(";")
for loc in locs:
loc = loc.strip()
if not loc:
continue
# 分割：seq_id 和 坐标
try:
seq_part, pos_part = loc.rsplit(":", 1)
except:
continue
# 🔥 关键：去掉 path，只保留 TRINITY ID
seq_id = seq_part.split("_len=")[0]
# 提取坐标
try:
start, end = map(int, pos_part.split("-"))
except:
continue
records.append([sample, ssr, seq_id, start, end])
loc_df = pd.DataFrame(records, columns=[
"Sample", "SSR", "SeqID", "Start", "End"
])
# ==========================
# Step2：加载所有fasta
# ==========================
fasta_dict = {}
samples = loc_df["Sample"].unique()
for s in samples:
fasta_file = f"{s}_Trinity_assemble.Trinity.cdhit95.fasta"
if not os.path.exists(fasta_file):
print(f"[WARN] missing fasta: {fasta_file}")
continue
seq_dict = SeqIO.to_dict(SeqIO.parse(fasta_file, "fasta"))
# 🔥 去掉 header 后面的乱七八糟信息
fasta_dict[s] = {
k.split(" ")[0]: v for k, v in seq_dict.items()
}
# ==========================
# Step3：辅助函数
# ==========================
def gc(seq):
seq = seq.upper()
return (seq.count("G") + seq.count("C")) / len(seq) * 100 if len(seq) > 0 else 0
def tm(seq):
seq = seq.upper()
return 2*(seq.count("A")+seq.count("T")) + 4*(seq.count("G")+seq.count("C"))
def qc(f, r):
if "NA" in (f, r):
return False
if not (40 <= gc(f) <= 60 and 40 <= gc(r) <= 60):
return False
if abs(tm(f) - tm(r)) > 2:
return False
return True
# ==========================
# Step4：设计引物
# ==========================
results = []
for _, row in loc_df.iterrows():
sample = row["Sample"]
ssr = row["SSR"]
seq_id = row["SeqID"]
start = int(row["Start"])
end = int(row["End"])
if sample not in fasta_dict:
continue
if seq_id not in fasta_dict[sample]:
# 🔥 debug关键
print(f"[MISS] {sample} {seq_id}")
continue
full_seq = str(fasta_dict[sample][seq_id].seq)
# 取flanking区域
left = max(0, start - FLANK - 1)
right = min(len(full_seq), end + FLANK)
flank = full_seq[left:right]
# SSR在flank中的位置
tstart = start - left - 1
tlen = end - start + 1
try:
p = primer3.bindings.designPrimers(
{
'SEQUENCE_TEMPLATE': flank,
'SEQUENCE_TARGET': [tstart, tlen]
},
{
'PRIMER_OPT_SIZE': 20,
'PRIMER_MIN_SIZE': 18,
'PRIMER_MAX_SIZE': 25,
'PRIMER_PRODUCT_SIZE_RANGE': [[100, 300]],
'PRIMER_NUM_RETURN': 1
}
)
f = p.get("PRIMER_LEFT_0_SEQUENCE", "NA")
r = p.get("PRIMER_RIGHT_0_SEQUENCE", "NA")
prod = p.get("PRIMER_PAIR_0_PRODUCT_SIZE", "NA")
pass_qc = qc(f, r)
except Exception as e:
f, r, prod = "NA", "NA", "NA"
pass_qc = False
results.append([
sample, ssr, seq_id, start, end,
f, r, prod, pass_qc
])
# ==========================
# Step5：输出
# ==========================
out_df = pd.DataFrame(results, columns=[
"Sample", "SSR", "SeqID", "Start", "End",
"Primer_F", "Primer_R", "Product_Size", "PASS"
])
out_df.to_csv("SSR_primers.tsv", sep="\t", index=False)
print("Done -> SSR_primers.tsv")只保留计算通过的标记和简单的SSR重复
awk 'NR==1 || ($9=="True" && split($2,a,"(")==2)' SSR_primers.tsv > SSR_primers.filtered.tsv去重复，去掉两侧引物完全相同，SSR也相同的行
awk 'NR==1 {print; next} !seen[$2 FS $6 FS $7]++' SSR_primers.filtered.tsv > SSR_primers.unique.tsv26年5月29日记录
为了节约空间，删除了所有trinity组装文件