commit e1b31567545330a9f0fc69034d00ff805abcf810 Author: Lihatoo <1747565629@gmail.com> Date: Thu May 28 01:04:44 2026 +0800 Optimize nupack service startup and compute diff --git a/.dockerignore b/.dockerignore new file mode 100644 index 0000000..0cd8183 --- /dev/null +++ b/.dockerignore @@ -0,0 +1,19 @@ +.git +.gitignore +__pycache__ +*.pyc +*.pyo +*.pyd +.pytest_cache +.mypy_cache +.DS_Store +other +nupack-history-*.json +NUPACK.html +NUPACK_files +1.png +nupack/vendor/nupack-4.0.2.0/source +nupack/vendor/nupack-4.0.2.0/package/*.whl +!nupack/vendor/nupack-4.0.2.0/package/nupack-4.0.2.0-cp312-cp312-linux_x86_64.whl +rna/* +!rna/ViennaRNA-2.7.2.tar.gz diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..25bfe7d --- /dev/null +++ b/.gitignore @@ -0,0 +1,18 @@ +__pycache__/ +*.py[cod] +.pytest_cache/ +.mypy_cache/ +.DS_Store + +# Local exports and captured reference site assets are not needed to rebuild the service. +nupack-history-*.json +NUPACK.html +NUPACK_files/ +1.png + +# The Docker build uses the vendored cp312 Linux x86_64 wheel plus the ViennaRNA tarball. +rna/ViennaRNA-2.7.2/ +rna/nupack-4.0.2.0.zip +nupack/vendor/nupack-4.0.2.0/source/ +nupack/vendor/nupack-4.0.2.0/package/*.whl +!nupack/vendor/nupack-4.0.2.0/package/nupack-4.0.2.0-cp312-cp312-linux_x86_64.whl diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..21ad0d6 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,48 @@ +FROM python:3.12-slim + +ENV PYTHONDONTWRITEBYTECODE=1 +ENV PYTHONUNBUFFERED=1 +ENV NP_HOST=0.0.0.0 +ENV NP_PORT=18765 +ENV ENABLE_RNAPLOT=1 + +WORKDIR /app + +RUN apt-get update && apt-get install -y --no-install-recommends \ + autoconf \ + automake \ + build-essential \ + libgsl-dev \ + libtool \ + pkg-config \ + perl \ + && rm -rf /var/lib/apt/lists/* + +COPY nupack/vendor/nupack-4.0.2.0/package /tmp/nupack-package +COPY rna/ViennaRNA-2.7.2.tar.gz /tmp/ViennaRNA-2.7.2.tar.gz + +RUN cd /tmp \ + && tar -xzf ViennaRNA-2.7.2.tar.gz \ + && cd ViennaRNA-2.7.2 \ + && ./configure --without-python --without-perl \ + && make -j"$(nproc)" \ + && make install \ + && ldconfig \ + && rm -rf /tmp/ViennaRNA-2.7.2 /tmp/ViennaRNA-2.7.2.tar.gz + +RUN python -m pip install --no-cache-dir \ + numpy \ + scipy \ + pandas \ + pyyaml \ + jinja2 \ + /tmp/nupack-package/nupack-4.0.2.0-cp312-cp312-linux_x86_64.whl + +RUN python -m pip install --no-cache-dir redis + +WORKDIR /app/service +COPY service /app/service + +EXPOSE 18765 + +CMD ["python", "server.py"] diff --git a/README.md b/README.md new file mode 100644 index 0000000..5e48d1a --- /dev/null +++ b/README.md @@ -0,0 +1,56 @@ +# NP Replica + +This project provides a local Chinese web replica of the core NUPACK `analysis/input` workflow. + +## Run + +```bash +docker compose up --build -d +``` + +Open `http://127.0.0.1:18765`. + +## Reverse Proxy + +The service listens on port `18765`. You can later reverse proxy `np.lihato.icu` to this port from your public server. + +## Endpoints + +- `GET /` +- `GET /health` +- `GET /api/example` +- `POST /api/analyze` +- `POST /api/jobs` +- `GET /api/jobs/` + +## Notes + +- `tube` mode supports strand concentrations and `SetSpec(max_size=...)`. +- `complex` mode supports explicit complex definitions such as `A+B`. +- Pair matrices are returned as a preview block to keep responses manageable. +- The Docker image compiles ViennaRNA `RNAplot` from `rna/ViennaRNA-2.7.2.tar.gz` and embeds SVG structure plots into MFE results. +- Multistrand MFE plots now use a split-strand layout by default when a structure contains `+`, with `RNAplot` retained as the fallback path. +- The web UI now supports Chinese and English switching, browser-side history, result export, and async polling for submitted jobs. +- Docker now starts three services: `np-replica`, `np-worker`, and `redis`. +- Async jobs use Redis-backed storage and queueing when `NP_REDIS_URL` is configured; otherwise the app falls back to the older in-memory mode. +- The browser-facing API and page behavior stay the same after enabling Redis-backed jobs. +- Worker concurrency is controlled by `NP_WORKER_CONCURRENCY` and defaults to `2`. +- Per-job native thread usage is capped by `NP_PER_JOB_THREAD_LIMIT` and defaults to `1` to prevent a single analysis from monopolizing the machine. +- `NP_STRUCTURE_PLOT_MODE` supports `auto` (default), `split`, and `rnaplot`. +- Static pages are served with ETag-based browser caching; JSON API responses remain uncached. +- The force-directed structure viewer lazy-loads D3 only when that view is opened, so the initial page load is not blocked by the external CDN. +- `/health` uses Redis queue length plus a running-job set instead of scanning every historical job on each poll. + + + + 现在 `docker-compose.yml` 面向 64 核 / 64G WSL 服务器的并发策略是： + + - 最多同时跑 16 个任务 + - 每个任务最多用 4 个计算线程 + - 后续任务进入队列等待 + - 理论上最多占用约 64 个 native 计算线程，避免 16 × 8 这类过度超卖 + + 如果后面你发现机器还会被压满，最直接的调法就是在 docker-compose.yml 里继续压： + + - 把 NP_WORKER_CONCURRENCY 改小 + - 或保持并发不变，把 NP_PER_JOB_THREAD_LIMIT 改成 2 或 1 diff --git a/design(2).txt b/design(2).txt new file mode 100644 index 0000000..0ba5259 --- /dev/null +++ b/design(2).txt @@ -0,0 +1,30 @@ +from nupack import * + +my_model = Model(material='dna', celsius=37,sodium=0.1, magnesium=0.02) + +c = Domain('N8', name='c') +cc = Domain('N3', name='cc') +cc2 = Domain('N3', name='cc2') +d = Domain('N8', name='d') +e = Domain('N7R1', name='e') +f = Domain('N10', name='f') +g = Domain('N8', name='g') +k = Domain('N8', name='k') +h = Domain('N15R1', name='h') +hh = Domain('N6', name='hh') +E = Domain('GGCTAGCTACAACGA', name='E') + +E1 = TargetStrand([h, E, ~e, ~c], name='Strand E1') +F1 = TargetStrand([~f, ~d, ~g, ~k, f, c, e, hh], name='Strand F1') + +E1F1 = TargetComplex([E1, F1], '.10(6.15(16+(10.24)32', name='E1F1') + +t1 = TargetTube(on_targets={E1F1: 1e-8}, name='t1', + off_targets=SetSpec(max_size=4)) + +pattern = Pattern(['A4', 'C4', 'G4', 'U4', 'T4']) +my_tubes = [t1,t2] +my_design = tube_design(tubes=my_tubes, + hard_constraints=[pattern],model=my_model) +my_result = my_design.run(trials=1) +print(my_result) \ No newline at end of file diff --git a/docker-compose.yml b/docker-compose.yml new file mode 100644 index 0000000..a27f57a --- /dev/null +++ b/docker-compose.yml @@ -0,0 +1,56 @@ +services: + redis: + image: redis:7.4.8-bookworm + container_name: np-redis + restart: always + mem_limit: 2g + command: ["redis-server", "--save", "", "--appendonly", "no"] + + np-replica: + image: np-app:v4.2 + build: + context: . + dockerfile: Dockerfile + container_name: np-replica + restart: always + mem_limit: 4g + depends_on: + - redis + ports: + - "18765:18765" + environment: + NP_HOST: 0.0.0.0 + NP_PORT: 18765 + NP_RUN_MODE: server + NP_REDIS_URL: redis://redis:6379/0 + NP_JOB_TTL_SECONDS: 3600 + NP_WORKER_CONCURRENCY: 16 + NP_PER_JOB_THREAD_LIMIT: 4 + NP_NUPACK_CACHE_GB: 8.0 + OMP_NUM_THREADS: 4 + OPENBLAS_NUM_THREADS: 4 + MKL_NUM_THREADS: 4 + NUMEXPR_NUM_THREADS: 4 + VECLIB_MAXIMUM_THREADS: 4 + GOTO_NUM_THREADS: 4 + + np-worker: + image: np-app:v4.2 + container_name: np-worker + restart: always + mem_limit: 56g + depends_on: + - redis + environment: + NP_RUN_MODE: worker + NP_REDIS_URL: redis://redis:6379/0 + NP_JOB_TTL_SECONDS: 3600 + NP_WORKER_CONCURRENCY: 16 + NP_PER_JOB_THREAD_LIMIT: 4 + NP_NUPACK_CACHE_GB: 8.0 + OMP_NUM_THREADS: 4 + OPENBLAS_NUM_THREADS: 4 + MKL_NUM_THREADS: 4 + NUMEXPR_NUM_THREADS: 4 + VECLIB_MAXIMUM_THREADS: 4 + GOTO_NUM_THREADS: 4 diff --git a/need.md b/need.md new file mode 100644 index 0000000..c0dd604 --- /dev/null +++ b/need.md @@ -0,0 +1,39 @@ +请修改当前 RNA 二级结构可视化中的“配对概率”视图，重点不是简单美化，而是修正可视化语义和可读性。 +如图：[[1.png]] +当前问题： +1. 现在用同一种绿色的不同透明度表示概率，用户几乎无法区分不同概率值。 +2. 当前视图容易让人误解“配对概率”到底表示节点概率还是边概率。 +3. 小尺寸视图下，透明度编码效果很差，截图后更不清楚。 +4. 右侧图例是 0~1，但图中节点颜色变化不够清晰。 + +修改目标： +1. 不要再使用“同色 + alpha透明度”作为主编码。 +2. 改成“固定不透明度 + 明确的颜色梯度”来表示数值大小。 +3. 优先使用单调、感知一致的 colormap（例如 viridis / plasma / magma），不要用当前这种浅绿色透明度方案。 +4. 节点颜色必须在 0~1 范围内有清楚区分，低值和高值要一眼能看出。 +5. 右侧 colorbar 要和图中实际颜色完全一致。 +6. 保持当前布局基本不变，不要先大改 UI。 + +语义要求： +1. 如果当前数据实际上是“每个碱基的边际配对倾向 / per-base probability”，那么按钮或图例文字不要再直接写 Pair probabilities，而要改成更准确的名称，例如： + - Per-base pairing probability + - Base-wise pairing score + - Marginal pairing probability +2. 如果当前数据确实是配对矩阵 P(i,j)，那就不要只给节点上色，后续需要支持把概率映射到配对边上。 + +这一步先做的事情： +1. 先保留当前节点着色模式。 +2. 去掉 alpha 映射，改成纯颜色映射。 +3. 调整图例标题，使其准确表达“节点值”而不是“边值”。 +4. 检查节点、文字、主链、配对线在高低概率下是否仍然清晰可读。 + +验收标准： +1. 同一颜色不同透明度的问题被彻底移除。 +2. 在小图（类似当前截图尺寸）下，0.2、0.5、0.8 三档能明显区分。 +3. 用户不会再把这个图误解为“边概率图”。 +4. 代码尽量少改，先保证现有功能可用。 + +请直接修改代码，并说明： +1. 改了哪些绘图参数 +2. 改了哪些 label / title / legend 文案 +3. 如果当前数据语义和 Pair probabilities 不一致，请明确指出并给出更合适命名 diff --git a/nupack/app.py b/nupack/app.py new file mode 100644 index 0000000..35b3132 --- /dev/null +++ b/nupack/app.py @@ -0,0 +1,64 @@ +from fastapi import FastAPI +from pydantic import BaseModel +from typing import List +from nupack import Model, Strand, Tube, SetSpec, tube_analysis + +app = FastAPI() + +class ModelInput(BaseModel): + material: str = "rna" + ensemble: str = "stacking" + celsius: float = 37 + sodium: float = 1.0 + magnesium: float = 0.0 + +class StrandInput(BaseModel): + name: str + sequence: str + concentration: float + unit: str = "uM" + +class TubeInput(BaseModel): + name: str = "tube1" + max_size: int = 2 + +class AnalysisPayload(BaseModel): + model: ModelInput + strands: List[StrandInput] + tube: TubeInput + +def unit_to_molar(value: float, unit: str) -> float: + mapping = {"M": 1, "mM": 1e-3, "uM": 1e-6, "nM": 1e-9, "pM": 1e-12} + return value * mapping[unit] + +@app.get("/health") +def health(): + return {"status": "ok"} + +@app.post("/run") +def run(payload: AnalysisPayload): + model = Model( + material=payload.model.material, + ensemble=payload.model.ensemble, + celsius=payload.model.celsius, + sodium=payload.model.sodium, + magnesium=payload.model.magnesium, + ) + + strands = {} + for s in payload.strands: + obj = Strand(s.sequence, name=s.name) + strands[obj] = unit_to_molar(s.concentration, s.unit) + + tube = Tube( + strands=strands, + complexes=SetSpec(max_size=payload.tube.max_size), + name=payload.tube.name, + ) + + result = tube_analysis(tubes=[tube], model=model) + + return { + "status": "success", + "result": str(result) + } diff --git a/nupack/requirements.txt b/nupack/requirements.txt new file mode 100644 index 0000000..16de0d0 --- /dev/null +++ b/nupack/requirements.txt @@ -0,0 +1,8 @@ +fastapi +uvicorn[standard] +numpy +scipy +pip +matplotlib +pandas +jupyterlab diff --git a/nupack/vendor/nupack-4.0.2.0/LICENSE.txt b/nupack/vendor/nupack-4.0.2.0/LICENSE.txt new file mode 100644 index 0000000..ac3c03d --- /dev/null +++ b/nupack/vendor/nupack-4.0.2.0/LICENSE.txt @@ -0,0 +1,19 @@ +NUPACK Software License Agreement for Non-Commercial Academic Use + +Copyright (c) 2003–2022. California Institute of Technology. All rights reserved. + +Use of the NUPACK Python module and/or source code (“Software”) in source form and/or binary form, with or without modification, is permitted for non-commercial academic purposes only, subject to the conditions and disclaimer stated below. + +Conditions + +1. Redistribution of the Software in source form and/or binary form is not permitted. +2. Web applications that use the Software in source form and/or binary form are not permitted. +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote derivative works without specific prior written permission. + +Disclaimer + +THE SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +Contact + +For any questions about this Software License Agreement please contact info@nupack.org \ No newline at end of file diff --git a/nupack/vendor/nupack-4.0.2.0/README.txt b/nupack/vendor/nupack-4.0.2.0/README.txt new file mode 100644 index 0000000..fdc97d9 --- /dev/null +++ b/nupack/vendor/nupack-4.0.2.0/README.txt @@ -0,0 +1,9 @@ +README file for NUPACK 4.0 +Copyright (c) 2003-2022. California Institute of Technology. All Rights Reserved. + +See LICENSE.txt file + +See NUPACK User Guide for installation instructions: +https://docs.nupack.org + +Technical support: support@nupack.org \ No newline at end of file diff --git a/nupack/vendor/nupack-4.0.2.0/package/LICENSE.txt b/nupack/vendor/nupack-4.0.2.0/package/LICENSE.txt new file mode 100644 index 0000000..ac3c03d --- /dev/null +++ b/nupack/vendor/nupack-4.0.2.0/package/LICENSE.txt @@ -0,0 +1,19 @@ +NUPACK Software License Agreement for Non-Commercial Academic Use + +Copyright (c) 2003–2022. California Institute of Technology. All rights reserved. + +Use of the NUPACK Python module and/or source code (“Software”) in source form and/or binary form, with or without modification, is permitted for non-commercial academic purposes only, subject to the conditions and disclaimer stated below. + +Conditions + +1. Redistribution of the Software in source form and/or binary form is not permitted. +2. Web applications that use the Software in source form and/or binary form are not permitted. +3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote derivative works without specific prior written permission. + +Disclaimer + +THE SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +Contact + +For any questions about this Software License Agreement please contact info@nupack.org \ No newline at end of file diff --git a/nupack/vendor/nupack-4.0.2.0/package/nupack-4.0.2.0-cp312-cp312-linux_x86_64.whl b/nupack/vendor/nupack-4.0.2.0/package/nupack-4.0.2.0-cp312-cp312-linux_x86_64.whl new file mode 100644 index 0000000..91b519b Binary files /dev/null and b/nupack/vendor/nupack-4.0.2.0/package/nupack-4.0.2.0-cp312-cp312-linux_x86_64.whl differ diff --git a/rna/TrnaStructureBeautifier.py b/rna/TrnaStructureBeautifier.py new file mode 100644 index 0000000..979c67c --- /dev/null +++ b/rna/TrnaStructureBeautifier.py @@ -0,0 +1,460 @@ +# _*_ coding:utf-8 _*_ +# +# @Version : 1.1 +# @Project : https://github.com/shueho/BioDataTools +# @Time : 2025/7/23 20:00 +# @Update : 2025/8/13 14:00 +# @Author : Hao Xue +# @E-mail : studid@163.com +# @File : TrnaStructureBeautifier.py +# +# Enhancement of tRNA secondary structure diagrams generated using the ViennaRNA package. +import os +import re +import argparse + + +# old = sys.argv[1] +def parse_arguments(): + parser = argparse.ArgumentParser(description="RNA结构图美化参数配置") + + # === 输入控制 === + parser.add_argument( + "-i", "--input", + type=str, + required=True, + help="SVG文件或者SVG文件存放的文件夹路径" + ) + + # === 基础布局 === + parser.add_argument( + "-s", "--size-weight", + type=float, + default=1.4, + help="图形缩放比例（默认 1.4）" + ) + parser.add_argument( + "-p", "--per-row", + type=int, + default=4, + help="每行图片数量（默认 4）" + ) + parser.add_argument( + "-hg", "--horizontal-gap", + type=int, + default=8, + help="图片水平间隔（默认 8）" + ) + parser.add_argument( + "-vg", "--vertical-gap", + type=int, + default=5, + help="图片垂直间隔（默认 5）" + ) + + # === 碱基连线颜色 === + parser.add_argument( + "-ac", "--adjacent-color", + type=str, + default="blue", + help='相邻碱基连线颜色（支持名称或 HEX，如 "blue" 或 "#00FF00"，默认 "blue"）' + ) + parser.add_argument( + "-pc", "--pair-color", + type=str, + default="red", + help='配对碱基连线颜色（支持名称或 HEX，如 "red" 或 "#FF0000"，默认 "red"）' + ) + + # === 碱基圆圈样式 === + parser.add_argument( + "-bf", "--base-fill", + type=str, + default="white", + help='碱基圆圈填充色（默认 "white"）' + ) + parser.add_argument( + "-bs", "--base-stroke", + type=str, + default="black", + help='碱基圆圈轮廓色（默认 "black"）' + ) + + # === 碱基美化图 === + parser.add_argument( + "-A", "--base-a", + type=str, + default="red", + help='美化A碱基圆圈填充色（默认 "red"）' + ) + parser.add_argument( + "-U", "--base-u", + type=str, + default="blue", + help='美化U/T碱基圆圈填充色（默认 "blue"）' + ) + parser.add_argument( + "-G", "--base-g", + type=str, + default="green", + help='美化G碱基圆圈填充色（默认 "green"）' + ) + parser.add_argument( + "-C", "--base-c", + type=str, + default="yellow", + help='美化C碱基圆圈填充色（默认 "yellow"）' + ) + + # === 反密码子控制 === + # parser.add_argument( + # "-af", "--anticodon-file", + # type=str, + # default=None, + # help="反密码子位点文件路径（文本文件，每行一个位置）" + # ) + # parser.add_argument( + # "-aF", "--anti-fill", + # type=str, + # default="red", + # help='反密码子圆圈填充色（默认 "red"）' + # ) + # parser.add_argument( + # "-aS", "--anti-stroke", + # type=str, + # default="black", + # help='反密码子圆圈轮廓色（默认 "black"）' + # ) + + args = parser.parse_args() + + return args + + +args = parse_arguments() +base_color = {"A": args.base_a, "U": args.base_u, "G": args.base_g, "C": args.base_c, } +# 设置目标文件夹名称 +modified_folder = "modified" +# 创建文件夹（如果已存在则忽略） +os.makedirs(modified_folder, exist_ok=True) +print("已创建文件夹: modified") + + +def read_svg(spath): + with open(spath) as f: + return f.read() + + +def set_new_text(x_y_char, fill_=args.base_fill, stroke_=args.base_stroke): + x = x_y_char[0] + y = x_y_char[1] + t = x_y_char[2].upper().replace("T", "U") + # 核心修改：保持圆圈半径r=5不变，调整dy=0使字母在垂直方向也居中（结合text-anchor=middle实现正中心） + return '\n {}\n'.format( + x, y, fill_, stroke_, x, y, t) +# r=圆圈半径，dy=字体垂直位移，font-size=字体大小 + +def set_new_lines(x_y_char_1, x_y_char_2, stroke_): + x1, y1, _ = x_y_char_1 + x2, y2, _ = x_y_char_2 + xy = (x1, y1, x2, y2) + return '\n'.format(*xy, stroke_) + + +# 保留原有 get_add_line 函数，不修改（保持代码秩序） +def get_add_line(svg_text): + # pattern = r']*?y1="([^"]+)"[^>]*?x2="([^"]+)"[^>]*?y2="([^"]+)"' + pattern = r']*\s*x1="([^"]+)"[^>]*\s*y1="([^"]+)"[^>]*\s*x2="([^"]+)"[^>]*\s*y2="([^"]+)"' + matches = re.findall(pattern, svg_text) + result = [] + for x1_str, y1_str, x2_str, y2_str in matches: + try: + x1 = float(x1_str) + y1 = float(y1_str) + x2 = float(x2_str) + y2 = float(y2_str) + result.append(set_new_lines((x1, y1, ""), (x2, y2, ""), args.pair_color)) + except ValueError: + # 如果坐标无法转换为 float，跳过该元素 + continue + return "".join(result) + + +# 核心修改1：extract_text_info 补充返回坐标边界（用于居中计算），生成内容逻辑完全不变 +def extract_text_info(svg_text): + # 正则表达式匹配 标签中的 x、y 和字符内容 + # pattern = r']*>([^<]+)' + pattern = r']*\s*x="([^"]+)"[^>]*\s*y="([^"]+)"[^>]*>(.*?)<\/text>' + matches = re.findall(pattern, svg_text) + result = [] + for x_str, y_str, char in matches: + try: + x = float(x_str) + y = float(y_str) + result.append((x, y, char)) + except ValueError: + # 如果 x 或 y 无法转换为 float，跳过该元素 + continue + + # 计算坐标边界（新增：用于后续居中） + if not result: + min_x = max_x = min_y = max_y = 0 + else: + xs = [p[0] for p in result] + ys = [p[1] for p in result] + min_x, max_x = min(xs), max(xs) + min_y, max_y = min(ys), max(ys) + + # 原有生成内容逻辑完全不变（保证字母在圆圈内） + content = '\n' + content += get_add_line(svg_text) + for i in range(1, len(result)): + # content += set_new_text(result[i]) + content += set_new_lines(result[i - 1], result[i], args.adjacent_color) + content += set_new_text(result[i - 1]) + + content += set_new_text(result[len(result) - 1]) + '\n' + + # 新增返回坐标边界，用于居中计算 + return content, min_x, max_x, min_y, max_y + + +def save_svg(content, outpath): + with open(outpath, "w") as f: + f.write(content) + + +# 核心修改2：modi_svg 改为超大画布 + 自动居中，保留原有逻辑框架 +def modi_svg(path): + t = read_svg(path) + name_ = os.path.basename(path) + if "-" in name_: + name_ = name_.split("-")[0] + else: + name_ = name_.split(".")[0] + + # 提取内容 + 坐标边界（新增） + mid_, min_x, max_x, min_y, max_y = extract_text_info(t) + scale = args.size_weight + + # 1. 超大画布（2000×2000，可自行调大） + canvas_w = 2000 + canvas_h = 2000 + + # 2. 计算结构中心（用于居中） + if min_x == max_x and min_y == max_y: + cx, cy = 0, 0 + else: + cx = (min_x + max_x) / 2 # 结构水平中心 + cy = (min_y + max_y) / 2 # 结构垂直中心 + + # 3. 计算居中偏移量（保证结构在画布正中间） + translate_x = (canvas_w / 2) / scale - cx # 水平居中 + translate_y = (canvas_h / 2) / scale - cy # 垂直居中 + + # 保留原有 head_ 结构，仅替换画布尺寸和偏移量，同时修改名称位置避免重叠 + head_ = ''' + + + + + {} + + '''.format(canvas_h, canvas_w, canvas_w, canvas_h, canvas_h, canvas_w, scale, scale, translate_x, translate_y, cx, min_y - 20, name_) + + tail_ = ' \n' + save_svg(head_ + mid_ + tail_, "modified/modified_" + name_ + ".svg") + return name_ + + +# save_svg(modi_svg(old), "new.svg") +def remove_xml_declaration(content): + """Remove XML declaration from SVG content.""" + return re.sub(r'<\?xml[^>]*\?>\s*', '', content) + + +def extract_dimensions(svg_content): + """Extract width and height from SVG content.""" + width_match = re.search(r'width="([^"]+)"', svg_content) + height_match = re.search(r'height="([^"]+)"', svg_content) + if not width_match or not height_match: + raise ValueError("Could not extract width or height") + return width_match.group(1), height_match.group(1) + + +def merge_horizontal_group(svg_files, output_file): + """Merge SVGs in a group horizontally.""" + svgs = [] + total_width = 0 + max_height = 0 + gap = args.horizontal_gap + + for filename in svg_files: + with open(filename, 'r') as f: + content = f.read() + content = remove_xml_declaration(content) + width, height = extract_dimensions(content) + try: + width_px = float(width.strip('px')) + height_px = float(height.strip('px')) + except ValueError: + raise ValueError(f"Invalid unit in file {filename}") + + svgs.append({ + 'filename': filename, + 'content': content, + 'width': width, + 'height': height, + 'width_px': width_px, + 'height_px': height_px + }) + + total_width += sum(svg_info['width_px'] for svg_info in svgs) + gap * (len(svgs) - 1) + max_height = max(max_height, height_px) + + new_svg = [ + '', + ''.format(total_width, max_height) + ] + + current_x = 0.0 + for svg_info in svgs: + inner_content = re.sub(r']*>', '', svg_info['content']) + inner_content = re.sub(r'', '', inner_content) + new_svg.append(''.format(current_x)) + new_svg.append(inner_content) + new_svg.append('') + current_x += svg_info['width_px'] + gap + + new_svg.append('') + + with open(output_file, 'w') as f: + f.write('\n'.join(new_svg)) + + +def merge_vertical_groups(group_files, output_file): + """Merge all horizontal groups vertically.""" + svgs = [] + max_width = 0 + total_height = 0 + gap = args.vertical_gap + + for filename in group_files: + with open(filename, 'r') as f: + content = f.read() + content = remove_xml_declaration(content) + width, height = extract_dimensions(content) + try: + width_px = float(width.strip('px')) + height_px = float(height.strip('px')) + except ValueError: + raise ValueError(f"Invalid unit in file {filename}") + + svgs.append({ + 'filename': filename, + 'content': content, + 'width': width, + 'height': height, + 'width_px': width_px, + 'height_px': height_px + }) + + max_width = max(max_width, width_px) + total_height += sum(svg_info['height_px'] for svg_info in svgs) + gap * (len(svgs) - 1) + + new_svg = [ + '', + ''.format(max_width, total_height) + ] + + current_y = 0.0 + for svg_info in svgs: + inner_content = re.sub(r']*>', '', svg_info['content']) + inner_content = re.sub(r'', '', inner_content) + new_svg.append(''.format(current_y)) + new_svg.append(inner_content) + new_svg.append('') + current_y += svg_info['height_px'] + gap + + new_svg.append('') + + with open(output_file, 'w') as f: + f.write('\n'.join(new_svg)) + + +def update_circle_colors(svg_content, base_color_map): + new_ = "" + ls = svg_content.split("\n") + for i in range(len(ls)): + if ")[A-Z](?=<\/text>)', ls[i + 2])[0] + if flag not in base_color_map: + new_ += ls[i] + "\n" + continue + new_ += ls[i].replace(args.base_fill, base_color_map[flag]) + "\n" + else: + new_ += ls[i] + "\n" + return new_ + + +if __name__ == "__main__": + # 示例调用 + print("参数解析结果:") + # print(f"主数据文件: {args.input}") + print(f"大小权重: {args.size_weight}") + print(f"每行图片数量: {args.per_row}") + print(f"图片水平间隔: {args.horizontal_gap}") + print(f"图片垂直间隔: {args.vertical_gap}") + print(f"相邻碱基连线颜色: {args.adjacent_color}") + print(f"配对碱基连线颜色: {args.pair_color}") + print(f"碱基填充色: {args.base_fill}") + print(f"碱基轮廓色: {args.base_stroke}") + # print(f"反密码子文件: {args.anticodon_file}") + # print(f"反密码子填充色: {args.anti_fill}") + # print(f"反密码子轮廓色: {args.anti_stroke}") + # 如果是文件，直接添加到列表 + file_list = [] + group_size = args.per_row + if os.path.isfile(args.input): + file_list.append(args.input) + print(f"已添加文件: {args.input}") + + # 如果是文件夹，添加文件夹内所有文件（不递归子文件夹） + elif os.path.isdir(args.input): + for filename in os.listdir(args.input): + file_path = os.path.join(args.input, filename) + if os.path.isfile(file_path) and "modified_" not in file_path and ".svg" == file_path[-4:]: # 只处理文件 + file_list.append(file_path) + file_list.sort() + print(f"已添加文件夹内所有文件: {args.input}") + modi_files = [] + for i in file_list: + n = modi_svg(i) + modi_files.append("modified/modified_" + n + ".svg") + + # Step 2: Merge into horizontal groups + group_index = 0 + group_files = [] + for i in range(0, len(modi_files), group_size): + group = modi_files[i:i + group_size] + if not group: + break + group_output = f"modified/group_{group_index}.svg" + merge_horizontal_group(group, str(group_output)) + group_files.append(group_output) + group_index += 1 + + # Step 3: Merge all groups vertically + final_output = "modified/final.svg" + merge_vertical_groups(group_files, str(final_output)) + print(f"未上色SVG保存到: {final_output}") + print(f"美化A碱基填充色: {args.base_a}") + print(f"美化U碱基填充色: {args.base_u}") + print(f"美化G碱基填充色: {args.base_g}") + print(f"美化C碱基填充色: {args.base_c}") + fisvg = read_svg(str(final_output)) + cosvg = update_circle_colors(fisvg, base_color) + save_svg(cosvg, "modified/final_color.svg") + # 核心修改3：修复打印路径错误（原写成 final.svg，改为 final_color.svg） + print(f"上色SVG保存到: modified/final_color.svg") +# print(modi_files) \ No newline at end of file diff --git a/rna/ViennaRNA-2.7.2.tar.gz b/rna/ViennaRNA-2.7.2.tar.gz new file mode 100644 index 0000000..b892282 Binary files /dev/null and b/rna/ViennaRNA-2.7.2.tar.gz differ diff --git a/rna/nupack代码示例.txt b/rna/nupack代码示例.txt new file mode 100644 index 0000000..850b1e8 --- /dev/null +++ b/rna/nupack代码示例.txt @@ -0,0 +1,10 @@ +model1 = Model(material='dna', celsius=37,sodium=0.1, magnesium=0.02) +Y1 = Strand('CGTTAACGCAGTGAGGACGGTAGTTTGTCGTTCCATCGCACC', name='Y1') +Y2 = Strand('CGTTAACGGGTGCGATGGAACGACTTTCGACAGGCCTGGTGTAATTTCACCCATGTTAGTCGA', name='Y2') +Y3 = Strand('ATACGGAAAATGGAGATAGGAAGAGTACAATGTCAGCGATAAATTCCGTATACGACACCAGGCCTGTCGATTACTACCGTCCTCACTG', name='Y3') + +t1 = Tube({Y1: 1e-8, Y2: 1e-8, Y3:1e-8 },complexes=SetSpec(max_size=4), name='Tube 1') + +my_result = tube_analysis(tubes=[t1], model=model1, + compute=['pfunc', 'pairs', 'mfe', 'sample', 'subopt']) +print(my_result) diff --git a/rna/工具调用.txt b/rna/工具调用.txt new file mode 100644 index 0000000..52bcb81 --- /dev/null +++ b/rna/工具调用.txt @@ -0,0 +1,139 @@ +from nupack import * +import subprocess +import os + +# ---------- 独立的辅助函数 ---------- +def get_sequence(Complex, result, skip_errors=True): + """ + 获取复合物的序列。 + - 若 Complex 为 None，返回字典 {复合物名: 序列}，跳过解析失败的复合物。 + - 若 Complex 为字符串或 Complex 对象，返回该复合物的序列字符串；失败时抛出异常（除非 skip_errors=False 时由调用者处理）。 + """ + # 统一转换为字符串名称 + if Complex is not None and not isinstance(Complex, str): + Complex = Complex.name if hasattr(Complex, 'name') else str(Complex) + + def _get_seq_from_strands(strand_names, comp_name=None): + """从链名列表获取拼接后的序列，如果失败则抛出异常""" + seq_parts = [] + for s in strand_names: + if not s.isidentifier(): + raise NameError( + f"复合物 '{comp_name}' 中的链名 '{s}' 不是有效的 Python 标识符，" + "请检查链的 name 属性是否包含空格或特殊字符。" + ) + try: + strand_obj = eval(s) + if not isinstance(strand_obj, Strand): + raise TypeError(f"'{s}' 不是 Strand 对象") + seq_parts.append(str(strand_obj)) + except NameError: + raise NameError(f"链名 '{s}' 未定义，请确保该变量已创建为 Strand 对象。") + return ''.join(seq_parts) + + if Complex is None: + # 批量模式 + sequences = {} + for comp in result: + name = comp.name if hasattr(comp, 'name') else str(comp) + try: + # 去除括号，按 '+' 分割 + inp = ''.join([ch for ch in name if ch not in '()']) + strands = inp.split('+') + sequences[name] = _get_seq_from_strands(strands, comp_name=name) + except Exception as e: + if skip_errors: + print(f"⚠️ 跳过复合物 '{name}'，原因：{e}") + else: + raise # 如果不跳过，则抛出异常 + return sequences + else: + # 单复合物模式 + inp = ''.join([ch for ch in Complex if ch not in '()']) + strands = inp.split('+') + return _get_seq_from_strands(strands, comp_name=Complex) + +def get_structure(Complex, result, skip_errors=True): + """ + 获取复合物的 MFE 结构（已去除 '+'）。 + - 若 Complex 为 None，返回字典 {复合物名: 结构}，跳过 result 中不存在的复合物。 + - 若 Complex 为字符串或 Complex 对象，返回该复合物的结构字符串；失败时抛出异常。 + """ + if Complex is not None and not isinstance(Complex, str): + Complex = Complex.name if hasattr(Complex, 'name') else str(Complex) + + if Complex is None: + structures = {} + for comp in result: + name = comp.name if hasattr(comp, 'name') else str(comp) + try: + structure = result[comp].mfe[0].structure + structures[name] = str(structure).replace('+', '') + except Exception as e: + if skip_errors: + print(f"⚠️ 跳过复合物 '{name}'，获取结构失败：{e}") + else: + raise + return structures + else: + structure = result[Complex].mfe[0].structure + return str(structure).replace('+', '') + +# ---------- 绘图函数（带错误跳过） ---------- +def danlian(Complex, result, output_folder="/home/zsy"): + """ + 绘制单个复合物的二级结构图。如果处理过程中出现错误，打印错误信息并返回，不中断程序。 + Complex 可为复合物名称字符串或 Complex 对象。 + """ + try: + sequence = get_sequence(Complex, result, skip_errors=False) # 这里不跳过，让异常抛出以便捕获 + best_structure = get_structure(Complex, result, skip_errors=False) + except Exception as e: + name = Complex if isinstance(Complex, str) else (Complex.name if hasattr(Complex, 'name') else str(Complex)) + print(f"❌ 处理复合物 '{name}' 时出错：{e}，已跳过绘图") + return + + # 获取用于文件名的字符串名称 + if isinstance(Complex, str): + seq_name = Complex + else: + seq_name = Complex.name if hasattr(Complex, 'name') else str(Complex) + + temp_file_path = os.path.join(output_folder, f"{seq_name}.seq") + try: + with open(temp_file_path, "w", newline='\n') as f: + f.write(f">{seq_name}\n") + f.write(f"{sequence}\n") + f.write(f"{best_structure}\n") + except Exception as e: + print(f"❌ 写入临时文件失败：{e}") + return + + try: + input_filename = f"{seq_name}.seq" + process = subprocess.run(["RNAplot", "--o", "svg", input_filename], cwd=output_folder, capture_output=True, text=True) + if process.returncode == 0: + print(f"🎉 绘图成功！请前往 {output_folder} 文件夹查看 {seq_name}_ss.svg") + os.remove(temp_file_path) + else: + print(f"❌ RNAplot 运行失败，报错信息：\n{process.stderr}") + except FileNotFoundError: + print("❌ 找不到 RNAplot 命令，请确认环境变量。") + except Exception as e: + print(f"❌ 绘图过程中出现未知错误：{e}") + + +使用方法 +算全部结果 +for comp in my_result: + danlian(comp, my_result, output_folder="/home/zsy") + +算单个结果 +# 假设您的计算结果保存在 my_result 中 +target = '(Y2+Y3)' # 替换为您的目标复合物名称 +# 输出结构 +print(get_structure(target, my_result)) +# 输出序列 +print(get_sequence(target, my_result)) +# 绘制二级结构图 +danlian(target, my_result, output_folder="/home/zsy") \ No newline at end of file diff --git a/rna/操作流程.txt b/rna/操作流程.txt new file mode 100644 index 0000000..f6731be --- /dev/null +++ b/rna/操作流程.txt @@ -0,0 +1,44 @@ +准备工作： +在 Windows 界面中搜索 控制面板 打开 程序-程序和功能-启动或关闭 Windows 功能-勾选‘适用于 Linux 的 Windows 子系统’与‘虚拟机平台’后点击确定后重启电脑 +在 Microsoft Store 中下载 Ubunut +按 Win+X 键 点击 ‘命令提示符（管理员）’ 输入 wsl --install 后回车。下载完毕后关闭即可。 如果遇到链接失败就是网络问题，链接手机热点后重新输入+回车即可 + +安装：打开 Ubunut +将群内的所有压缩包直接保存在D盘 +nupack 安装 +mkdir nupack-latest +cd nupack-latest +cp /mnt/d/nupack-4.0.2.0.zip ./ +sudo apt install unzip +unzip nupack-4.0.2.0.zip +cd .. +wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh +/bin/bash Miniconda3-latest-Linux-x86_64.sh -b +miniconda3/bin/conda update -n base -c defaults conda +rm Miniconda3-latest-Linux-x86_64.sh +export PATH=$HOME/miniconda3/bin:$PATH +echo 'export PATH=$HOME/miniconda3/bin:$PATH' >> ~/.bashrc +conda install numpy scipy pip matplotlib pandas jupyterlab +pip install -U nupack -f ./nupack-latest/nupack-4.0.2.0/package +jupyter lab + +RNAplot 安装 +mkdir rna +cd rna +cp /mnt/d/ViennaRNA-2.7.2.tar.gz ./ +tar -zxvf ViennaRNA-2.7.2.tar.gz +cd ViennaRNA-2.7.2 +（可能出现报错的解决方法 +sudo apt update +sudo apt install build-essentia +sudo apt install pkg-config） +./configure +make +sudo make install + +图片美化方法 + +将生成后的svg图片保存在一个文件夹内，并将该文件夹与群内的 py 文件保存在一个文件夹。后 Shift+鼠标右键 点击在此处打开 Powershell 窗口后 +输入命令 +python TrnaStructureBeautifier.py -i 文件夹名称 用于指定一个文件夹 +python TrnaStructureBeautifier.py -i 文件名称 用于指定一个文件 \ No newline at end of file diff --git a/rna/调用RNAplot的代码.txt b/rna/调用RNAplot的代码.txt new file mode 100644 index 0000000..73187a5 --- /dev/null +++ b/rna/调用RNAplot的代码.txt @@ -0,0 +1,46 @@ +from nupack import * +import subprocess +import os + +def huatu(Complex,result):#Complex是你想查看的复合物，eg: '(StrandA+StrandB)' 引号和括号都要有; + #result是nupack计算（tube_analysis）的结果，eg: nupack指南里用过my_result、my_results和tube_results等, 或者你自己有其他命名 + def jiegou(Complex,result): + m=str(result[Complex].mfe[0].structure) + mm="" + for i in m: + if i != "+": + mm=mm+i + if i=="": + continue + return mm + def xulie(Complex): + inp='' + for i in Complex: + if i != '(' and i != ')': + inp=inp+i + else: + continue + l=inp.split('+',-1) + xulie='' + for i in l: + xulie=xulie+str(eval(i)) + return xulie + sequence=xulie(Complex) + best_structure=jiegou(Complex,result) + seq_name =Complex + output_folder="/home/liusiye/桌面/nupack-4.0.2.0"#指定图片保存在哪个文件夹，需自己更改 + temp_file_path = os.path.join(output_folder, f"{seq_name}.seq") + with open(temp_file_path, "w", newline='\n') as f: + f.write(f">{seq_name}\n") + f.write(f"{sequence}\n") + f.write(f"{best_structure}\n") + try: + input_filename = f"{seq_name}.seq" + process = subprocess.run(["RNAplot", "--o", "svg", input_filename], cwd=output_folder, capture_output=True, text=True) + if process.returncode == 0: + print(f"🎉 绘图成功！请前往 {output_folder} 文件夹查看 {seq_name}_ss.svg") + os.remove(temp_file_path) + else: + print("❌ RNAplot 运行失败，报错信息：\n", process.stderr) + except FileNotFoundError: + print("❌ 找不到 RNAplot 命令，请确认环境变量。") \ No newline at end of file diff --git a/service/design-guide.html b/service/design-guide.html new file mode 100644 index 0000000..d850a25 --- /dev/null +++ b/service/design-guide.html @@ -0,0 +1,191 @@ + + + + + + NUPACK Design 使用笔记 + + + +
+

NUPACK replica / Design workflow notes

+

Design 使用笔记

+

这份笔记只讲当前网页里已经接入的功能：Design domains、target complexes、target tubes、hard constraints、soft constraints 和 defect weights。

+ + + +
+

1. 快速跑通

+

主界面选择 Design，点击 载入示例，会得到一个可直接运行的双链设计：

+ 
domain:
+  a = N10
+
+strands:
+  A = a
+  B = ~a
+
+target complex:
+  AB_target = A+B
+  structure = (10+)10
+
+target tube:
+  AB_target = 1 uM
+  off-target max size = 2
+

~a 表示 domain a 的反向互补，因此 (10+)10 这个全双链目标结构是合法的。

+
+ +
+

2. Targets 怎么填

+ + + + + + + + + + +
区域含义示例
Design Domains给每段可设计区域命名，并写 IUPAC 约束。a = N10
链输入Design 模式下，序列框可以写 domain composition。A = a, B = ~a
Target Complex指定哪些链组成目标复合物，以及目标二级结构。A+B, (10+)10
Target TubeTube design 才需要。指定 on-target 浓度，off-target 由 max size 自动补全。AB_target = 1 uM
+
+ +
+

3. Hard constraints 怎么用

+

Hard constraint 是“必须满足”的约束，写错会让设计空间为空，或直接报错。

+ + + + + + + + + + + +
类型必填字段当前推荐示例
Match左侧、右侧a 与 a
Complementarity左侧、右侧A 与 B
Similarityscope、reference、上下限scope=a, reference=R10
Windowscope、source sequencesscope=a, 每行一个 10 nt source
Patternpatterns；scope 可空表示全局AAAA, UUUU
+

如果看到 “Constraint scope cannot be empty”，就是当前 constraint 类型需要填写 scope，但该字段为空。新版界面会在提交前拦截这类错误。

+
+ +
+

4. Soft constraints 和 weights

+

Soft constraint 不会替代 ensemble defect，只是给优化目标增加加权惩罚。权重越大，设计器越偏向满足它，但也可能变慢。

+

Defect weights 用来告诉设计器哪些 domain、strand、complex 或 tube 更重要。初学时建议先不填，确认 target 能跑通后再逐步添加。

+
+ +
+

5. 为什么 Design 会慢

+

NUPACK design 是优化问题，不是一次性的结构分析。影响耗时的主要因素：

+
    +
  • off-target max size 越大，需要考虑的非目标复合物越多。
    • +
  • f_stop 越小，停止条件越严格，通常越慢。
    • +
  • trials 大于 1 时会跑多个随机种子，NUPACK 源码里会并行提交多个 trial。
    • +
  • hard constraints 太多或互相矛盾，会反复搜索甚至失败。
    • +
+

当前服务已显式设置 NUPACK 的 config.threads，不再只依赖 OMP_NUM_THREADS。本地 NUPACK 源码里有 NUPACK_CUDA 编译选项，但当前 wheel/镜像不是 CUDA 构建；直接“打开 GPU”不能生效，除非重新编译 NUPACK 的 CUDA 版本并替换镜像里的 wheel。

+
+
+ + diff --git a/service/favicon.svg b/service/favicon.svg new file mode 100644 index 0000000..2efd9fd --- /dev/null +++ b/service/favicon.svg @@ -0,0 +1,4 @@ + + + + diff --git a/service/index.html b/service/index.html new file mode 100644 index 0000000..0b0fb88 --- /dev/null +++ b/service/index.html @@ -0,0 +1,5625 @@ + + + + + + + NUPACK + + + +
+
+
+
+
+ + + + +
+
+

+

+
+
-
+
-
+
-
+
-
+
+
+ +
+
+
+
+

+
+ + +
+
+ +
+

+
+ + + + + +
+
+ +
+

+
+
+ +
+
+ +
+

+
+ + + + + + +
+
+ +
+

+
+ + + + + + +
+

+
+ +
+

+
+
+ + + +
+
+
+ +
+
+
+ +
+

+
+ + +
+
+ + + + + +
+ + +
+ +
+

+
+ + + + +
+
+
+
+
+ +
+
+
+
+
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/service/server.py b/service/server.py new file mode 100644 index 0000000..8f0d065 --- /dev/null +++ b/service/server.py @@ -0,0 +1,2092 @@ +import json +import mimetypes +import os +import re +import threading +import time +import traceback +from decimal import Decimal +from concurrent.futures import FIRST_COMPLETED, ThreadPoolExecutor, wait +from http import HTTPStatus +from http.server import BaseHTTPRequestHandler, ThreadingHTTPServer +from pathlib import Path +import subprocess +import tempfile +from urllib.parse import parse_qs, urlparse +from uuid import uuid4 + +import numpy +from nupack import ( + Complex, + Complementarity, + Domain, + Diversity, + EnergyMatch, + DesignOptions, + Library, + Match, + Model, + Pattern, + Similarity, + SSM, + SetSpec, + Strand, + TargetComplex, + TargetStrand, + TargetTube, + Tube, + Weights, + Window, + complex_analysis, + complex_design, + tube_analysis, + tube_design, +) +from nupack import config as nupack_config +from split_strand_svg import render_split_strands_svg + +try: + import redis +except ImportError: + redis = None + + +ROOT = Path(__file__).resolve().parent +INDEX_PATH = ROOT / "index.html" +FAVICON_PATH = ROOT / "favicon.svg" +GUIDE_PATH = ROOT / "design-guide.html" + +HOST = os.environ.get("NP_HOST", "0.0.0.0") +PORT = int(os.environ.get("NP_PORT", "18765")) +RNA_PLOT_CMD = os.environ.get("RNA_PLOT_CMD", "RNAplot") +ENABLE_RNAPLOT = os.environ.get("ENABLE_RNAPLOT", "1") != "0" +STRUCTURE_PLOT_MODE = os.environ.get("NP_STRUCTURE_PLOT_MODE", "auto").strip().lower() +RUN_MODE = os.environ.get("NP_RUN_MODE", "server") +REDIS_URL = os.environ.get("NP_REDIS_URL", "").strip() +JOB_QUEUE_KEY = os.environ.get("NP_JOB_QUEUE_KEY", "np_replica:jobs") +JOB_RUNNING_KEY = os.environ.get("NP_JOB_RUNNING_KEY", "np_replica:jobs:running") +WORKER_CONCURRENCY = int(os.environ.get("NP_WORKER_CONCURRENCY", "2")) +PER_JOB_THREAD_LIMIT = int(os.environ.get("NP_PER_JOB_THREAD_LIMIT", "1")) +NUPACK_CACHE_GB = float(os.environ.get("NP_NUPACK_CACHE_GB", "2.0")) + +UNIT_SCALE = { + "M": 1.0, + "mM": 1e-3, + "uM": 1e-6, + "nM": 1e-9, + "pM": 1e-12, +} +IUPAC_CODES = "ACGTUWSMKRYBDHVN" +IUPAC_CONSTRAINT_TOKEN = re.compile(rf"[{IUPAC_CODES}](?:\d+)?") + +VALID_COMPUTE = {"pfunc", "pairs", "mfe", "sample", "subopt", "ensemble_size"} +JOB_STORE = {} +JOB_LOCK = threading.Lock() +JOB_TTL_SECONDS = int(os.environ.get("NP_JOB_TTL_SECONDS", "3600")) +JOB_MAX_COUNT = int(os.environ.get("NP_JOB_MAX_COUNT", "64")) +SHARE_STORE = {} +SHARE_LOCK = threading.Lock() +SHARE_MAX_COUNT = int(os.environ.get("NP_SHARE_MAX_COUNT", "100")) +SHARE_MAX_BYTES = int(os.environ.get("NP_SHARE_MAX_BYTES", str(12 * 1024 * 1024))) +SHARE_KEY_PREFIX = os.environ.get("NP_SHARE_KEY_PREFIX", "np_replica:shares") +SHARE_INDEX_KEY = f"{SHARE_KEY_PREFIX}:index" +REDIS_CLIENT = None + + +def json_bytes(payload, status=HTTPStatus.OK): + return ( + status, + "application/json; charset=utf-8", + json.dumps(payload, ensure_ascii=False, indent=2).encode("utf-8"), + ) + + +def log_event(message): + print(f"[np-replica] {time.strftime('%Y-%m-%d %H:%M:%S')} {message}", flush=True) + + +def apply_thread_limits(): + # Keep native math libraries and NUPACK's own executor from using all CPU cores inside one job. + thread_limit_int = max(1, PER_JOB_THREAD_LIMIT) + thread_limit = str(thread_limit_int) + for key in ( + "OMP_NUM_THREADS", + "OPENBLAS_NUM_THREADS", + "MKL_NUM_THREADS", + "NUMEXPR_NUM_THREADS", + "VECLIB_MAXIMUM_THREADS", + "GOTO_NUM_THREADS", + ): + os.environ[key] = thread_limit + nupack_config.threads = thread_limit_int + nupack_config.cache = max(0.1, NUPACK_CACHE_GB) + + +def html_bytes(path): + body = path.read_bytes() + mime, _ = mimetypes.guess_type(str(path)) + return HTTPStatus.OK, mime or "text/html; charset=utf-8", body + + +def file_etag(path): + stat = path.stat() + return f'W/"{stat.st_mtime_ns:x}-{stat.st_size:x}"' + + +def redis_enabled(): + return bool(REDIS_URL and redis is not None) + + +def redis_client(): + global REDIS_CLIENT + if not redis_enabled(): + return None + if REDIS_CLIENT is None: + REDIS_CLIENT = redis.Redis.from_url(REDIS_URL, decode_responses=True) + return REDIS_CLIENT + + +def queue_size(): + if not redis_enabled(): + with JOB_LOCK: + return sum(1 for job in JOB_STORE.values() if job.get("status") in {"queued", "running"}) + try: + return int(redis_client().llen(JOB_QUEUE_KEY)) + except Exception: + return None + + +def prune_redis_running_jobs(): + if not redis_enabled(): + return 0 + client = redis_client() + removed = 0 + try: + for job_id in client.smembers(JOB_RUNNING_KEY): + raw = client.get(job_key(job_id)) + if raw is None: + client.srem(JOB_RUNNING_KEY, job_id) + removed += 1 + continue + try: + if json.loads(raw).get("status") != "running": + client.srem(JOB_RUNNING_KEY, job_id) + removed += 1 + except Exception: + client.srem(JOB_RUNNING_KEY, job_id) + removed += 1 + except Exception: + return removed + return removed + + +def job_stats(): + stats = { + "queued": 0, + "running": 0, + "success": 0, + "error": 0, + } + if redis_enabled(): + prune_redis_running_jobs() + client = redis_client() + stats["queued"] = int(client.llen(JOB_QUEUE_KEY)) + stats["running"] = int(client.scard(JOB_RUNNING_KEY)) + return stats + + with JOB_LOCK: + prune_jobs() + for job in JOB_STORE.values(): + status = job.get("status") + if status in stats: + stats[status] += 1 + return stats + + +def job_key(job_id): + return f"np_replica:job:{job_id}" + + +def share_key(share_id): + return f"{SHARE_KEY_PREFIX}:item:{share_id}" + + +def normalize_sequence(sequence): + return "".join(sequence.upper().split()) + + +def normalize_design_sequence(sequence): + return re.sub(r"\s+", "", str(sequence or "").upper()) + + +def is_valid_iupac_constraint(sequence): + seq = normalize_design_sequence(sequence) + if not seq: + return False + index = 0 + while index < len(seq): + match = IUPAC_CONSTRAINT_TOKEN.match(seq, index) + if not match: + return False + index = match.end() + return True + + +def ensure_unit_interval_limits(lower, upper, label): + if not (0 <= lower <= upper <= 1): + raise ValueError(f"{label} limits must satisfy 0 <= lower <= upper <= 1.") + + +def unit_to_molar(value, unit): + if unit not in UNIT_SCALE: + raise ValueError(f"Unsupported concentration unit: {unit}") + return float(value) * UNIT_SCALE[unit] + + +def build_model(model_input): + return Model( + material=model_input.get("material", "rna"), + ensemble=model_input.get("ensemble", "stacking"), + celsius=float(model_input.get("celsius", 37.0)), + sodium=float(model_input.get("sodium", 1.0)), + magnesium=float(model_input.get("magnesium", 0.0)), + ) + + +def build_model_summary(model_input): + return { + "material": model_input.get("material", "rna"), + "ensemble": model_input.get("ensemble", "stacking"), + "celsius": float(model_input.get("celsius", 37.0)), + "sodium": float(model_input.get("sodium", 1.0)), + "magnesium": float(model_input.get("magnesium", 0.0)), + } + + +def parse_compute(payload): + compute = payload.get("compute") or ["pfunc", "mfe"] + compute = [item for item in compute if item in VALID_COMPUTE] + if not compute: + raise ValueError("At least one compute option is required.") + return compute + + +def parse_options(payload): + raw = payload.get("options") or {} + options = { + "num_sample": int(raw.get("num_sample", 20)), + "energy_gap": float(raw.get("energy_gap", 1.0)), + "sparsity_fraction": float(raw.get("sparsity_fraction", 1.0)), + "sparsity_threshold": float(raw.get("sparsity_threshold", 0.0)), + "single_mfe": bool(raw.get("single_mfe", False)), + "result_limit": int(raw.get("result_limit", 25)), + "pairs_preview_size": int(raw.get("pairs_preview_size", 24)), + } + if options["num_sample"] < 0 or options["num_sample"] > 1000: + raise ValueError("num_sample must be between 0 and 1000.") + if options["energy_gap"] < 0: + raise ValueError("energy_gap must be non-negative.") + if not 0 <= options["sparsity_fraction"] <= 1: + raise ValueError("sparsity_fraction must be between 0 and 1.") + if options["sparsity_threshold"] < 0: + raise ValueError("sparsity_threshold must be non-negative.") + if options["result_limit"] < 1: + raise ValueError("result_limit must be at least 1.") + if options["pairs_preview_size"] < 8 or options["pairs_preview_size"] > 128: + raise ValueError("pairs_preview_size must be between 8 and 128.") + return options + + +def parse_design_options(payload): + raw = payload.get("design") or {} + options = { + "trials": int(raw.get("trials", 1)), + "result_limit": int(raw.get("result_limit", 25)), + "off_target_max_size": int(raw.get("off_target_max_size", payload.get("tube", {}).get("max_size", 2))), + "stop_condition": float(raw.get("stop_condition", 0.02)), + "seed": int(raw.get("seed", 0)), + "wobble_mutations": bool(raw.get("wobble_mutations", False)), + "max_time_seconds": int(raw.get("max_time_seconds", 0)), + } + if options["trials"] < 1 or options["trials"] > 8: + raise ValueError("design trials must be between 1 and 8.") + if options["result_limit"] < 1: + raise ValueError("design result_limit must be at least 1.") + if options["off_target_max_size"] < 1 or options["off_target_max_size"] > 8: + raise ValueError("design off_target_max_size must be between 1 and 8.") + if not 0 < options["stop_condition"] < 1: + raise ValueError("design stop_condition must be between 0 and 1.") + if options["max_time_seconds"] < 0: + raise ValueError("design max_time_seconds must be non-negative.") + return options + + +def build_strands(strand_payload): + if not strand_payload: + raise ValueError("At least one strand is required.") + + strand_map = {} + ordered = [] + + for row in strand_payload: + name = (row.get("name") or "").strip() + sequence = normalize_sequence(row.get("sequence") or "") + if not name: + raise ValueError("Every strand needs a name.") + if not sequence: + raise ValueError(f"Strand {name} is missing a sequence.") + if name in strand_map: + raise ValueError(f"Duplicate strand name: {name}") + if not re.fullmatch(r"[ACGTUWSMKRYBDHVN]+", sequence): + raise ValueError(f"Strand {name} contains unsupported characters.") + concentration = float(row.get("concentration", 0)) + if concentration < 0: + raise ValueError(f"Strand {name} concentration must be non-negative.") + + strand = Strand(sequence, name=name) + strand_map[name] = strand + ordered.append( + { + "name": name, + "sequence": sequence, + "concentration": concentration, + "unit": row.get("unit", "uM"), + "object": strand, + } + ) + + return strand_map, ordered + + +def build_design_domains(domain_payload): + domain_map = {} + ordered = [] + for row in domain_payload: + name = (row.get("name") or "").strip() + sequence_constraint = normalize_design_sequence(row.get("sequence") or "") + if not name: + raise ValueError("Every design domain needs a name.") + if not sequence_constraint: + raise ValueError(f"Design domain {name} is missing a sequence constraint.") + if name in domain_map: + raise ValueError(f"Duplicate domain name: {name}") + if not is_valid_iupac_constraint(sequence_constraint): + raise ValueError(f"Design domain {name} contains unsupported constraint characters.") + + domain = Domain(sequence_constraint, name=name) + domain_map[name] = domain + ordered.append( + { + "name": name, + "constraint": sequence_constraint, + "object": domain, + } + ) + return domain_map, ordered + + +def parse_domain_composition(text, domain_map): + tokens = [token.strip() for token in re.split(r"[\s,]+", str(text or "").strip()) if token.strip()] + if not tokens: + raise ValueError("Domain composition cannot be empty.") + + domains = [] + for token in tokens: + complement = False + domain_name = token + if token.startswith("~"): + complement = True + domain_name = token[1:] + elif token.endswith("*"): + complement = True + domain_name = token[:-1] + + if not domain_name: + raise ValueError(f"Invalid domain token: {token}") + if domain_name not in domain_map: + raise ValueError(f"Unknown domain name in strand composition: {domain_name}") + + domain = domain_map[domain_name] + domains.append(~domain if complement else domain) + return domains + + +def build_design_strands(strand_payload, domain_map=None): + if not strand_payload: + raise ValueError("At least one design strand is required.") + + domain_map = domain_map or {} + use_domain_composition = bool(domain_map) + target_strand_map = {} + ordered = [] + + for row in strand_payload: + name = (row.get("name") or "").strip() + raw_definition = str(row.get("sequence") or "").strip() + sequence_constraint = normalize_design_sequence(raw_definition) + if not name: + raise ValueError("Every design strand needs a name.") + if not raw_definition: + raise ValueError(f"Design strand {name} is missing a strand definition.") + if name in target_strand_map: + raise ValueError(f"Duplicate strand name: {name}") + + if use_domain_composition: + try: + strand_domains = parse_domain_composition(raw_definition, domain_map) + except ValueError: + if (" " in raw_definition) or ("," in raw_definition) or ("~" in raw_definition) or ("*" in raw_definition): + raise + if not is_valid_iupac_constraint(sequence_constraint): + raise + inline_domain = Domain(sequence_constraint, name=name) + strand_domains = [inline_domain] + target_strand = TargetStrand(strand_domains, name=name) + constraint_kind = "sequence_constraint" + constraint_value = sequence_constraint + else: + inline_domain = None + target_strand = TargetStrand(strand_domains, name=name) + constraint_kind = "domain_composition" + constraint_value = raw_definition + else: + if not is_valid_iupac_constraint(sequence_constraint): + raise ValueError(f"Design strand {name} contains unsupported constraint characters.") + inline_domain = Domain(sequence_constraint, name=name) + strand_domains = [inline_domain] + target_strand = TargetStrand(strand_domains, name=name) + constraint_kind = "sequence_constraint" + constraint_value = sequence_constraint + + target_strand_map[name] = target_strand + ordered.append( + { + "name": name, + "definition": raw_definition, + "constraint": constraint_value, + "constraint_kind": constraint_kind, + "object": target_strand, + "domains": strand_domains, + "domain": inline_domain, + } + ) + + return target_strand_map, ordered + + +def resolve_scope_reference(expr, domain_map, strand_map, *, allow_global=False): + text = str(expr or "").strip() + if not text: + if allow_global: + return None + raise ValueError("Constraint scope cannot be empty.") + + if text in strand_map: + return strand_map[text] + return parse_domain_composition(text, domain_map) + + +def resolve_domain_scope(expr, domain_map): + text = str(expr or "").strip() + if not text: + raise ValueError("Domain scope cannot be empty.") + return parse_domain_composition(text, domain_map) + + +def resolve_target_complex_scope(expr, target_complex_map): + text = str(expr or "").strip() + if not text: + raise ValueError("Target complex scope cannot be empty.") + tokens = [token.strip() for token in re.split(r"[\s,]+", text) if token.strip()] + if not tokens: + raise ValueError("Target complex scope cannot be empty.") + missing = [token for token in tokens if token not in target_complex_map] + if missing: + raise ValueError(f"Unknown target complex name(s): {', '.join(missing)}") + return [target_complex_map[token] for token in tokens] + + +def parse_sequence_list(text): + if isinstance(text, list): + values = text + else: + normalized = str(text or "").replace("|", "\n").replace(";", "\n") + values = normalized.splitlines() + output = [item.strip() for item in values if item and item.strip()] + if not output: + raise ValueError("Constraint sequence list cannot be empty.") + for seq in output: + if not re.fullmatch(rf"[{IUPAC_CODES}]+", seq.upper()): + raise ValueError(f"Invalid sequence source: {seq!r}") + return output + + +def parse_pattern_list(text): + if isinstance(text, list): + values = text + else: + values = re.split(r"[\n,]+", str(text or "")) + output = [item.strip().upper() for item in values if item and item.strip()] + if not output: + raise ValueError("Pattern constraint requires at least one pattern.") + for pattern in output: + if not is_valid_iupac_constraint(pattern): + raise ValueError(f"Invalid pattern: {pattern!r}") + return output + + +def parse_catalog_list(text): + raw_lines = [line.strip() for line in str(text or "").splitlines() if line.strip()] + if not raw_lines: + raise ValueError("Library constraint requires at least one catalog row.") + catalog = [] + for line in raw_lines: + library = [item.strip().upper() for item in re.split(r"[,|]+", line) if item.strip()] + if not library: + raise ValueError("Library constraint contains an empty catalog row.") + for seq in library: + if not re.fullmatch(rf"[{IUPAC_CODES}]+", seq): + raise ValueError(f"Invalid library sequence: {seq!r}") + catalog.append(library) + return catalog + + +def build_hard_constraints(payload, domain_map, strand_map): + constraints = [] + for index, row in enumerate(payload.get("hard_constraints") or [], start=1): + constraint_type = (row.get("type") or "").strip().lower() + if not constraint_type: + continue + + try: + def require_text(field, label): + value = str(row.get(field) or "").strip() + if not value: + raise ValueError(f"{label} is required.") + return value + + if constraint_type == "match": + constraints.append( + Match( + resolve_scope_reference(require_text("left", "Match left scope"), domain_map, strand_map), + resolve_scope_reference(require_text("right", "Match right scope"), domain_map, strand_map), + ) + ) + elif constraint_type == "complementarity": + constraints.append( + Complementarity( + resolve_scope_reference(require_text("left", "Complementarity left scope"), domain_map, strand_map), + resolve_scope_reference(require_text("right", "Complementarity right scope"), domain_map, strand_map), + wobble_mutations=bool(row.get("wobble_mutations", False)), + ) + ) + elif constraint_type == "diversity": + word = int(row.get("word", 4)) + types = int(row.get("types", 2)) + if word < 1: + raise ValueError("Diversity word must be >= 1.") + if types < 1 or types > 4: + raise ValueError("Diversity types must be between 1 and 4.") + kwargs = { + "word": word, + "types": types, + } + scope = resolve_scope_reference(row.get("scope"), domain_map, strand_map, allow_global=True) + if scope is not None: + kwargs["scope"] = scope + constraints.append(Diversity(**kwargs)) + elif constraint_type == "similarity": + lower = float(row.get("min_fraction", 0.0)) + upper = float(row.get("max_fraction", 1.0)) + ensure_unit_interval_limits(lower, upper, "Similarity") + reference = normalize_design_sequence(row.get("reference") or "") + if not is_valid_iupac_constraint(reference): + raise ValueError("Similarity reference must be a valid IUPAC constraint.") + constraints.append( + Similarity( + resolve_scope_reference(require_text("scope", "Similarity scope"), domain_map, strand_map), + reference, + limits=[lower, upper], + ) + ) + elif constraint_type == "window": + constraints.append( + Window( + resolve_scope_reference(require_text("scope", "Window scope"), domain_map, strand_map), + sources=parse_sequence_list(row.get("sources") or ""), + ) + ) + elif constraint_type == "library": + constraints.append( + Library( + resolve_scope_reference(require_text("scope", "Library scope"), domain_map, strand_map), + catalog=parse_catalog_list(row.get("catalog") or ""), + ) + ) + elif constraint_type == "pattern": + kwargs = {"patterns": parse_pattern_list(row.get("patterns") or "")} + scope = resolve_scope_reference(row.get("scope"), domain_map, strand_map, allow_global=True) + if scope is not None: + kwargs["scope"] = scope + constraints.append(Pattern(**kwargs)) + else: + raise ValueError(f"Unsupported hard constraint type: {constraint_type}") + except Exception as exc: + raise ValueError(f"Invalid hard constraint #{index}: {exc}") from exc + + return constraints + + +def build_soft_constraints(payload, domain_map, strand_map, target_complex_map): + constraints = [] + for index, row in enumerate(payload.get("soft_constraints") or [], start=1): + constraint_type = (row.get("type") or "").strip().lower() + if not constraint_type: + continue + + try: + def require_text(field, label): + value = str(row.get(field) or "").strip() + if not value: + raise ValueError(f"{label} is required.") + return value + + if constraint_type == "pattern": + weight = float(row.get("weight", 1.0)) + if weight < 0: + raise ValueError("Pattern weight must be non-negative.") + kwargs = { + "patterns": parse_pattern_list(row.get("patterns") or ""), + "weight": weight, + } + scope = resolve_scope_reference(row.get("scope"), domain_map, strand_map, allow_global=True) + if scope is not None: + kwargs["scope"] = scope + constraints.append(Pattern(**kwargs)) + elif constraint_type == "similarity": + lower = float(row.get("min_fraction", 0.0)) + upper = float(row.get("max_fraction", 1.0)) + ensure_unit_interval_limits(lower, upper, "Similarity") + reference = normalize_design_sequence(row.get("reference") or "") + if not is_valid_iupac_constraint(reference): + raise ValueError("Similarity reference must be a valid IUPAC constraint.") + weight = float(row.get("weight", 1.0)) + if weight < 0: + raise ValueError("Similarity weight must be non-negative.") + constraints.append( + Similarity( + resolve_scope_reference(require_text("scope", "Similarity scope"), domain_map, strand_map), + reference, + limits=[lower, upper], + weight=weight, + ) + ) + elif constraint_type == "ssm": + word = int(row.get("word", 4)) + if word < 1: + raise ValueError("SSM word must be >= 1.") + weight = float(row.get("weight", 1.0)) + if weight < 0: + raise ValueError("SSM weight must be non-negative.") + kwargs = { + "word": word, + "weight": weight, + } + scope = resolve_target_complex_scope(row.get("scope"), target_complex_map) if row.get("scope") else None + if scope is not None: + kwargs["scope"] = scope + constraints.append(SSM(**kwargs)) + elif constraint_type == "energy_match": + weight = float(row.get("weight", 1.0)) + if weight < 0: + raise ValueError("EnergyMatch weight must be non-negative.") + kwargs = { + "domains": resolve_domain_scope(require_text("scope", "EnergyMatch domain scope"), domain_map), + "weight": weight, + } + if row.get("energy_ref") not in {None, ""}: + kwargs["energy_ref"] = float(row.get("energy_ref")) + constraints.append(EnergyMatch(**kwargs)) + else: + raise ValueError(f"Unsupported soft constraint type: {constraint_type}") + except Exception as exc: + raise ValueError(f"Invalid soft constraint #{index}: {exc}") from exc + + return constraints + + +def parse_design_complexes(payload, target_strand_map): + target_rows = payload.get("design_complexes") or payload.get("design_targets") or [] + if not target_rows: + raise ValueError("At least one design target complex is required.") + + targets = [] + target_complex_map = {} + + for idx, row in enumerate(target_rows, start=1): + name = (row.get("name") or f"target_{idx}").strip() or f"target_{idx}" + if name in target_complex_map: + raise ValueError(f"Duplicate design target complex name: {name}") + strand_text = (row.get("strands") or "").strip() + structure = (row.get("structure") or "").strip() + if not strand_text: + raise ValueError(f"Design target complex {name} is missing strand membership.") + if not structure: + raise ValueError(f"Design target complex {name} is missing a target structure.") + + tokens = [token.strip() for token in strand_text.split("+") if token.strip()] + if not tokens: + raise ValueError(f"Design target complex {name} has an invalid strand list.") + missing = [token for token in tokens if token not in target_strand_map] + if missing: + raise ValueError( + f"Unknown strand name(s) in design target complex {name}: {', '.join(missing)}" + ) + + target_complex = TargetComplex( + [target_strand_map[token] for token in tokens], + structure, + name=name, + ) + target_payload = { + "name": name, + "strands": tokens, + "structure": structure, + "object": target_complex, + } + targets.append(target_payload) + target_complex_map[name] = target_complex + + return targets, target_complex_map + + +def parse_design_tubes(payload, target_rows, target_complex_map, default_max_size): + tube_rows = payload.get("design_tubes") or [] + if not tube_rows: + legacy_targets = payload.get("design_targets") or [] + if legacy_targets: + on_targets = [] + for row in legacy_targets: + if row.get("concentration") in {None, ""}: + continue + on_targets.append( + { + "complex": row.get("name"), + "concentration": row.get("concentration"), + "unit": row.get("unit", "uM"), + } + ) + if on_targets: + tube_rows = [{ + "name": ((payload.get("tube") or {}).get("name") or "design_tube"), + "max_size": default_max_size, + "on_targets": on_targets, + }] + + if not tube_rows: + raise ValueError("Tube design requires at least one target tube.") + + tubes = [] + ordered_rows = [] + tube_name_set = set() + for index, row in enumerate(tube_rows, start=1): + tube_name = (row.get("name") or f"tube_{index}").strip() or f"tube_{index}" + if tube_name in tube_name_set: + raise ValueError(f"Duplicate target tube name: {tube_name}") + tube_name_set.add(tube_name) + max_size = int(row.get("max_size", default_max_size)) + if max_size < 1 or max_size > 8: + raise ValueError(f"Target tube {tube_name} max_size must be between 1 and 8.") + + on_targets_input = row.get("on_targets") or [] + if not on_targets_input: + raise ValueError(f"Target tube {tube_name} must include at least one on-target complex.") + + on_targets = {} + serialized_on_targets = [] + seen_complex_names = set() + for entry in on_targets_input: + complex_name = (entry.get("complex") or "").strip() + if complex_name not in target_complex_map: + raise ValueError(f"Unknown target complex {complex_name} in tube {tube_name}.") + if complex_name in seen_complex_names: + raise ValueError(f"Duplicate on-target complex {complex_name} in tube {tube_name}.") + seen_complex_names.add(complex_name) + concentration = float(entry.get("concentration", 0)) + unit = entry.get("unit", "uM") + if concentration <= 0: + raise ValueError(f"Target tube {tube_name} requires positive target concentrations.") + concentration_M = unit_to_molar(concentration, unit) + complex_obj = target_complex_map[complex_name] + on_targets[complex_obj] = concentration_M + serialized_on_targets.append( + { + "complex": complex_name, + "concentration": concentration, + "unit": unit, + "target_concentration_M": concentration_M, + } + ) + + tube = TargetTube( + on_targets=on_targets, + off_targets=SetSpec(max_size=max_size), + name=tube_name, + ) + ordered_rows.append( + { + "name": tube_name, + "max_size": max_size, + "on_targets": serialized_on_targets, + "object": tube, + } + ) + tubes.append(tube) + + return ordered_rows, tubes + + +def validate_design_object_names(design_domains, design_strands, target_rows, tube_rows): + name_map = {} + for kind, rows in ( + ("domain", design_domains), + ("strand", design_strands), + ("target complex", target_rows), + ("target tube", tube_rows), + ): + for row in rows: + name = (row.get("name") or "").strip() + if not name: + continue + prev = name_map.get(name) + if prev and prev != kind: + raise ValueError( + f"Object name {name!r} is reused across {prev} and {kind}. " + "Design object names must be unique in a job." + ) + name_map[name] = kind + + +def build_defect_weights(payload, design_mode, design_domains, design_strands, target_rows, tube_rows): + rows = payload.get("defect_weights") or [] + if not rows: + return None + + domain_index = {item["name"]: item["object"] for item in design_domains} + strand_index = {item["name"]: item["object"] for item in design_strands} + complex_index = {item["name"]: item["object"] for item in target_rows} + tube_index = {item["name"]: item["object"] for item in tube_rows} + weight_base = [item["object"] for item in tube_rows] if design_mode == "tube" else [item["object"] for item in target_rows] + weights = Weights(weight_base) + + for index, row in enumerate(rows, start=1): + kind = (row.get("kind") or "").strip().lower() + value = float(row.get("weight", 1.0)) + if value < 0: + raise ValueError(f"Invalid defect weight #{index}: weight must be non-negative") + name = (row.get("name") or "").strip() + try: + if kind == "global": + weights[:] = value + elif kind == "domain": + weights[domain_index[name]] = value + elif kind == "strand": + weights[:, strand_index[name]] = value + elif kind == "complex": + weights[:, :, complex_index[name]] = value + elif kind == "tube": + if design_mode != "tube": + raise ValueError("tube-level defect weights are only valid in tube design mode") + weights[:, :, :, tube_index[name]] = value + else: + raise ValueError(f"Unsupported defect weight kind: {kind}") + except KeyError as exc: + raise ValueError(f"Unknown entity for defect weight #{index}: {name}") from exc + except Exception as exc: + raise ValueError(f"Invalid defect weight #{index}: {exc}") from exc + + return weights + + +def parse_complex_lines(text, strand_map): + complexes = [] + raw_lines = [line.strip() for line in text.splitlines() if line.strip()] + if not raw_lines: + raise ValueError("Complex mode requires at least one complex definition.") + + for idx, line in enumerate(raw_lines, start=1): + tokens = [token.strip() for token in line.split("+") if token.strip()] + if not tokens: + raise ValueError(f"Invalid complex definition on line {idx}: {line}") + missing = [token for token in tokens if token not in strand_map] + if missing: + raise ValueError( + f"Unknown strand name(s) in complex line {idx}: {', '.join(missing)}" + ) + strands = [strand_map[token] for token in tokens] + complexes.append(Complex(strands, name=f"complex_{idx}")) + + return complexes + + +def stringify_complex(complex_obj): + return " + ".join(strand.name for strand in complex_obj.strands) + + +def flatten_sequence(complex_obj): + return "".join(str(strand) for strand in complex_obj.strands) + + +def sanitize_structure_for_rnaplot(structure): + return str(structure).replace("+", "") + + +def safe_name(name): + return re.sub(r"[^A-Za-z0-9_.-]+", "_", name).strip("_") or "structure" + + +def render_rnaplot_svg(complex_obj, structure, suffix): + if not ENABLE_RNAPLOT: + return None + + seq_name = safe_name(f"{stringify_complex(complex_obj)}_{suffix}") + sequence = flatten_sequence(complex_obj) + structure_text = sanitize_structure_for_rnaplot(structure) + + with tempfile.TemporaryDirectory(prefix="rnaplot-") as tmpdir: + tmp_path = Path(tmpdir) + input_path = tmp_path / f"{seq_name}.seq" + output_path = tmp_path / f"{seq_name}_ss.svg" + + input_path.write_text( + f">{seq_name}\n{sequence}\n{structure_text}\n", + encoding="utf-8", + newline="\n", + ) + + process = subprocess.run( + [RNA_PLOT_CMD, "-f", "svg", input_path.name], + cwd=tmpdir, + capture_output=True, + text=True, + timeout=20, + ) + + if process.returncode != 0: + raise RuntimeError(process.stderr.strip() or "RNAplot failed.") + + if not output_path.exists(): + raise RuntimeError("RNAplot finished without producing an SVG file.") + + return output_path.read_text(encoding="utf-8") + + +def render_structure_svg(complex_obj, structure, suffix): + structure_text = str(structure) + multistrand = len(complex_obj.strands) > 1 and "+" in structure_text + should_try_split = STRUCTURE_PLOT_MODE == "split" or ( + STRUCTURE_PLOT_MODE != "rnaplot" and multistrand + ) + split_error = None + + if should_try_split: + try: + return render_split_strands_svg( + [str(strand) for strand in complex_obj.strands], + structure_text, + title=safe_name(f"{stringify_complex(complex_obj)}_{suffix}"), + ) + except Exception as exc: + split_error = exc + if STRUCTURE_PLOT_MODE == "split": + raise RuntimeError(f"Split-strand layout failed: {exc}") from exc + + if STRUCTURE_PLOT_MODE == "split": + if split_error is not None: + raise RuntimeError(f"Split-strand layout failed: {split_error}") from split_error + return None + + try: + return render_rnaplot_svg(complex_obj, structure, suffix) + except Exception as exc: + if split_error is not None: + raise RuntimeError( + f"Split-strand layout failed: {split_error}; RNAplot fallback failed: {exc}" + ) from exc + raise + + +def serialize_structures(complex_obj, items, plot_kind=None): + output = [] + for index, item in enumerate(items or []): + row = { + "structure": str(item.structure), + "energy": round(float(item.energy), 6), + "stack_energy": round(float(item.stack_energy), 6), + } + if plot_kind and index == 0: + try: + row["rnaplot_svg"] = render_structure_svg(complex_obj, item.structure, plot_kind) + except Exception as exc: + row["rnaplot_error"] = str(exc) + output.append(row) + return output + + +def serialize_samples(items): + return [str(item) for item in (items or [])] + + +def serialize_pairs(pair_matrix, preview_limit=24): + if pair_matrix is None: + return None + + dense = numpy.asarray(pair_matrix.to_array(), dtype=float) + preview = dense[:preview_limit, :preview_limit] + return { + "shape": list(dense.shape), + "preview": numpy.round(preview, 6).tolist(), + "preview_size": int(preview_limit), + "preview_truncated": bool(dense.shape[0] > preview_limit), + } + + +def parse_structure_pairs(structure): + stack = [] + pairs = [] + compact_index = -1 + index_map = [] + + for char in str(structure): + if char == "+": + continue + compact_index += 1 + index_map.append(compact_index) + if char == "(": + stack.append(compact_index) + elif char == ")": + if not stack: + raise ValueError("Unbalanced structure: missing opening bracket.") + left = stack.pop() + pairs.append((left, compact_index)) + + if stack: + raise ValueError("Unbalanced structure: missing closing bracket.") + + pairs.sort() + return pairs, compact_index + 1 + + +def build_structure_probabilities(structure, pair_matrix): + if pair_matrix is None: + return None + + dense = numpy.asarray(pair_matrix.to_array(), dtype=float) + mfe_pairs, structure_length = parse_structure_pairs(structure) + + if dense.shape[0] != dense.shape[1]: + return None + if structure_length != dense.shape[0]: + return None + + pair_partner = {} + for left, right in mfe_pairs: + pair_partner[left] = right + pair_partner[right] = left + + residue_probabilities = [] + for index in range(structure_length): + partner = pair_partner.get(index) + if partner is None: + probability = float(dense[index, index]) + else: + probability = float(dense[index, partner]) + residue_probabilities.append(round(probability, 6)) + + pair_probabilities = [ + { + "i": left + 1, + "j": right + 1, + "probability": round(float(dense[left, right]), 6), + } + for left, right in mfe_pairs + ] + + return { + "length": structure_length, + "residue_probabilities": residue_probabilities, + "pair_probabilities": pair_probabilities, + } + + +def serialize_complex_result(complex_obj, data, pairs_preview_size=24): + payload = { + "name": getattr(complex_obj, "name", None), + "display": stringify_complex(complex_obj), + "strand_names": [strand.name for strand in complex_obj.strands], + "sequence": flatten_sequence(complex_obj), + "strand_lengths": [len(str(strand)) for strand in complex_obj.strands], + } + + if data.pfunc is not None: + payload["pfunc"] = format(data.pfunc, "g") if isinstance(data.pfunc, Decimal) else str(data.pfunc) + if data.free_energy is not None: + payload["free_energy_kcal_mol"] = round(float(data.free_energy), 6) + if data.ensemble_size is not None: + payload["ensemble_size"] = int(data.ensemble_size) + if data.mfe_stack is not None: + payload["mfe_stack_kcal_mol"] = round(float(data.mfe_stack), 6) + if data.mfe is not None: + payload["mfe"] = serialize_structures(complex_obj, data.mfe, plot_kind="mfe") + if data.subopt is not None: + payload["subopt"] = serialize_structures(complex_obj, data.subopt) + if data.sample is not None: + payload["sample"] = serialize_samples(data.sample) + if data.pairs is not None: + payload["pairs"] = serialize_pairs(data.pairs, preview_limit=pairs_preview_size) + if payload.get("mfe"): + payload["mfe"][0]["pair_probability_annotations"] = build_structure_probabilities( + payload["mfe"][0]["structure"], data.pairs + ) + + return payload + + +def attach_pairs_result(payload, pair_data, pairs_preview_size=24): + if pair_data is None or pair_data.pairs is None: + return payload + + payload["pairs"] = serialize_pairs(pair_data.pairs, preview_limit=pairs_preview_size) + if payload.get("mfe"): + payload["mfe"][0]["pair_probability_annotations"] = build_structure_probabilities( + payload["mfe"][0]["structure"], pair_data.pairs + ) + return payload + + +def serialize_target_defect_rows(df, limit=None): + if df is None: + return [] + rows = [] + for record in df.to_dict("records"): + row = {} + for key, value in record.items(): + if key in {"tube", "complex"}: + continue + if isinstance(value, numpy.floating): + row[key] = float(value) + else: + row[key] = value + rows.append(row) + if limit is not None: + rows = rows[:limit] + return rows + + +def serialize_design_result( + design_result, + design_mode, + model_summary, + design_options, + target_rows, + tube_rows, + ordered_domains, + ordered_strands, +): + designed_domains = [] + designed_domain_map = getattr(design_result, "domains", {}) or {} + for item in ordered_domains: + designed_domain = designed_domain_map.get(item["object"]) + designed_domains.append( + { + "name": item["name"], + "constraint": item["constraint"], + "sequence": str(designed_domain) if designed_domain is not None else None, + "length": len(str(designed_domain)) if designed_domain is not None else None, + } + ) + + designed_strands = [] + for item in ordered_strands: + target_strand = item["object"] + analysis_strand = design_result.to_analysis[target_strand] + designed_strands.append( + { + "name": item["name"], + "constraint": item["constraint"], + "constraint_kind": item.get("constraint_kind", "sequence_constraint"), + "definition": item.get("definition", item["constraint"]), + "sequence": str(analysis_strand), + "length": len(str(analysis_strand)), + } + ) + + target_complexes = [] + for target in target_rows: + target_complex = target["object"] + analysis_complex = design_result.to_analysis[target_complex] + target_complexes.append( + { + "name": target["name"], + "display": stringify_complex(analysis_complex), + "strand_names": list(target["strands"]), + "structure": target["structure"], + "sequence": flatten_sequence(analysis_complex), + "target_concentration_M": target.get("target_concentration_M"), + } + ) + + concentration_rows = [] + concentration_table = getattr(getattr(design_result, "concentrations", None), "table", None) + if design_mode == "tube" and concentration_table is not None: + for record in concentration_table.to_dict("records"): + concentration_rows.append( + { + "complex_name": record["complex_name"], + "concentration_M": float(record["concentration"]), + "target_concentration_M": float(record["target_concentration"]), + "nucleotides": int(record["nucleotides"]), + } + ) + concentration_rows.sort(key=lambda row: row["concentration_M"], reverse=True) + + return { + "workflow": "design", + "mode": design_mode, + "model": model_summary, + "compute": ["design"], + "options": design_options, + "strands": designed_strands, + "complexes": target_complexes, + "total_complex_count": len(target_complexes), + "displayed_complex_count": len(target_complexes), + "tube": { + "name": tube_rows[0]["name"] if tube_rows else "design_tube", + "max_size": tube_rows[0]["max_size"] if tube_rows else design_options["off_target_max_size"], + "complex_concentrations": concentration_rows[: design_options["result_limit"]], + "total_complex_concentrations": len(concentration_rows), + } if design_mode == "tube" else None, + "design": { + "ensemble_defect": float(design_result.ensemble_defect), + "domains": designed_domains, + "stats": { + key: float(value) if isinstance(value, (int, float, numpy.floating)) else value + for key, value in design_result.stats.items() + }, + "job_options": { + "f_stop": design_options["stop_condition"], + "seed": design_options["seed"], + "wobble_mutations": design_options["wobble_mutations"], + "max_time_seconds": design_options["max_time_seconds"], + }, + "target_tubes": [ + { + "name": row["name"], + "max_size": row["max_size"], + "on_targets": row["on_targets"], + } + for row in tube_rows + ], + "targets": target_complexes, + "defects": { + "tubes": serialize_target_defect_rows(getattr(design_result.defects, "tubes", None)), + "complexes": serialize_target_defect_rows(getattr(design_result.defects, "complexes", None)), + "tube_complexes": serialize_target_defect_rows(getattr(design_result.defects, "tube_complexes", None)), + }, + "concentrations": concentration_rows[: design_options["result_limit"]], + }, + } + + +def sort_concentrations(complex_concentrations): + rows = [] + for complex_obj, value in complex_concentrations.items(): + rows.append( + { + "display": stringify_complex(complex_obj), + "strand_names": [strand.name for strand in complex_obj.strands], + "concentration_M": float(value), + } + ) + rows.sort(key=lambda row: row["concentration_M"], reverse=True) + return rows + + +def run_job_payload(payload): + workflow = payload.get("workflow", "analysis") + mode = payload.get("mode", "tube") + model_input = payload.get("model") or {} + model = build_model(model_input) + model_summary = build_model_summary(model_input) + + if workflow == "design": + design_options = parse_design_options(payload) + design_job_options = DesignOptions( + f_stop=design_options["stop_condition"], + seed=design_options["seed"], + wobble_mutations=design_options["wobble_mutations"], + max_time=design_options["max_time_seconds"], + ) + design_domain_map, design_domains = build_design_domains(payload.get("design_domains") or []) + target_strand_map, design_strands = build_design_strands( + payload.get("strands") or [], + domain_map=design_domain_map, + ) + target_rows, target_complex_map = parse_design_complexes(payload, target_strand_map) + hard_constraints = build_hard_constraints(payload, design_domain_map, target_strand_map) + soft_constraints = build_soft_constraints( + payload, + design_domain_map, + target_strand_map, + target_complex_map, + ) + tube_rows = [] + design_tubes = [] + if mode == "tube": + tube_rows, design_tubes = parse_design_tubes( + payload, + target_rows, + target_complex_map, + design_options["off_target_max_size"], + ) + validate_design_object_names(design_domains, design_strands, target_rows, tube_rows) + defect_weights = build_defect_weights( + payload, + mode, + design_domains, + design_strands, + target_rows, + tube_rows, + ) + + if mode == "tube": + tube_name = ((payload.get("tube") or {}).get("name") or "design_tube").strip() or "design_tube" + design_job = tube_design( + tubes=design_tubes, + model=model, + options=design_job_options, + hard_constraints=hard_constraints, + soft_constraints=soft_constraints, + defect_weights=defect_weights, + ) + elif mode == "complex": + design_job = complex_design( + complexes=[row["object"] for row in target_rows], + model=model, + options=design_job_options, + hard_constraints=hard_constraints, + soft_constraints=soft_constraints, + defect_weights=defect_weights, + ) + else: + raise ValueError(f"Unsupported design mode: {mode}") + + results = design_job.run(trials=design_options["trials"]) + best_result = min(results, key=lambda item: float(item.ensemble_defect)) + return serialize_design_result( + best_result, + mode, + model_summary, + design_options, + target_rows, + tube_rows, + design_domains, + design_strands, + ) + + if workflow != "analysis": + raise ValueError(f"Unsupported workflow: {workflow}") + + compute = parse_compute(payload) + options = parse_options(payload) + nupack_options = { + key: value + for key, value in options.items() + if key not in {"result_limit", "pairs_preview_size"} + } + strand_map, strands = build_strands(payload.get("strands") or []) + + if mode == "tube": + wants_pairs = "pairs" in compute + tube_compute = [item for item in compute if item != "pairs"] + if not tube_compute: + tube_compute = ["pfunc"] + tube_cfg = payload.get("tube") or {} + max_size = int(tube_cfg.get("max_size", 2)) + if max_size < 1 or max_size > 8: + raise ValueError("max_size must be between 1 and 8.") + strand_concentrations = { + item["object"]: unit_to_molar(item["concentration"], item["unit"]) + for item in strands + } + tube = Tube( + strands=strand_concentrations, + complexes=SetSpec(max_size=max_size), + name=(tube_cfg.get("name") or "tube1").strip() or "tube1", + ) + result = tube_analysis([tube], model=model, compute=tube_compute, options=nupack_options) + concentration_rows = sort_concentrations(result[tube].complex_concentrations) + total_complex_count = len(concentration_rows) + result_limit = options["result_limit"] + displayed_concentration_rows = concentration_rows[:result_limit] + concentration_by_display = { + row["display"]: row["concentration_M"] for row in displayed_concentration_rows + } + pair_data_by_display = {} + if wants_pairs: + displayed_complex_objects = [ + complex_obj + for complex_obj in result.complexes + if stringify_complex(complex_obj) in concentration_by_display + ] + if displayed_complex_objects: + pair_result = complex_analysis( + displayed_complex_objects, + model=model, + compute=["pairs"], + options=nupack_options, + ) + pair_data_by_display = { + stringify_complex(complex_obj): data + for complex_obj, data in pair_result.complexes.items() + } + complexes = [] + for complex_obj, data in result.complexes.items(): + display_name = stringify_complex(complex_obj) + if display_name not in concentration_by_display: + continue + row = serialize_complex_result( + complex_obj, + data, + pairs_preview_size=options["pairs_preview_size"], + ) + if wants_pairs: + attach_pairs_result( + row, + pair_data_by_display.get(display_name), + pairs_preview_size=options["pairs_preview_size"], + ) + row["concentration_M"] = concentration_by_display[row["display"]] + complexes.append(row) + complexes.sort(key=lambda row: row["concentration_M"], reverse=True) + return { + "workflow": "analysis", + "mode": "tube", + "model": model_summary, + "compute": compute, + "options": options, + "strands": [ + { + "name": item["name"], + "sequence": item["sequence"], + "concentration": item["concentration"], + "unit": item["unit"], + "concentration_M": unit_to_molar(item["concentration"], item["unit"]), + } + for item in strands + ], + "tube": { + "name": tube.name, + "max_size": max_size, + "complex_concentrations": displayed_concentration_rows, + "total_complex_concentrations": total_complex_count, + }, + "complexes": complexes, + "total_complex_count": total_complex_count, + "displayed_complex_count": len(complexes), + } + + if mode == "complex": + complexes = parse_complex_lines(payload.get("complexes_text", ""), strand_map) + result = complex_analysis(complexes, model=model, compute=compute, options=nupack_options) + rows = [ + serialize_complex_result( + complex_obj, + data, + pairs_preview_size=options["pairs_preview_size"], + ) + for complex_obj, data in result.complexes.items() + ] + rows.sort(key=lambda row: row["display"]) + total_complex_count = len(rows) + rows = rows[: options["result_limit"]] + return { + "workflow": "analysis", + "mode": "complex", + "model": model_summary, + "compute": compute, + "options": options, + "strands": [ + { + "name": item["name"], + "sequence": item["sequence"], + } + for item in strands + ], + "complexes": rows, + "total_complex_count": total_complex_count, + "displayed_complex_count": len(rows), + } + + raise ValueError(f"Unsupported analysis mode: {mode}") + + +def prune_jobs(now=None): + now = now or time.time() + expired_ids = [] + for job_id, job in JOB_STORE.items(): + updated_at = job.get("updated_at", job.get("created_at", now)) + if now - updated_at > JOB_TTL_SECONDS: + expired_ids.append(job_id) + for job_id in expired_ids: + JOB_STORE.pop(job_id, None) + + if len(JOB_STORE) > JOB_MAX_COUNT: + keep_ids = sorted( + JOB_STORE, + key=lambda job_id: JOB_STORE[job_id].get("created_at", 0), + reverse=True, + )[:JOB_MAX_COUNT] + keep_set = set(keep_ids) + for job_id in list(JOB_STORE): + if job_id not in keep_set: + JOB_STORE.pop(job_id, None) + + +def set_job_data(job): + if redis_enabled(): + redis_client().setex(job_key(job["job_id"]), JOB_TTL_SECONDS, json.dumps(job, ensure_ascii=False)) + return + with JOB_LOCK: + prune_jobs(job.get("updated_at")) + JOB_STORE[job["job_id"]] = dict(job) + + +def update_job_data(job_id, **updates): + if redis_enabled(): + client = redis_client() + current = get_job(job_id) + if current is None: + current = {"job_id": job_id, "created_at": time.time()} + current.update(updates) + current["updated_at"] = time.time() + client.setex(job_key(job_id), JOB_TTL_SECONDS, json.dumps(current, ensure_ascii=False)) + return current + + with JOB_LOCK: + current = JOB_STORE.get(job_id) + if current is None: + current = {"job_id": job_id, "created_at": time.time()} + current.update(updates) + current["updated_at"] = time.time() + prune_jobs(current["updated_at"]) + JOB_STORE[job_id] = current + return dict(current) + + +def create_job(payload): + job_id = uuid4().hex + now = time.time() + job = { + "job_id": job_id, + "status": "queued", + "error": None, + "result": None, + "created_at": now, + "updated_at": now, + "payload": payload, + } + set_job_data(job) + log_event(f"accepted job_id={job_id} mode={payload.get('mode', 'tube')}") + + if redis_enabled(): + redis_client().lpush(JOB_QUEUE_KEY, job_id) + else: + thread = threading.Thread(target=_run_job, args=(job_id, payload), daemon=True) + thread.start() + return job_id + + +def _run_job(job_id, payload): + started_at = time.time() + log_event(f"running job_id={job_id}") + update_job_data(job_id, status="running") + if redis_enabled(): + redis_client().sadd(JOB_RUNNING_KEY, job_id) + try: + result = run_job_payload(payload) + elapsed = round(time.time() - started_at, 3) + update_job_data(job_id, status="success", result=result, payload=None, elapsed_seconds=elapsed) + log_event(f"success job_id={job_id} elapsed={elapsed}s") + except Exception as exc: + elapsed = round(time.time() - started_at, 3) + update_job_data( + job_id, + status="error", + error={ + "message": str(exc), + "traceback": traceback.format_exc(), + }, + payload=None, + elapsed_seconds=elapsed, + ) + log_event(f"error job_id={job_id} elapsed={elapsed}s message={exc}") + finally: + if redis_enabled(): + redis_client().srem(JOB_RUNNING_KEY, job_id) + + +def get_job(job_id): + if redis_enabled(): + raw = redis_client().get(job_key(job_id)) + if raw is None: + return None + job = json.loads(raw) + job.pop("payload", None) + return job + with JOB_LOCK: + prune_jobs() + job = JOB_STORE.get(job_id) + if job is None: + return None + output = dict(job) + output.pop("payload", None) + return output + + +def prune_shares(now=None): + if SHARE_MAX_COUNT < 1: + return + if redis_enabled(): + client = redis_client() + extra = client.zcard(SHARE_INDEX_KEY) - SHARE_MAX_COUNT + if extra <= 0: + return + stale_ids = client.zrange(SHARE_INDEX_KEY, 0, extra - 1) + if stale_ids: + client.delete(*(share_key(share_id) for share_id in stale_ids)) + client.zrem(SHARE_INDEX_KEY, *stale_ids) + return + + keep_ids = sorted( + SHARE_STORE, + key=lambda share_id: SHARE_STORE[share_id].get("created_at", 0), + reverse=True, + )[:SHARE_MAX_COUNT] + keep_set = set(keep_ids) + for share_id in list(SHARE_STORE): + if share_id not in keep_set: + SHARE_STORE.pop(share_id, None) + + +def create_share(record): + payload = record.get("payload") + result = record.get("result") + if not isinstance(payload, dict): + raise ValueError("Share payload must contain an input payload object.") + if not isinstance(result, dict): + raise ValueError("Share payload must contain a result object.") + + now = time.time() + share_id = uuid4().hex[:16] + item = { + "id": share_id, + "created_at": now, + "created_at_iso": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime(now)), + "payload": payload, + "result": result, + "result_summary": record.get("result_summary") or {}, + } + + if redis_enabled(): + client = redis_client() + client.set(share_key(share_id), json.dumps(item, ensure_ascii=False)) + client.zadd(SHARE_INDEX_KEY, {share_id: now}) + prune_shares(now) + return item + + with SHARE_LOCK: + SHARE_STORE[share_id] = item + prune_shares(now) + return dict(item) + + +def get_share(share_id): + share_id = str(share_id or "").strip() + if not re.fullmatch(r"[0-9a-fA-F]{8,64}", share_id): + return None + if redis_enabled(): + raw = redis_client().get(share_key(share_id)) + return json.loads(raw) if raw else None + with SHARE_LOCK: + item = SHARE_STORE.get(share_id) + return dict(item) if item else None + + +def run_worker_loop(): + if not redis_enabled(): + raise RuntimeError("Worker mode requires NP_REDIS_URL and the redis package.") + + client = redis_client() + worker_count = max(1, WORKER_CONCURRENCY) + log_event( + f"Starting worker loop on Redis queue {JOB_QUEUE_KEY} " + f"(concurrency={worker_count}, per_job_threads={max(1, PER_JOB_THREAD_LIMIT)})" + ) + + with ThreadPoolExecutor(max_workers=worker_count, thread_name_prefix="np-job") as executor: + inflight = set() + + while True: + finished = {future for future in inflight if future.done()} + if finished: + inflight.difference_update(finished) + for future in finished: + future.result() + + if len(inflight) >= worker_count: + done, _ = wait(inflight, return_when=FIRST_COMPLETED, timeout=1) + inflight.difference_update(done) + for future in done: + future.result() + continue + + item = client.brpop(JOB_QUEUE_KEY, timeout=2) + if not item: + continue + _, job_id = item + raw = client.get(job_key(job_id)) + if raw is None: + continue + job = json.loads(raw) + payload = job.get("payload") + if payload is None: + continue + inflight.add(executor.submit(_run_job, job_id, payload)) + + +EXAMPLE_PAYLOAD = { + "workflow": "analysis", + "mode": "tube", + "model": { + "material": "rna", + "ensemble": "stacking", + "celsius": 37, + "sodium": 1.0, + "magnesium": 0.0, + }, + "compute": ["pfunc", "mfe", "pairs"], + "options": { + "num_sample": 20, + "energy_gap": 1.0, + "sparsity_fraction": 1.0, + "sparsity_threshold": 0.0, + "single_mfe": False, + "result_limit": 25, + "pairs_preview_size": 24, + }, + "strands": [ + {"name": "A", "sequence": "AGUCUAGGAU", "concentration": 1.0, "unit": "uM"}, + {"name": "B", "sequence": "UUAACCCACG", "concentration": 2.0, "unit": "uM"}, + ], + "tube": {"name": "tube1", "max_size": 2}, + "complexes_text": "A\nB\nA+B", + "design": { + "trials": 1, + "result_limit": 25, + "off_target_max_size": 2, + "stop_condition": 0.02, + "seed": 0, + "wobble_mutations": False, + "max_time_seconds": 0, + }, + "design_domains": [ + {"name": "a", "sequence": "N10"}, + {"name": "b", "sequence": "N10"}, + ], + "design_complexes": [ + { + "name": "AB_target", + "strands": "A+B", + "structure": "(10+)10", + } + ], + "design_tubes": [ + { + "name": "Tube 1", + "max_size": 2, + "on_targets": [ + {"complex": "AB_target", "concentration": 1.0, "unit": "uM"}, + ], + } + ], + "hard_constraints": [], + "soft_constraints": [], + "defect_weights": [], +} + +# Use A=a and B=~a so the full-duplex target "(10+)10" is base-pair consistent. +DESIGN_TUBE_EXAMPLE_PAYLOAD = { + "workflow": "design", + "mode": "tube", + "model": { + "material": "rna", + "ensemble": "stacking", + "celsius": 37, + "sodium": 1.0, + "magnesium": 0.0, + }, + "compute": ["pfunc", "mfe", "pairs"], + "options": { + "num_sample": 20, + "energy_gap": 1.0, + "sparsity_fraction": 1.0, + "sparsity_threshold": 0.0, + "single_mfe": False, + "result_limit": 25, + "pairs_preview_size": 24, + }, + "strands": [ + {"name": "A", "sequence": "a", "concentration": 1.0, "unit": "uM"}, + {"name": "B", "sequence": "~a", "concentration": 1.0, "unit": "uM"}, + ], + "tube": {"name": "Design Tube", "max_size": 2}, + "complexes_text": "A+B", + "design": { + "trials": 1, + "result_limit": 25, + "off_target_max_size": 2, + "stop_condition": 0.05, + "seed": 1, + "wobble_mutations": False, + "max_time_seconds": 0, + }, + "design_domains": [ + {"name": "a", "sequence": "N10"}, + ], + "design_complexes": [ + { + "name": "AB_target", + "strands": "A+B", + "structure": "(10+)10", + } + ], + "design_tubes": [ + { + "name": "Design Tube", + "max_size": 2, + "on_targets": [ + {"complex": "AB_target", "concentration": 1.0, "unit": "uM"}, + ], + } + ], + "hard_constraints": [], + "soft_constraints": [], + "defect_weights": [], +} + +DESIGN_COMPLEX_EXAMPLE_PAYLOAD = { + **DESIGN_TUBE_EXAMPLE_PAYLOAD, + "mode": "complex", + "tube": {"name": "Design Tube", "max_size": 1}, + "design": { + **DESIGN_TUBE_EXAMPLE_PAYLOAD["design"], + "off_target_max_size": 1, + }, + "design_tubes": [], +} + + +def get_example_payload(query): + params = parse_qs(query) + workflow = (params.get("workflow") or ["analysis"])[0] + mode = (params.get("mode") or ["tube"])[0] + if workflow == "design": + if mode == "complex": + return DESIGN_COMPLEX_EXAMPLE_PAYLOAD + return DESIGN_TUBE_EXAMPLE_PAYLOAD + return EXAMPLE_PAYLOAD + + +class AppHandler(BaseHTTPRequestHandler): + server_version = "NPReplica/0.1" + + def do_GET(self): + parsed = urlparse(self.path) + + if parsed.path == "/": + self._respond_file(INDEX_PATH, cache_control="public, max-age=60, stale-while-revalidate=86400") + return + + if parsed.path in {"/favicon.svg", "/favicon.ico"}: + self._respond_file(FAVICON_PATH, cache_control="public, max-age=86400") + return + + if parsed.path == "/design-guide.html": + self._respond_file(GUIDE_PATH, cache_control="public, max-age=3600") + return + + if parsed.path == "/health": + stats = job_stats() + self._respond( + *json_bytes( + { + "status": "ok", + "job_backend": "redis" if redis_enabled() else "memory", + "run_mode": RUN_MODE, + "queue_depth": queue_size(), + "jobs_running": stats["running"], + "jobs_queued": stats["queued"], + "job_ttl_seconds": JOB_TTL_SECONDS, + "worker_concurrency": max(1, WORKER_CONCURRENCY), + "per_job_thread_limit": max(1, PER_JOB_THREAD_LIMIT), + "nupack_threads": int(getattr(nupack_config, "threads", PER_JOB_THREAD_LIMIT)), + "nupack_cache_gb": float(getattr(nupack_config, "cache", NUPACK_CACHE_GB)), + } + ) + ) + return + + if parsed.path == "/api/example": + self._respond(*json_bytes(get_example_payload(parsed.query))) + return + + if parsed.path.startswith("/api/jobs/"): + job_id = parsed.path.rsplit("/", 1)[-1] + job = get_job(job_id) + if job is None: + self._respond(*json_bytes({"error": "Job not found"}, status=HTTPStatus.NOT_FOUND)) + return + self._respond(*json_bytes(job)) + return + + if parsed.path.startswith("/api/shares/"): + share_id = parsed.path.rsplit("/", 1)[-1] + share = get_share(share_id) + if share is None: + self._respond(*json_bytes({"error": "Share not found"}, status=HTTPStatus.NOT_FOUND)) + return + self._respond(*json_bytes({"status": "success", "share": share})) + return + + self._respond(*json_bytes({"error": "Not found"}, status=HTTPStatus.NOT_FOUND)) + + def do_POST(self): + parsed = urlparse(self.path) + + if parsed.path != "/api/analyze": + if parsed.path == "/api/jobs": + try: + length = int(self.headers.get("Content-Length", "0")) + raw = self.rfile.read(length) + payload = json.loads(raw.decode("utf-8")) + job_id = create_job(payload) + self._respond(*json_bytes({"status": "accepted", "job_id": job_id}, status=HTTPStatus.ACCEPTED)) + except Exception as exc: + self._respond( + *json_bytes( + { + "status": "error", + "error": str(exc), + "traceback": traceback.format_exc(), + }, + status=HTTPStatus.BAD_REQUEST, + ) + ) + return + if parsed.path == "/api/shares": + try: + length = int(self.headers.get("Content-Length", "0")) + if length > SHARE_MAX_BYTES: + raise ValueError(f"Share payload is too large. Limit is {SHARE_MAX_BYTES} bytes.") + raw = self.rfile.read(length) + record = json.loads(raw.decode("utf-8")) + share = create_share(record) + self._respond( + *json_bytes( + { + "status": "success", + "share_id": share["id"], + "url": f"/?share={share['id']}", + "max_count": SHARE_MAX_COUNT, + }, + status=HTTPStatus.CREATED, + ) + ) + except Exception as exc: + self._respond( + *json_bytes( + { + "status": "error", + "error": str(exc), + "traceback": traceback.format_exc(), + }, + status=HTTPStatus.BAD_REQUEST, + ) + ) + return + self._respond(*json_bytes({"error": "Not found"}, status=HTTPStatus.NOT_FOUND)) + return + + try: + length = int(self.headers.get("Content-Length", "0")) + raw = self.rfile.read(length) + payload = json.loads(raw.decode("utf-8")) + result = run_job_payload(payload) + self._respond(*json_bytes({"status": "success", "result": result})) + except Exception as exc: + self._respond( + *json_bytes( + { + "status": "error", + "error": str(exc), + "traceback": traceback.format_exc(), + }, + status=HTTPStatus.BAD_REQUEST, + ) + ) + + def log_message(self, format_, *args): + print(f"{self.address_string()} - {format_ % args}") + + def _respond_file(self, path, cache_control): + etag = file_etag(path) + if self.headers.get("If-None-Match") == etag: + self.send_response(HTTPStatus.NOT_MODIFIED) + self.send_header("ETag", etag) + self.send_header("Cache-Control", cache_control) + self.end_headers() + return + status, content_type, body = html_bytes(path) + self._respond(status, content_type, body, cache_control=cache_control, extra_headers={"ETag": etag}) + + def _respond(self, status, content_type, body, cache_control=None, extra_headers=None): + self.send_response(status) + self.send_header("Content-Type", content_type) + self.send_header("Content-Length", str(len(body))) + self.send_header("Cache-Control", cache_control or "no-store, max-age=0") + if cache_control is None: + self.send_header("Pragma", "no-cache") + self.send_header("Expires", "0") + for key, value in (extra_headers or {}).items(): + self.send_header(key, value) + self.end_headers() + self.wfile.write(body) + + +def main(): + apply_thread_limits() + if RUN_MODE == "worker": + run_worker_loop() + return + server = ThreadingHTTPServer((HOST, PORT), AppHandler) + log_event( + f"Serving NP replica on http://{HOST}:{PORT} " + f"(jobs via {'redis' if redis_enabled() else 'memory'}, " + f"per_job_threads={max(1, PER_JOB_THREAD_LIMIT)})" + ) + server.serve_forever() + + +if __name__ == "__main__": + main() diff --git a/service/split_strand_svg.py b/service/split_strand_svg.py new file mode 100644 index 0000000..e799ecc --- /dev/null +++ b/service/split_strand_svg.py @@ -0,0 +1,345 @@ +from __future__ import annotations + +import json +import math +from dataclasses import dataclass, field +from html import escape + + +NODE_RADIUS = 10.0 +PRIMARY_SPACE = 20.0 +PAIR_SPACE = 20.0 +PADDING = 52.0 + + +def _fmt(value): + return f"{float(value):.2f}" + + +def _fmt3(value): + return f"{float(value):.3f}" + + +def _pairmap_from_structure(structure): + pair_stack = [] + dangling_stack = [] + pairs = [-1] * len(structure) + + for index, char in enumerate(structure): + if char == "(": + pair_stack.append(index) + elif char == ")": + if pair_stack: + partner = pair_stack.pop() + pairs[index] = partner + pairs[partner] = index + else: + dangling_stack.append(index) + + if pair_stack: + if len(pair_stack) != len(dangling_stack): + raise ValueError("Unbalanced structure for split-strand layout.") + for left, right in zip(pair_stack, reversed(dangling_stack)): + pairs[left] = right + pairs[right] = left + + return pairs + + +@dataclass +class _TreeNode: + children: list["_TreeNode"] = field(default_factory=list) + is_pair: bool = False + index_a: int = -1 + index_b: int = -1 + x: float = 0.0 + y: float = 0.0 + go_x: float = 0.0 + go_y: float = 0.0 + + +def _add_nodes_recursive(pairmap, root, start, end): + if start > end: + raise ValueError(f"Invalid recursive span for split-strand layout: {start}>{end}") + + if pairmap[start] == end: + child = _TreeNode(is_pair=True, index_a=start, index_b=end) + _add_nodes_recursive(pairmap, child, start + 1, end - 1) + root.children.append(child) + return + + child = _TreeNode() + cursor = start + while cursor <= end: + partner = pairmap[cursor] + if partner > cursor: + _add_nodes_recursive(pairmap, child, cursor, partner) + cursor = partner + 1 + continue + child.children.append(_TreeNode(index_a=cursor)) + cursor += 1 + root.children.append(child) + + +def _setup_coords_recursive(node, parent, start_x, start_y, go_x, go_y, flipped=False): + cross_x = -go_y + cross_y = go_x + node.go_x = go_x + node.go_y = go_y + + if len(node.children) == 1: + node.x = start_x + node.y = start_y + child = node.children[0] + next_x = start_x + go_x * PRIMARY_SPACE + next_y = start_y + (-1 if flipped else 1) * go_y * PRIMARY_SPACE + if child.is_pair or (not child.is_pair and child.index_a >= 0): + _setup_coords_recursive(child, node, next_x, next_y, go_x, go_y, flipped=flipped) + else: + _setup_coords_recursive(child, node, start_x, start_y, go_x, go_y, flipped=flipped) + return + + if not node.children: + node.x = start_x + node.y = start_y + return + + pair_count = sum(1 for child in node.children if child.is_pair) + circle_length = (len(node.children) + 1) * PRIMARY_SPACE + (pair_count + 1) * PAIR_SPACE + circle_radius = circle_length / (2 * math.pi) + length_walker = PAIR_SPACE / 2.0 + + if parent is None: + node.x = go_x * circle_radius + node.y = go_y * circle_radius + else: + node.x = parent.x + go_x * circle_radius + node.y = parent.y + (-1 if flipped else 1) * go_y * circle_radius + + for child in node.children: + length_walker += PRIMARY_SPACE + if child.is_pair: + length_walker += PAIR_SPACE / 2.0 + + rad_angle = length_walker / circle_length * 2 * math.pi - math.pi / 2.0 + if parent is None: + rad_angle -= math.pi / 2.0 + + child_x = node.x + math.cos(rad_angle) * cross_x * circle_radius + math.sin(rad_angle) * go_x * circle_radius + child_y = node.y + (-1 if flipped else 1) * math.cos(rad_angle) * cross_y * circle_radius + (-1 if flipped else 1) * math.sin(rad_angle) * go_y * circle_radius + + child_go_x = child_x - node.x + child_go_y = child_y - node.y + child_go_len = math.hypot(child_go_x, child_go_y) or 1.0 + + _setup_coords_recursive( + child, + node, + child_x, + child_y, + child_go_x / child_go_len, + (-1 if flipped else 1) * child_go_y / child_go_len, + flipped=flipped, + ) + + if child.is_pair: + length_walker += PAIR_SPACE / 2.0 + + +def _collect_coords_recursive(node, xs, ys, flipped=False): + if node.is_pair: + cross_x = -node.go_y + cross_y = node.go_x + xs[node.index_a] = node.x + cross_x * PAIR_SPACE / 2.0 + xs[node.index_b] = node.x - cross_x * PAIR_SPACE / 2.0 + ys[node.index_a] = node.y + (-1 if flipped else 1) * cross_y * PAIR_SPACE / 2.0 + ys[node.index_b] = node.y + (1 if flipped else -1) * cross_y * PAIR_SPACE / 2.0 + elif node.index_a >= 0: + xs[node.index_a] = node.x + ys[node.index_a] = node.y + + for child in node.children: + _collect_coords_recursive(child, xs, ys, flipped=flipped) + + +def _layout_positions(display_structure): + pairmap = _pairmap_from_structure(display_structure) + root = _TreeNode() + cursor = 0 + while cursor < len(pairmap): + partner = pairmap[cursor] + if partner > cursor: + _add_nodes_recursive(pairmap, root, cursor, partner) + cursor = partner + 1 + continue + root.children.append(_TreeNode(index_a=cursor)) + cursor += 1 + + xs = [0.0] * len(display_structure) + ys = [0.0] * len(display_structure) + _setup_coords_recursive(root, None, 0.0, 0.0, 0.0, 1.0, flipped=False) + _collect_coords_recursive(root, xs, ys, flipped=False) + + min_x = min(x - NODE_RADIUS for x in xs) + min_y = min(y - NODE_RADIUS for y in ys) + xs = [x - min_x for x in xs] + ys = [y - min_y for y in ys] + return pairmap, xs, ys + + +def _strand_spans(display_sequence): + spans = [] + start = None + for index, char in enumerate(display_sequence): + if char == " ": + if start is not None: + spans.append((start, index - 1)) + start = None + continue + if start is None: + start = index + if start is not None: + spans.append((start, len(display_sequence) - 1)) + return spans + + +def render_split_strands_svg(strand_sequences, structure, title=None): + if not strand_sequences: + raise ValueError("Split-strand layout requires at least one strand.") + + display_sequence = " ".join(str(sequence) for sequence in strand_sequences) + display_structure = str(structure).replace("+", " ") + + if len(display_sequence) != len(display_structure): + raise ValueError("Sequence and structure length mismatch for split-strand layout.") + + pairmap, xs, ys = _layout_positions(display_structure) + visible_indices = [index for index, char in enumerate(display_sequence) if char != " "] + if not visible_indices: + raise ValueError("Split-strand layout produced no visible residues.") + + min_x = min(xs[index] for index in visible_indices) - NODE_RADIUS - PADDING + min_y = min(ys[index] for index in visible_indices) - NODE_RADIUS - PADDING + max_x = max(xs[index] for index in visible_indices) + NODE_RADIUS + PADDING + max_y = max(ys[index] for index in visible_indices) + NODE_RADIUS + PADDING + + width = max_x - min_x + height = max_y - min_y + + def tx(index): + return xs[index] - min_x + + def ty(index): + return ys[index] - min_y + + compact_index = {} + compact_counter = 1 + for index in visible_indices: + compact_index[index] = compact_counter + compact_counter += 1 + + svg_lines = [ + '', + ( + f'' + ), + ] + if title: + svg_lines.append(f"{escape(title)}") + svg_lines.extend( + [ + ' ", + ' ", + "", + f' ', + ' ', + ] + ) + + for strand_index, (start, end) in enumerate(_strand_spans(display_sequence), start=1): + points = " ".join( + f"{_fmt3(tx(index))},{_fmt3(ty(index))}" + for index in range(start, end + 1) + ) + svg_lines.append(f' ') + + svg_lines.append(' ') + for left, right in enumerate(pairmap): + if right <= left or display_sequence[left] == " " or display_sequence[right] == " ": + continue + svg_lines.append( + f' ' + ) + svg_lines.append(" ") + + svg_lines.append(' ') + for index in visible_indices: + base = display_sequence[index] + svg_lines.append( + f' {escape(base)}' + ) + svg_lines.append(" ") + + compact_sequence = display_sequence.replace(" ", "") + basepair_rows = [ + f' {{ i: {compact_index[left]}, j: {compact_index[right]}, type: "cWW" }}' + for left, right in enumerate(pairmap) + if right > left and display_sequence[left] != " " and display_sequence[right] != " " + ] + coord_rows = [ + f' {{ x: {_fmt3(tx(index))}, y: {_fmt3(ty(index))} }}' + for index in visible_indices + ] + svg_lines.extend( + [ + " ", + '", + "", + ] + ) + return "\n".join(svg_lines) diff --git a/绘图自动分割.txt b/绘图自动分割.txt new file mode 100644 index 0000000..94abb75 --- /dev/null +++ b/绘图自动分割.txt @@ -0,0 +1,94 @@ + 在ununtu上使用该命令安装draw_rna +pip install matplotlib numpy draw_rna --upgrade -i https://pypi.tuna.tsinghua.edu.cn/simple + +-------------------------------打开原本步骤中lab中的记事本-------------------------- +# 依赖导入 +from nupack import * +from draw_rna.ipynb_draw import draw_struct +import os +import matplotlib.pyplot as plt + +# 开启显示（恢复正常绘图展示） +plt.ion() + +# ====================================== +# 仅保留两个核心函数，完全用 draw_rna +# ====================================== +def get_mfe_structure(complex_name: str, nupack_result, strands: list): + """获取复合物MFE结构，支持同源二聚体 Y1+Y1""" + strand_names = complex_name.strip('()').split('+') + strand_map = {s.name: s for s in strands} + strand_objs = [strand_map[name] for name in strand_names] + complex_obj = Complex(strand_objs) + return str(nupack_result[complex_obj].mfe[0].structure) + +def draw_complex(complex_name: str, nupack_result, strands: list, output_folder="./"): + """draw_rna绘图：无页面显示，直接保存图片到目标路径""" + os.makedirs(output_folder, exist_ok=True) + + # 获取结构和序列（完全不变） + struct = get_mfe_structure(complex_name, nupack_result, strands) + strand_names = complex_name.strip('()').split('+') + strand_map = {s.name: s for s in strands} + seq = " ".join([str(strand_map[name]) for name in strand_names]) # 序列用空格分隔 + draw_structure = struct.replace('+', ' ') # 结构用空格分隔 + # 设置全局字体大小，可以根据需要调整数字 (例如 14, 18, 20) + plt.rcParams['font.size'] = 12 + # ---------- 修正点：显式创建画布并传递给 draw_struct ---------- + fig, ax = plt.subplots(figsize=(15, 15)) # 可调尺寸，保证图形清晰 + draw_struct(seq, draw_structure, ax=ax) # 在指定 ax 上绘图 + + # 保存并关闭 + save_name = complex_name.strip('()').replace('+', '_') + '.png' + save_path = os.path.join(output_folder, save_name) + plt.savefig(save_path, dpi=300, bbox_inches='tight') + + print(f"\n✅ 图片已保存至：{save_path}") + print(f"MFE 结构：{struct}") +def draw_all_complexes(nupack_result, strands: list, output_folder="./"): + print("\n========== 开始批量绘制所有复合物 ==========") + + # 过滤出 Complex 对象 + complex_objs = [obj for obj in nupack_result.keys() if not isinstance(obj, Tube)] + total = len(complex_objs) + success_count = 0 + + for idx, complex_obj in enumerate(complex_objs, 1): + complex_name = complex_obj.name # ✅ 修正点 + try: + draw_complex(complex_name, nupack_result, strands, output_folder) + success_count += 1 + print(f"进度：{idx}/{total} 完成") + except Exception as e: + print(f"❌ 绘制复合物 {complex_name} 失败: {e}") + + print("\n========== 所有复合物绘制完成 ==========") + print(f"✅ 批量绘制完成！成功：{success_count} / 总计：{total}") + print(f"📁 图片保存在：{os.path.abspath(output_folder)}") +-------------------------------- nupack计算-------------------------------------- + +model1 = Model(material='dna', celsius=37,sodium=0.1, magnesium=0.02) +Y1 = Strand('CGTTAACGCAGTGAGGACGGTAGTTTGTCGTTCCATCGCACC', name='Y1') +Y2 = Strand('CGTTAACGGGTGCGATGGAACGACTTTCGACAGGCCTGGTGTAATTTCACCCATGTTAGTCGA', name='Y2') +Y3 = Strand('ATACGGAAAATGGAGATAGGAAGAGTACAATGTCAGCGATAAATTCCGTATACGACACCAGGCCTGTCGATTACTACCGTCCTCACTG', name='Y3') + +t1 = Tube({Y1: 1e-8, Y2: 1e-8, Y3:1e-8 },complexes=SetSpec(max_size=4), name='Tube 1') + +my_result = tube_analysis(tubes=[t1], model=model1, + compute=['pfunc', 'pairs', 'mfe', 'sample', 'subopt'], + options={'num_sample': 1, 'energy_gap': 0.5}) +print(my_result) + +--------------------------- --------------- 单复合物绘图------------------------------ +target = '(Y1+Y1+Y2+Y3)' # 替换为你的目标复合物名称 +all_strands = [Y1, Y2, Y3] # 所有链的列表 + +# 1. 获取并打印MFE结构 +mfe_struct = get_mfe_structure(target, my_result, all_strands) +#print("MFE结构：", mfe_struct) + +# 2. 绘制二级结构 +draw_complex(target, my_result, all_strands, output_folder="/home/zsy") +--------------------------------------------批量绘图----------------------------------- +strands_list = [Y1, Y2, Y3] +draw_all_complexes(my_result, strands_list, output_folder="/home/zsy/jihe")