Varicella-Zoster Virus (VZV) in complex skin organoids is about using advanced 3D human skin models to understand how VZV infects, spreads, and interacts with host immunity in its most relevant organ: the skin. 这些类器官尽量重建了人皮肤的层次结构和细胞多样性，使得对VZV的研究比传统二维细胞培养更接近真实感染情境。¹²³⁴

为什么皮肤和类器官重要

VZV 的致病性高度依赖其在表皮角质形成细胞中的复制，这最终形成充满高滴度游离病毒的水疱，是人际传播和建立潜伏感染的关键步骤。 传统小鼠模型难以完全模拟人特异性VZV感染，因此发展成人皮肤组织模型、皮肤器官培养以及皮肤类器官对于研究病毒复制、皮损形成和疫苗株减毒机制至关重要。⁵⁶³⁷¹

复杂皮肤类器官是什么

复杂皮肤类器官通常来源于人多能干细胞或皮肤前体细胞，经过定向分化形成多层表皮、真皮样结构，有时还包括毛囊、皮脂腺甚至神经和免疫细胞成分。 这些3D结构提供了立体组织环境和细胞间相互作用，可以更真实地研究VZV在不同表皮层的扩散、跨层传播以及向神经末梢的进入过程。²³⁸⁴¹

在类器官中研究VZV复制与传播

类器官允许模拟病毒从基底层角质形成细胞开始复制，再向上层分化细胞扩展的过程，这与患者皮损中观察到的VZV复制模式高度一致。 通过在类器官中接种标记VZV，研究者可以动态监测病毒如何在组织内扩散、形成局灶性病灶，以及是否存在类似体内那样的“离散水疱”样病变。⁶³⁹⁷⁵¹

免疫控制与固有限制因子

皮肤类器官为研究VZV如何与局部固有免疫和限制因子（如干扰素通路、自噬、PML核小体等）互作提供平台。 通过在类器官中操纵这些信号（例如阻断或增强I型干扰素、改变角质形成细胞分化相关通路），可以评估这些因素如何影响病毒复制强度、病灶大小以及病毒向神经末梢的传播能力。³¹⁰⁹⁴¹

向神经系统的连接与潜伏相关模型

VZV最终在感觉神经节建立潜伏，因此将皮肤类器官与人源神经元或神经类器官连接，是当前研究中的一个重要发展方向。 一些工作利用人神经元或神经类器官模型，已经显示可以在体外建立VZV潜伏和再激活系统，未来与皮肤类器官耦合可更好模拟“皮肤–神经–神经节”轴上的完整生命周期。^1112413

与其他模型相比的优势和局限

与单层细胞相比，皮肤类器官在细胞多样性、分化梯度和屏障特性方面更接近真实皮肤，因此对药物筛选、毒力基因鉴定以及疫苗株特性分析更具生理相关性。 不过，当前类器官通常缺乏完整血管系统和成熟适应性免疫成分，这限制了对全身免疫反应和T细胞介导控制的模拟，需要与动物模型或临床数据结合解读。¹⁴⁴⁵¹

如果你愿意，可以继续围绕这个主题细化几个环节，例如：

• “如何在皮肤类器官中设计VZV基因功能筛选实验？”
• “如何把类器官数据嵌入到系统层面的病毒–宿主互作网络分析？” ^151617181920

模型系统与方法学

• 您所使用的复杂皮肤类器官在多大程度上重建了体内人皮肤的关键特征，例如分层结构、屏障形成、神经支配以及驻留免疫细胞？这些模型在解释VZV传播和免疫逃逸数据时目前主要的局限在哪里？
• 您是如何将经典病毒学方法与对受感染皮肤和神经类器官的多组学分析相结合，用于刻画VZV–宿主相互作用的时间动态？在区分“早期限制”与“成功免疫逃逸”方面，哪些读出指标被证明最具信息量？
• 在稳定感染复杂皮肤和周围神经系统类器官时，使用报告病毒时遇到的主要技术挑战是什么？您如何验证类器官中观察到的传播模式以及潜伏/再激活特征能够真实反映人皮肤和感觉神经节中的情况？

皮肤与神经中的传播

• 在您的比较分析中，皮肤类器官与神经类器官之间是否存在不同的细胞间传播方式或动力学？是否有特定的病毒基因或宿主通路在这两个解剖部位对传播发挥了差异化的调控作用？
• 在类器官不同表皮层中，病毒基因表达谱（例如基底层的立即早期基因与上层角化细胞的晚期糖蛋白表达）与体外或体内人皮肤模型中的描述有多接近？这对确定最有效的干预靶点有何启示？
• 与以往的皮肤器官培养或SCID-hu小鼠模型相比，类器官模型是否帮助澄清了VZV何时以及通过何种方式进入皮肤中的感觉神经末梢并建立潜伏感染？

先天限制因子与固有免疫

• 您提到已有工具用于研究POL III和PML核小体等限制因子；在皮肤微环境的类器官感染中，这些因子与VZV之间的相互作用有何新发现？VZV采用了哪些机制在局部对这些固有防御进行反制？
• 在您的筛选中，是否发现了特异于表皮或特异于神经元的关键抗病毒蛋白？它们的表达模式是否能解释体内水痘皮损呈离散分布且单个皮损大小受限的现象？
• 在类器官中阻断I型干扰素信号后，皮损大小及病毒产生量与既往异种移植实验相比有何差异？这一结果是否提示了VZV针对某些新的候选通路进行免疫逃逸？

免疫逃逸与免疫调控

• 在皮肤和神经类器官中，VZV是否同样能重现hiPSC来源神经球模型中报道的那种对干扰素信号和抗原呈递的深度抑制，例如ISG诱导减弱以及MHC II相关基因（如CD74）下调？
• 有研究显示VZV可以重塑T细胞表面标志促进其向皮肤归巢；在您的类器官微环境中，是否观察到类似的趋化因子、整合素或黏附分子改变，从而有利于组织内传播？
• 与传统二维角质形成细胞培养相比，在三维皮肤类器官中，VZV对表皮分化通路（例如K10降解以及对kallikrein家族的调控）的影响有何不同？这是否改变了我们对水疱形成及病毒排出的认识？

病毒基因与分子机制

• 通过类器官感染，对于特定病毒蛋白（如与PML核小体互作的ORF61）在决定皮肤趋向性、皮损形态以及溶菌性感染与免疫识别之间平衡中的作用，有哪些新的机制性见解？
• 近期多蛋白组学研究绘制了大量VZV–宿主蛋白互作网络；在类器官中，您验证了其中哪些互作具有功能意义？是否有出乎意料的宿主靶点被发现是调控传播或免疫逃逸的关键节点？
• 是否存在这样的病毒功能：在二维细胞培养中似乎可有可无，但在复杂皮肤类器官中却被证明对有效复制或传播是必需的，从而提示存在皮肤特异的毒力因子？

转化应用与未来方向

• 您如何展望利用皮肤和神经类器官作为VZV下一代疫苗或抗病毒药物的前临床评估平台，特别是那些旨在增强局部先天免疫应答或阻断特定免疫逃逸机制的策略？
• 鉴于人群中VZV相关疾病存在显著个体差异，您是否计划构建患者来源或基因分型分层的类器官面板，用以研究宿主遗传背景如何塑造限制因子表达以及对免疫逃逸的易感性？
• 在您看来，就VZV皮肤发病机制和免疫控制的若干悬而未决问题而言，类器官模型在未来几年最有优势回答的是哪些问题？哪些问题仍然必须依赖体内或临床研究来解决？

如果你希望，我也可以把这些问题进一步压缩成适合现场快速提问的精简版要点。 ^11023456789

Here are example seminar questions you could ask, focused on newly advanced work using complex skin (and neuronal) organoids for VZV spread, immune control, and evasion.

Model system and methodology

• How do the complex skin organoids used in your project recapitulate key in vivo features of human skin, such as stratification, barrier formation, innervation and resident immune cells, and what are the main limitations that still affect interpretation of VZV spread and immune evasion data?¹³
• Could you elaborate on how you combined classical virology with multi‑omics in infected skin and neuronal organoids to map VZV–host interactions over time, and which readouts turned out to be most informative for distinguishing early restriction from successful immune evasion?¹⁰¹
• What were the main technical challenges in stably infecting complex skin and PNS organoids with reporter VZV, and how did you validate that the observed spread patterns and latency/reactivation features reflect the situation in human skin and sensory ganglia?³¹

Viral spread in skin and neurons

• In your comparative analysis, did you observe distinct modes or kinetics of cell‑to‑cell spread in skin versus neuronal organoids, and are there specific viral genes or host pathways that differentially control spread in these two compartments?¹³
• How closely does the pattern of viral gene expression across epidermal layers in organoids (immediate‑early in basal cells versus late glycoproteins in suprabasal layers) mirror what has been described in ex vivo or in vivo human skin, and what does that tell us about where intervention would be most effective?⁸³
• Have organoid models helped to clarify when and how VZV gains access to sensory nerve endings in skin and then establishes latency, compared with earlier skin organ culture or SCID‑hu mouse data?⁴³⁸

Innate restriction factors and intrinsic immunity

• You mentioned established tools to study POL III and PML nuclear bodies as restriction factors; what have organoid infections revealed about how VZV counters these particular intrinsic defenses in the skin microenvironment?⁶¹
• Are there epidermis‑specific or neuron‑specific antiviral proteins that emerged as key restriction factors from your screens, and how do their expression patterns correlate with discrete lesion formation and limited lesion size seen in vivo?³¹
• How does blocking type I interferon signaling in organoids affect lesion size and viral yield compared with earlier xenograft experiments, and does this reveal any new candidate pathways that VZV targets for immune evasion?⁸³

Immune evasion and modulation

• To what extent do skin and neuronal organoids reproduce the profound suppression of interferon signaling and antigen presentation that has been reported in hiPSC‑derived neurospheroids, for example reduced ISG induction and MHC class II–associated genes such as CD74?⁵⁷¹
• Some studies show that VZV remodels T‑cell surface markers to promote skin homing; have similar VZV‑driven changes in chemokines, integrins or adhesion molecules been observed within the organoid microenvironment that might facilitate intra‑tissue spread?¹⁸
• How does VZV manipulation of epidermal differentiation pathways (e.g., degradation of K10 and modulation of kallikreins) look in your 3D organoids compared with traditional monolayer keratinocyte cultures, and does this change how we think about blister formation and virus shedding?³

Viral genes and molecular mechanisms

• What new insights have organoid infections given into the role of specific viral proteins such as ORF61 (and its interaction with PML nuclear bodies) in determining skin tropism, lesion morphology, and the balance between lytic replication and immune detection?⁶⁸
• Recent multi‑proteomic surveys mapped dozens of VZV–host protein interactions; which of these interactions have you been able to validate functionally in organoids, and did any unexpected host targets emerge as central nodes for controlling spread or immune evasion?¹⁰¹
• Have you identified viral functions that appear dispensable in 2D cell culture but are clearly required for efficient replication or spread in complex skin organoids, suggesting skin‑specific virulence determinants?⁶¹

Translational and future directions

• How do you envision using skin and neuronal organoids as preclinical platforms to test next‑generation VZV vaccines or antivirals, particularly those aimed at enhancing local innate responses or blocking specific immune‑evasive mechanisms?⁴⁸¹
• Given the inter‑individual variability in human VZV disease, do you plan to develop patient‑derived or genetically stratified organoid panels to study how host genetics shapes restriction factor expression and susceptibility to immune evasion?¹⁰¹
• Which open questions about VZV skin pathogenesis and immune control do you think organoid models are best positioned to answer in the next few years, and where will in vivo or clinical studies still be indispensable?⁸¹³

²⁹

根据德国《居留法》（AufenthG）第51条第2款，针对在德居住超过15年或配偶是德国人的长居持有者，存在豁免“离境6个月自动失效”的特别规定。

1 针对在德国合法居住超过15年的人士

如果您持有长居（Niederlassungserlaubnis），且同时满足以下两个条件，离开德国后长居不会失效：

居住时间： 您已在德国合法居住至少 15年（时间计算从首次获得合法居留算起，不仅是拿到长居后的时间）。

生活保障（最关键点）：
您的生活来源必须是有保障的（Lebensunterhalt gesichert）。这意味着即使您长期居住在国外，也必须证明有足够的资金支持（如养老金、资产收入、存款等）。

注意： 如果您依赖德国的社会救济金（如Bürgergeld），则不享受此豁免。

2 针对德国公民的配偶

如果您持有长居，且满足以下条件，长居在离境后同样不会失效：

配偶身份：您的配偶必须拥有德国国籍。 婚姻关系：双方必须处于婚姻存续状态（eheliche Lebensgemeinschaft）。

特别提醒： 如果婚姻关系破裂（离婚或正式分居），此豁免权可能会立即终止，届时将重新适用“6个月离境失效”的规定。

!!! 至关重要的操作步骤 申请“不失效证明” !!!

虽然法律有此豁免，但系统可能会因为您长期离境自动判定长居失效。为了避免将来入境被拒或引起不必要的麻烦，请务必在离境前完成以下步骤：

前往外管局：

在离开德国前，联系居住地的外管局（Ausländerbehörde）。

提交申请： 申请开具一份“长居许可不失效证明” (Bescheinigung über das Nichterlöschen der Niederlassungserlaubnis)。

携带材料： 通常需要护照、长居卡、居住满15年的证明（或结婚证及配偶德国护照）、以及资金/养老金证明。

获得凭证： 拿到书面证明后，您即可放心长期离境。将来返回德国时，请出示该证明和护照。

法律依据： § 51 Abs. 2 Satz 1 & 2 AufenthG

德语申请信模版

您的姓名/现居地址/电话/邮箱

An die Ausländerbehörde (外管局地址，如不知道具体部门，写具体城市名即可，例如 Stadt Frankfurt am Main)

(日期: TT.MM.JJJJ)

Betreff: Antrag auf Ausstellung einer Bescheinigung über das Nichterlöschen der Niederlassungserlaubnis gemäß § 51 Abs. 2 AufenthG Mein Aktenzeichen / Geburtsdatum: (您的档案号，如果有的话，或者写出生日期)

Sehr geehrte Damen und Herren,

hiermit beantrage ich eine Bescheinigung darüber, dass meine Niederlassungserlaubnis auch bei einem geplanten Auslandsaufenthalt von mehr als sechs Monaten nicht erlischt.

Ich plane, mich ab dem (预计离境日期) für längere Zeit im Ausland aufzuhalten.

Begründung: (请从以下 A 和 B 中二选一，删除不适用的那个)

(选项 A：如果您在德居住已超过15年) Gemäß § 51 Abs. 2 Satz 1 AufenthG erlischt die Niederlassungserlaubnis nicht, da ich mich seit mindestens 15 Jahren rechtmäßig im Bundesgebiet aufhalte (seit 入德年份) und mein Lebensunterhalt gesichert ist.

(选项 B：如果您的配偶是德国人) Gemäß § 51 Abs. 2 Satz 2 AufenthG erlischt die Niederlassungserlaubnis nicht, da ich mit einem deutschen Staatsangehörigen, (配偶姓名), in ehelicher Lebensgemeinschaft lebe.

Anbei übersende ich Ihnen folgende Unterlagen zum Nachweis:

Kopie meines Reisepasses und der Niederlassungserlaubnis

(针对选项A) Nachweis über die Sicherung des Lebensunterhalts (z.B. Rentenbescheid, Vermögensnachweis) (中文备注：生活来源证明，如退休金通知、资产证明)

(针对选项B) Kopie der Heiratsurkunde und des Personalausweises meines Ehepartners (中文备注：结婚证复印件及配偶身份证复印件)

Ich bitte um eine schriftliche Bestätigung. Für Rückfragen stehe ich Ihnen gerne zur Verfügung.

Mit freundlichen Grüßen

签名

姓名打印体

使用小贴士：

二选一： 模版中“Begründung”部分有两个选项，请务必删除不符合您情况的那一段，只保留一段。

附件 (Anlagen)：

如果是15年规定：重点在于证明您“有钱”，不需要领救济金。请附上养老金证明（Rentenbescheid）、银行存款证明或房产证明等。

如果是配偶规定：重点在于证明“婚姻关系”和“配偶国籍”。请附上结婚证和配偶的德国身份证/护照复印件。

发送方式： 建议先通过电子邮件发送给外管局，如果没有回复，再通过挂号信（Einschreiben）邮寄，以确保有据可查。

1.

    #Set the comparison condition in the R-script
    comparison <- "Moxi_18h_vs_Untreated_18h"
    # e.g.
    # comparison <- "Mitomycin_18h_vs_Untreated_18h"
    # comparison <- "Mitomycin_8h_vs_Untreated_8h"
    # comparison <- "Moxi_8h_vs_Untreated_8h"
    # comparison <- "Mitomycin_4h_vs_Untreated_4h"
    # comparison <- "Moxi_4h_vs_Untreated_4h"
    (r_env) Rscript make_circos_from_deseq.R

    (circos-env) cd circos_Moxi_18h_vs_Untreated_18h
    (circos-env) touch karyotype.txt (see the step2)
    (circos-env) touch circos.conf (see the step3)
    cp circos_Moxi_18h_vs_Untreated_18h/karyotype.txt circos_Moxi_8h_vs_Untreated_8h/
    cp circos_Moxi_18h_vs_Untreated_18h/karyotype.txt circos_Moxi_4h_vs_Untreated_4h/
    cp circos_Moxi_18h_vs_Untreated_18h/karyotype.txt circos_Mitomycin_18h_vs_Untreated_18h/
    cp circos_Moxi_18h_vs_Untreated_18h/karyotype.txt circos_Mitomycin_8h_vs_Untreated_8h/
    cp circos_Moxi_18h_vs_Untreated_18h/karyotype.txt circos_Mitomycin_4h_vs_Untreated_4h/
    cp circos_Moxi_18h_vs_Untreated_18h/circos.conf circos_Moxi_8h_vs_Untreated_8h/
    cp circos_Moxi_18h_vs_Untreated_18h/circos.conf circos_Moxi_4h_vs_Untreated_4h/
    cp circos_Moxi_18h_vs_Untreated_18h/circos.conf circos_Mitomycin_18h_vs_Untreated_18h/
    cp circos_Moxi_18h_vs_Untreated_18h/circos.conf circos_Mitomycin_8h_vs_Untreated_8h/
    cp circos_Moxi_18h_vs_Untreated_18h/circos.conf circos_Mitomycin_4h_vs_Untreated_4h/

    (circos-env) circos -conf circos.conf
    #or (circos-env) /home/jhuang/mambaforge/envs/circos-env/bin/circos  -conf /mnt/md1/DATA/Data_JuliaFuchs_RNAseq_2025/results/star_salmon/degenes/circos/circos.conf

    #(circos-env) jhuang@WS-2290C:/mnt/md1/DATA/Data_JuliaFuchs_RNAseq_2025/results/star_salmon/degenes$ find . -name "*circos.png"
    mv ./circos_Moxi_4h_vs_Untreated_4h/circos.png circos_Moxi_4h_vs_Untreated_4h.png
    mv ./circos_Moxi_18h_vs_Untreated_18h/circos.png circos_Moxi_18h_vs_Untreated_18h.png
    mv ./circos_Moxi_8h_vs_Untreated_8h/circos.png circos_Moxi_8h_vs_Untreated_8h.png
    mv ./circos_Mitomycin_18h_vs_Untreated_18h/circos.png circos_Mitomycin_18h_vs_Untreated_18h.png
    mv ./circos_Mitomycin_8h_vs_Untreated_8h/circos.png circos_Mitomycin_8h_vs_Untreated_8h.png
    mv ./circos_Mitomycin_4h_vs_Untreated_4h/circos.png circos_Mitomycin_4h_vs_Untreated_4h.png

2. touch karyotype.txt

    chr - CP052959.1 CP052959.1 0 2706926 grey

3. touch circos.conf

    karyotype = karyotype.txt

<ideogram>

<spacing>
    default = 0.005r
    </spacing>

    radius       = 0.80r
    thickness    = 20p
    fill         = yes

    show_label   = yes
    label_radius = 1.05r
    label_size   = 30p
    label_font   = bold
    label_parallel = yes
    </ideogram>

    # --- Wichtig: Schalter auf Top-Level, NICHT im 
<ticks>-Block ---
    show_ticks       = yes
    show_tick_labels = yes

<ticks>
    # Starte direkt an der äußeren Ideogramm-Kante
    radius      = dims(ideogram,radius_outer)
    orientation = out          # Ticks nach außen zeichnen (oder 'in' für nach innen)
    color       = black
    thickness   = 2p
    size        = 8p

    # kleine Ticks alle 100 kb, ohne Label

<tick>
    spacing      = 50000
    size         = 8p
    thickness    = 3p
    color        = black
    show_label   = no
    </tick>

    # große Ticks alle 500 kb, mit Label in Mb

<tick>
    spacing      = 100000
    size         = 12p
    thickness    = 5p
    color        = black

    show_label   = yes
    label_size   = 18p
    label_offset = 6p
    multiplier   = 0.000001     # in Mb
    format       = %.1f
    suffix       = " Mb"
    </tick>

    </ticks>

<plots>

    # Density of up-regulated genes

<plot>
    show      = yes
    type      = histogram
    file      = data/density_up.txt
    r0        = 0.88r
    r1        = 0.78r
    fill_color = red
    thickness = 1p
    </plot>

    # Density of down-regulated genes

<plot>
    show      = yes
    type      = histogram
    file      = data/density_down.txt
    r0        = 0.78r
    r1        = 0.68r
    fill_color = blue
    thickness = 1p
    </plot>

    # Scatter of individual significantly DE genes

<plot>
    show           = yes
    type           = scatter
    file           = data/genes_scatter.txt
    r0             = 0.46r
    r1             = 0.76r
    glyph          = circle
    glyph_size     = 5p
    stroke_thickness = 0

    min            = -15
    max            =  15

<rules>

<rule>
    condition = var(value) > 0
    color     = red
    </rule>

<rule>
    condition = var(value) < 0
    color     = blue
    </rule>
    </rules>

    </plot>

    </plots>

<image>
    <<include etc/image.conf>>
    </image>

    <<include etc/colors_fonts_patterns.conf>>
    <<include etc/housekeeping.conf>>

4. Rscript make_circos_from_deseq.R

   ############################################################
   # make_circos_from_deseq.R
   #
   # - Read DESeq2 results (annotated CSV)
   # - Read BED annotation
   # - Merge by GeneID
   # - Classify genes (up / down / ns)
   # - Create Circos input files:
   #     * genes_scatter.txt
   #     * density_up.txt
   #     * density_down.txt
   # - Create annotated versions of the above
   # - Export annotated tables into one Excel workbook
   ############################################################
   
   ###############################
   # 0. Single parameter to set
   ###############################
   
   # Just change this line for each comparison:
   #comparison <- "Moxi_18h_vs_Untreated_18h"
   # e.g.
   #comparison <- "Mitomycin_18h_vs_Untreated_18h"
   #comparison <- "Mitomycin_8h_vs_Untreated_8h"
   #comparison <- "Moxi_8h_vs_Untreated_8h"
   #comparison <- "Mitomycin_4h_vs_Untreated_4h"
   comparison <- "Moxi_4h_vs_Untreated_4h"
   
   ###############################
   # 0. Derived settings
   ###############################
   
   # DESeq2 result file (annotated)
   deseq_file <- paste0(comparison, "-all_annotated.csv")
   
   # BED file with gene coordinates (constant for this genome)
   bed_file   <- "CP052959_m.bed"
   
   # Create a filesystem-friendly directory name from comparison
   safe_comp  <- gsub("[^A-Za-z0-9_]+", "_", comparison)
   
   # Output directory for Circos files (e.g. circos_Moxi_18h_vs_Untreated_18h)
   out_dir    <- paste0("circos_", safe_comp)
   
   # Thresholds for significance
   padj_cutoff <- 0.05
   lfc_cutoff  <- 1      # |log2FC| >= 1
   
   # Bin size for density in bp (e.g. 10000 = 10 kb)
   bin_size <- 10000
   
   # Comparison label for Excel / plots (human-readable)
   comparison_label <- comparison
   
   options(scipen = 999) # turn off scientific notation
   
   ###############################
   # 1. Setup & packages
   ###############################
   
   if (!dir.exists(out_dir)) dir.create(out_dir, showWarnings = FALSE)
   data_dir <- file.path(out_dir, "data")
   if (!dir.exists(data_dir)) dir.create(data_dir, showWarnings = FALSE)
   
   # openxlsx for Excel export
   if (!requireNamespace("openxlsx", quietly = TRUE)) {
   stop("Package 'openxlsx' is required. Please install it with install.packages('openxlsx').")
   }
   library(openxlsx)
   
   ###############################
   # 2. Read data
   ###############################
   
   message("Reading DESeq2 results from: ", deseq_file)
   deseq <- read.csv(deseq_file, stringsAsFactors = FALSE)
   
   # Check required columns in DESeq2 table
   required_cols <- c("GeneID", "log2FoldChange", "padj")
   if (!all(required_cols %in% colnames(deseq))) {
   stop("DESeq2 table must contain columns: ", paste(required_cols, collapse = ", "))
   }
   
   message("Reading BED annotation from: ", bed_file)
   bed_cols <- c("chr","start","end","gene_id","score",
               "strand","thickStart","thickEnd",
               "itemRgb","blockCount","blockSizes","blockStarts")
   
   annot <- read.table(
   bed_file,
   header = FALSE,
   sep = "\t",
   stringsAsFactors = FALSE,
   col.names = bed_cols
   )
   
   ###############################
   # 3. Merge DESeq2 + annotation
   ###############################
   
   # Match DESeq2 GeneID (e.g. "gene-HJI06_09365") to BED gene_id
   merged <- merge(
   deseq,
   annot[, c("chr", "start", "end", "gene_id")],
   by.x = "GeneID",
   by.y = "gene_id"
   )
   
   if (nrow(merged) == 0) {
   stop("No overlap between DESeq2 GeneID and BED gene_id.")
   }
   
   # Midpoint of each gene
   merged$mid <- round((merged$start + merged$end) / 2)
   
   ###############################
   # 4. Classify genes
   ###############################
   
   merged$regulation <- "ns"
   merged$regulation[merged$padj < padj_cutoff & merged$log2FoldChange >=  lfc_cutoff]  <- "up"
   merged$regulation[merged$padj < padj_cutoff & merged$log2FoldChange <= -lfc_cutoff] <- "down"
   
   table_reg <- table(merged$regulation)
   message("Regulation counts: ", paste(names(table_reg), table_reg, collapse = " | "))
   
   ###############################
   # 5. Scatter files (per-gene)
   ###############################
   
   scatter <- merged[merged$regulation != "ns", ]
   
   # Circos scatter input: chr  start  end  value
   scatter_file <- file.path(data_dir, "genes_scatter.txt")
   scatter_out  <- scatter[, c("chr", "mid", "mid", "log2FoldChange")]
   
   write.table(scatter_out,
               scatter_file,
               quote = FALSE, sep = "\t",
               row.names = FALSE, col.names = FALSE)
   
   message("Written Circos scatter file: ", scatter_file)
   
   # Annotated scatter for Excel / inspection
   scatter_annot <- scatter[, c(
   "chr",
   "mid",              # start
   "mid",              # end
   "log2FoldChange",
   "GeneID",
   "padj",
   "regulation"
   )]
   colnames(scatter_annot)[1:4] <- c("chr", "start", "end", "log2FoldChange")
   
   scatter_annot_file <- file.path(data_dir, "genes_scatter_annotated.tsv")
   write.table(scatter_annot,
               scatter_annot_file,
               quote = FALSE, sep = "\t",
               row.names = FALSE, col.names = TRUE)
   
   message("Written annotated scatter file: ", scatter_annot_file)
   
   ###############################
   # 6. Density function (bins)
   ###############################
   
   bin_chr <- function(df_chr, bin_size, direction = c("up", "down")) {
   direction <- match.arg(direction)
   
   chr_name <- df_chr$chr[1]
   max_pos  <- max(df_chr$mid)
   
   # number of bins
   n_bins <- ceiling((max_pos + 1) / bin_size)
   
   starts <- seq(0, by = bin_size, length.out = n_bins)
   ends   <- starts + bin_size
   
   # init counts & gene lists
   counts    <- integer(n_bins)
   gene_list <- vector("list", n_bins)
   
   df_dir <- df_chr[df_chr$regulation == direction, ]
   
   if (nrow(df_dir) > 0) {
       # bin index for each gene
       bin_index <- floor(df_dir$mid / bin_size) + 1
       bin_index[bin_index < 1]        <- 1
       bin_index[bin_index > n_bins]   <- n_bins
   
       # accumulate counts and GeneIDs
       for (i in seq_len(nrow(df_dir))) {
       idx <- bin_index[i]
       counts[idx] <- counts[idx] + 1L
       gene_list[[idx]] <- c(gene_list[[idx]], df_dir$GeneID[i])
       }
   }
   
   gene_ids <- vapply(
       gene_list,
       function(x) {
       if (length(x) == 0) "" else paste(unique(x), collapse = ";")
       },
       character(1)
   )
   
   data.frame(
       chr      = chr_name,
       start    = as.integer(starts),
       end      = as.integer(ends),
       value    = as.integer(counts),
       gene_ids = gene_ids,
       stringsAsFactors = FALSE
   )
   }
   
   ###############################
   # 7. Density up/down for all chromosomes
   ###############################
   
   chr_list <- split(merged, merged$chr)
   
   density_up_list   <- lapply(chr_list, bin_chr, bin_size = bin_size, direction = "up")
   density_down_list <- lapply(chr_list, bin_size = bin_size, FUN = bin_chr, direction = "down")
   
   density_up   <- do.call(rbind, density_up_list)
   density_down <- do.call(rbind, density_down_list)
   
   # Plain Circos input (no gene_ids)
   density_up_file   <- file.path(data_dir, "density_up.txt")
   density_down_file <- file.path(data_dir, "density_down.txt")
   
   write.table(density_up[, c("chr", "start", "end", "value")],
               density_up_file,
               quote = FALSE, sep = "\t",
               row.names = FALSE, col.names = FALSE)
   
   write.table(density_down[, c("chr", "start", "end", "value")],
               density_down_file,
               quote = FALSE, sep = "\t",
               row.names = FALSE, col.names = FALSE)
   
   message("Written Circos density files: ",
           density_up_file, " and ", density_down_file)
   
   # Annotated density files (with gene_ids)
   density_up_annot_file   <- file.path(data_dir, "density_up_annotated.tsv")
   density_down_annot_file <- file.path(data_dir, "density_down_annotated.tsv")
   
   write.table(density_up,
               density_up_annot_file,
               quote = FALSE, sep = "\t",
               row.names = FALSE, col.names = TRUE)
   
   write.table(density_down,
               density_down_annot_file,
               quote = FALSE, sep = "\t",
               row.names = FALSE, col.names = TRUE)
   
   message("Written annotated density files: ",
           density_up_annot_file, " and ", density_down_annot_file)
   
   ###############################
   # 8. Export annotated tables to Excel
   ###############################
   
   excel_file <- file.path(
   out_dir,
   paste0("circos_annotations_", comparison_label, ".xlsx")
   )
   
   wb <- createWorkbook()
   
   addWorksheet(wb, "scatter_points")
   writeData(wb, "scatter_points", scatter_annot)
   
   addWorksheet(wb, "density_up")
   writeData(wb, "density_up", density_up)
   
   addWorksheet(wb, "density_down")
   writeData(wb, "density_down", density_down)
   
   saveWorkbook(wb, excel_file, overwrite = TRUE)
   
   message("Excel workbook written to: ", excel_file)
   
   message("Done.")

1. Installation and Database Setup

To install and prepare NCBI AMRFinderPlus in the bacto environment:

mamba activate bacto
mamba install ncbi-amrfinderplus
mamba update ncbi-amrfinderplus

mamba activate bacto
amrfinder -u

• This will:
  • Download and install the latest AMRFinderPlus version and its database.
  • Create /home/jhuang/mambaforge/envs/bacto/share/amrfinderplus/data/.
  • Symlink the latest database version for use.

Check available organism options for annotation:

amrfinder --list_organisms
#Available --organism options: Acinetobacter_baumannii, Burkholderia_cepacia, Burkholderia_mallei, Burkholderia_pseudomallei, Campylobacter, Citrobacter_freundii, Clostridioides_difficile, Corynebacterium_diphtheriae, Enterobacter_asburiae, Enterobacter_cloacae, Enterococcus_faecalis, Enterococcus_faecium, Escherichia, Klebsiella_oxytoca, Klebsiella_pneumoniae, Neisseria_gonorrhoeae, Neisseria_meningitidis, Pseudomonas_aeruginosa, Salmonella, Serratia_marcescens, Staphylococcus_aureus, Staphylococcus_pseudintermedius, Streptococcus_agalactiae, Streptococcus_pneumoniae, Streptococcus_pyogenes, Vibrio_cholerae, Vibrio_parahaemolyticus, Vibrio_vulnificus

• Supported values include species such as Escherichia, Klebsiella_pneumoniae, Enterobacter_cloacae, Pseudomonas_aeruginosa and many others.

2. Batch Analysis: Bash Script for Genome Screening

Use the following script to screen multiple genomes using AMRFinderPlus and output only β-lactam/beta-lactamase hits from a metadata table.

Input: genome_metadata.tsv — tab-separated columns: filename_TAB_organism, with header.

filename    organism
58.fasta    Escherichia coli
92.fasta    Klebsiella pneumoniae
125.fasta   Enterobacter cloacae complex
128.fasta   Enterobacter cloacae complex
130.fasta   Enterobacter cloacae complex
147.fasta   Citrobacter freundii
149.fasta   Citrobacter freundii
160.fasta   Citrobacter braakii
161.fasta   Citrobacter braakii
168.fasta   Providencia stuartii
184.fasta   Klebsiella aerogenes
65.fasta    Pseudomonas aeruginosa
201.fasta   Pseudomonas aeruginosa
209.fasta   Pseudomonas aeruginosa
167.fasta   Serratia marcescens

Run:

cd ~/DATA/Data_Patricia_AMRFinderPlus_2025/genomes
./run_amrfinder_and_summarize.sh genome_metadata.tsv
#./run_amrfinder_and_summarize.sh genome_metadata_149.tsv
#OR_DETECT_RUN: amrfinder -n 92.fasta -o amrfinder_results/92.amrfinder.tsv --plus --organism Klebsiella_pneumoniae --threads 1

python summarize_from_amrfinder_results.py amrfinder_results
# or, since that's the default:
# python summarize_from_amrfinder_results.py

Produce

• AMRFinder-wide outputs

  • amrfinder_all.tsv
  • amrfinder_summary_by_isolate_gene.tsv
  • amrfinder_summary_by_gene.tsv
  • amrfinder_summary_by_class.tsv (if a class column exists)
  • amrfinder_summary.xlsx (with multiple sheets)

• β-lactam-only outputs (if Class and Subclass are present)

  • beta_lactam_all.tsv
  • beta_lactam_summary_by_gene.tsv
  • beta_lactam_summary_by_isolate_gene.tsv
  • beta_lactam_all.xlsx
  • beta_lactam_summary.xlsx

Report

Please find attached the updated AMRFinderPlus summary files, now including isolate 167. For β-lactam–specific results, please see beta_lactam_all.xlsx and beta_lactam_summary.xlsx. In particular, beta_lactam_summary.xlsx contains two sheets:

• by_gene – aggregated counts and isolate lists for each β-lactam gene
• by_isolate_gene – per-isolate overview of detected β-lactam genes

Script:

• run_amrfinder_and_summarize.sh

      #!/usr/bin/env bash
      set -euo pipefail
  
      META_FILE="${1:-}"
  
      if [[ -z "$META_FILE" || ! -f "$META_FILE" ]]; then
      echo "Usage: $0 genome_metadata.tsv" >&2
      exit 1
      fi
  
      OUTDIR="amrfinder_results"
      mkdir -p "$OUTDIR"
  
      echo ">>> Checking AMRFinder installation..."
      amrfinder -V || { echo "ERROR: amrfinder not working"; exit 1; }
      echo
  
      echo ">>> Running AMRFinderPlus on all genomes listed in $META_FILE"
  
      # --- loop over metadata file ---
      # expected columns: filename
  <TAB>organism
      tail -n +2 "$META_FILE" | while IFS=$'\t' read -r fasta organism; do
      # skip empty lines
      [[ -z "$fasta" ]] && continue
  
      if [[ ! -f "$fasta" ]]; then
      echo "WARN: FASTA file '$fasta' not found, skipping."
      continue
      fi
  
      isolate_id="${fasta%.fasta}"
  
      # map free-text organism to AMRFinder --organism names (optional)
      org_opt=""
      case "$organism" in
      "Escherichia coli")              org_opt="--organism Escherichia" ;;
      "Klebsiella pneumoniae")         org_opt="--organism Klebsiella_pneumoniae" ;;
      "Enterobacter cloacae complex")  org_opt="--organism Enterobacter_cloacae" ;;
      "Citrobacter freundii")          org_opt="--organism Citrobacter_freundii" ;;
      "Citrobacter braakii")           org_opt="--organism Citrobacter_freundii" ;;
      "Pseudomonas aeruginosa")        org_opt="--organism Pseudomonas_aeruginosa" ;;
      "Serratia marcescens")           org_opt="--organism Serratia_marcescens" ;;
      # others (Providencia stuartii, Klebsiella aerogenes)
      # currently have no organism-specific rules in AMRFinder, so we omit --organism
      *)                               org_opt="" ;;
      esac
  
      out_tsv="${OUTDIR}/${isolate_id}.amrfinder.tsv"
  
      echo "  - ${fasta} (${organism}) -> ${out_tsv} ${org_opt}"
      amrfinder -n "$fasta" -o "$out_tsv" --plus $org_opt
      done
  
      echo ">>> AMRFinderPlus runs finished."
      echo ">>> All done."
      echo "   - Individual reports: ${OUTDIR}/*.amrfinder.tsv"

• summarize_from_amrfinder_results.py

      #!/usr/bin/env python3
      """
      summarize_from_amrfinder_results.py
  
      Usage:
      python summarize_from_amrfinder_results.py [amrfinder_results_dir]
  
      Default directory is "amrfinder_results" (relative to current working dir).
  
      This script:
      1) Reads all *.amrfinder.tsv in the given directory
      2) Merges them into a combined table
      3) Generates AMRFinder-wide summaries (amrfinder_* files)
      4) Applies a β-lactam filter:
  
              Element type == "AMR" (case-insensitive)
      AND Class or Subclass contains "beta-lactam" (case-insensitive)
  
      and generates β-lactam-only summaries (beta_lactam_* files).
  
      It NEVER re-runs AMRFinder; it only uses existing TSV files.
      """
  
      import sys
      import os
      import glob
      import re
  
      try:
      import pandas as pd
      except ImportError:
      sys.stderr.write(
              "ERROR: pandas is not installed.\n"
              "Install with something like:\n"
              "  mamba install pandas openpyxl -c conda-forge -c bioconda\n"
      )
      sys.exit(1)
  
      # ---------------------------------------------------------------------
      # Helpers
      # ---------------------------------------------------------------------
  
      def read_one(path):
      """Read one *.amrfinder.tsv and add an 'isolate_id' column from the filename."""
      df = pd.read_csv(path, sep="\t", dtype=str)
      df.columns = [c.strip() for c in df.columns]
      isolate_id = os.path.basename(path).replace(".amrfinder.tsv", "")
      df["isolate_id"] = isolate_id
      return df
  
      def pick(df, *candidates):
      """Return the first existing column name among candidates (normalized names)."""
      for c in candidates:
              if c in df.columns:
              return c
      return None
  
      # ---------------------------------------------------------------------
      # AMRFinder-wide summaries (no β-lactam filter)
      # ---------------------------------------------------------------------
  
      def make_amrfinder_summaries(
      df_all,
      col_gene,
      col_seq,
      col_class,
      col_subcls,
      col_ident,
      col_cov,
      col_iso,
      ):
      """Summaries for ALL AMRFinder hits (no β-lactam filter)."""
      if df_all.empty:
              print("[amrfinder] No rows in merged table, skipping summaries.")
              return
  
      # full merged table
      df_all.to_csv("amrfinder_all.tsv", sep="\t", index=False)
      print(">>> Full AMRFinder table written to: amrfinder_all.tsv")
  
      # ---- summary by isolate × gene ----
      rows = []
      for (iso, gene), sub in df_all.groupby([col_iso, col_gene], dropna=False):
              row = {
              "isolate_id": iso,
              "Gene_symbol": sub[col_gene].iloc[0],
              "n_hits": len(sub),
              }
              if col_seq is not None:
              row["Sequence_name"] = sub[col_seq].iloc[0]
              if col_class is not None:
              row["Class"] = sub[col_class].iloc[0]
              if col_subcls is not None:
              row["Subclass"] = sub[col_subcls].iloc[0]
              if col_ident is not None:
              vals = pd.to_numeric(sub[col_ident], errors="coerce")
              row["%identity_min"] = vals.min()
              row["%identity_max"] = vals.max()
              if col_cov is not None:
              vals = pd.to_numeric(sub[col_cov], errors="coerce")
              row["%coverage_min"] = vals.min()
              row["%coverage_max"] = vals.max()
              rows.append(row)
  
      summary_iso_gene = pd.DataFrame(rows)
      summary_iso_gene.to_csv(
              "amrfinder_summary_by_isolate_gene.tsv", sep="\t", index=False
      )
      print(">>> Isolate × gene summary written to: amrfinder_summary_by_isolate_gene.tsv")
  
      # ---- summary by gene ----
      def join(vals):
              uniq = sorted(set(vals.dropna().astype(str)))
              return ",".join(uniq)
  
      rows = []
      for gene, sub in df_all.groupby(col_gene, dropna=False):
              row = {
              "Gene_symbol": sub[col_gene].iloc[0],
              "n_isolates": sub[col_iso].nunique(),
              "isolates": ",".join(sorted(set(sub[col_iso].dropna().astype(str)))),
              "n_hits": len(sub),
              }
              if col_seq is not None:
              row["Sequence_name"] = join(sub[col_seq])
              if col_class is not None:
              row["Class"] = join(sub[col_class])
              if col_subcls is not None:
              row["Subclass"] = join(sub[col_subcls])
              rows.append(row)
  
      summary_gene = pd.DataFrame(rows)
      summary_gene = summary_gene.sort_values("n_isolates", ascending=False)
      summary_gene.to_csv("amrfinder_summary_by_gene.tsv", sep="\t", index=False)
      print(">>> Gene-level summary written to: amrfinder_summary_by_gene.tsv")
  
      # ---- summary by class/subclass ----
      summary_class = None
      if col_class is not None:
              group_cols = [col_class]
              if col_subcls is not None:
              group_cols.append(col_subcls)
  
              summary_class = (
              df_all.groupby(group_cols, dropna=False)
              .agg(
                      n_isolates=(col_iso, "nunique"),
                      n_hits=(col_iso, "size"),
              )
              .reset_index()
              )
              summary_class.to_csv("amrfinder_summary_by_class.tsv", sep="\t", index=False)
              print(">>> Class-level summary written to: amrfinder_summary_by_class.tsv")
      else:
              print(">>> No 'class' column found; amrfinder_summary_by_class.tsv not created.")
  
      # ---- Excel workbook ----
      try:
              with pd.ExcelWriter("amrfinder_summary.xlsx") as xw:
              df_all.to_excel(xw, sheet_name="amrfinder_all", index=False)
              summary_iso_gene.to_excel(xw, sheet_name="by_isolate_gene", index=False)
              summary_gene.to_excel(xw, sheet_name="by_gene", index=False)
              if summary_class is not None:
                      summary_class.to_excel(xw, sheet_name="by_class", index=False)
              print(">>> Excel workbook written: amrfinder_summary.xlsx")
      except Exception as e:
              print("WARNING: could not write amrfinder_summary.xlsx:", e)
  
      # ---------------------------------------------------------------------
      # β-lactam summaries
      # ---------------------------------------------------------------------
  
      def make_beta_lactam_summaries(
      df_beta,
      col_gene,
      col_seq,
      col_subcls,
      col_ident,
      col_cov,
      col_iso,
      ):
      """Summaries for β-lactam subset (after mask)."""
      if df_beta.empty:
              print("[beta_lactam] No β-lactam hits in subset, skipping.")
              return
  
      # full β-lactam table
      beta_all_tsv = "beta_lactam_all.tsv"
      df_beta.to_csv(beta_all_tsv, sep="\t", index=False)
      print(">>> β-lactam / β-lactamase hits written to: %s" % beta_all_tsv)
  
      # -------- summary by gene (with list of isolates) --------
      group_cols = [col_gene]
      if col_seq is not None:
              group_cols.append(col_seq)
      if col_subcls is not None:
              group_cols.append(col_subcls)
  
      def join_isolates(vals):
              uniq = sorted(set(vals.dropna().astype(str)))
              return ",".join(uniq)
  
      summary_gene = (
              df_beta.groupby(group_cols, dropna=False)
              .agg(
              n_isolates=(col_iso, "nunique"),
              isolates=(col_iso, join_isolates),
              n_hits=(col_iso, "size"),
              )
              .reset_index()
      )
  
      rename_map = {}
      if col_gene is not None:
              rename_map[col_gene] = "Gene_symbol"
      if col_seq is not None:
              rename_map[col_seq] = "Sequence_name"
      if col_subcls is not None:
              rename_map[col_subcls] = "Subclass"
      summary_gene.rename(columns=rename_map, inplace=True)
  
      sum_gene_tsv = "beta_lactam_summary_by_gene.tsv"
      summary_gene.to_csv(sum_gene_tsv, sep="\t", index=False)
      print(">>> Gene-level β-lactam summary written to: %s" % sum_gene_tsv)
      print("    (includes 'isolates' = comma-separated isolate_ids)")
  
      # -------- summary by isolate & gene (with annotation) --------
      rows = []
      for (iso, gene), sub in df_beta.groupby([col_iso, col_gene], dropna=False):
              row = {
              "isolate_id": iso,
              "Gene_symbol": sub[col_gene].iloc[0],
              "n_hits": len(sub),
              }
              if col_seq is not None:
              row["Sequence_name"] = sub[col_seq].iloc[0]
              if col_subcls is not None:
              row["Subclass"] = sub[col_subcls].iloc[0]
  
              if col_ident is not None:
              vals = pd.to_numeric(sub[col_ident], errors="coerce")
              row["%identity_min"] = vals.min()
              row["%identity_max"] = vals.max()
              if col_cov is not None:
              vals = pd.to_numeric(sub[col_cov], errors="coerce")
              row["%coverage_min"] = vals.min()
              row["%coverage_max"] = vals.max()
  
              rows.append(row)
  
      summary_iso_gene = pd.DataFrame(rows)
      sum_iso_gene_tsv = "beta_lactam_summary_by_isolate_gene.tsv"
      summary_iso_gene.to_csv(sum_iso_gene_tsv, sep="\t", index=False)
      print(">>> Isolate × gene β-lactam summary written to: %s" % sum_iso_gene_tsv)
      print("    (includes 'Gene_symbol' and 'Sequence_name' annotation columns)")
  
      # -------- optional Excel exports --------
      try:
              with pd.ExcelWriter("beta_lactam_all.xlsx") as xw:
              df_beta.to_excel(xw, sheet_name="beta_lactam_all", index=False)
              with pd.ExcelWriter("beta_lactam_summary.xlsx") as xw:
              summary_gene.to_excel(xw, sheet_name="by_gene", index=False)
              summary_iso_gene.to_excel(xw, sheet_name="by_isolate_gene", index=False)
              print(">>> Excel workbooks written: beta_lactam_all.xlsx, beta_lactam_summary.xlsx")
      except Exception as e:
              print("WARNING: could not write β-lactam Excel files:", e)
  
      # ---------------------------------------------------------------------
      # Main
      # ---------------------------------------------------------------------
  
      def main():
      outdir = sys.argv[1] if len(sys.argv) > 1 else "amrfinder_results"
  
      if not os.path.isdir(outdir):
              sys.stderr.write("ERROR: directory '%s' not found.\n" % outdir)
              sys.exit(1)
  
      files = sorted(glob.glob(os.path.join(outdir, "*.amrfinder.tsv")))
      if not files:
              sys.stderr.write("ERROR: no *.amrfinder.tsv files found in '%s'.\n" % outdir)
              sys.exit(1)
  
      print(">>> Found %d AMRFinder TSV files in: %s" % (len(files), outdir))
      for f in files:
              print("   -", os.path.basename(f))
  
      dfs = [read_one(f) for f in files]
      df = pd.concat(dfs, ignore_index=True)
  
      # normalize column names for internal use
      norm_cols = {c: c.strip().lower().replace(" ", "_") for c in df.columns}
      df.rename(columns=norm_cols, inplace=True)
  
      # locate columns (handles your Element type / subtype + older formats)
      col_gene       = pick(df, "gene_symbol", "genesymbol")
      col_seq        = pick(df, "sequence_name", "sequencename")
      col_elemtype   = pick(df, "element_type")
      col_elemsub    = pick(df, "element_subtype")
      col_class      = pick(df, "class")
      col_subcls     = pick(df, "subclass")
      col_ident      = pick(df, "%identity_to_reference_sequence", "identity")
      col_cov        = pick(df, "%coverage_of_reference_sequence", "coverage_of_reference_sequence")
      col_iso        = "isolate_id"
  
      print("\nDetected columns:")
      for label, col in [
              ("gene", col_gene),
              ("sequence", col_seq),
              ("element_type", col_elemtype),
              ("element_subtype", col_elemsub),
              ("class", col_class),
              ("subclass", col_subcls),
              ("%identity", col_ident),
              ("%coverage", col_cov),
              ("isolate_id", col_iso),
      ]:
              print("  %-14s: %s" % (label, col))
  
      if col_gene is None:
              sys.stderr.write(
              "ERROR: could not find a gene symbol column "
              "(expected something like 'Gene symbol' in the original AMRFinder output).\n"
              )
              sys.exit(1)
  
      print("\n=== Generating AMRFinder-wide summaries (all hits) ===")
      make_amrfinder_summaries(
              df_all=df,
              col_gene=col_gene,
              col_seq=col_seq,
              col_class=col_class,
              col_subcls=col_subcls,
              col_ident=col_ident,
              col_cov=col_cov,
              col_iso=col_iso,
      )
  
      # -----------------------------------------------------------------
      # β-lactam subset
      #
      # New logic for your current data:
      #   Element type == "AMR"
      #   AND Class or Subclass contains "beta-lactam"
      #
      # Falls back to just Class/Subclass if Element type not present.
      # -----------------------------------------------------------------
      if (col_elemtype is not None) or (col_class is not None or col_subcls is not None):
  
              # element type AMR (if column exists, otherwise all True)
              if col_elemtype is not None:
              mask_amr = df[col_elemtype].str.contains("AMR", case=False, na=False)
              else:
              mask_amr = pd.Series(True, index=df.index)
  
              # beta-lactam pattern (handles BETA-LACTAM, beta lactam, etc.)
              beta_pattern = re.compile(r"beta[- ]?lactam", re.IGNORECASE)
  
              mask_beta = pd.Series(False, index=df.index)
              if col_class is not None:
              mask_beta |= df[col_class].fillna("").str.contains(beta_pattern)
              if col_subcls is not None:
              mask_beta |= df[col_subcls].fillna("").str.contains(beta_pattern)
  
              mask = mask_amr & mask_beta
              df_beta = df.loc[mask].copy()
  
              if df_beta.empty:
              print(
                      "\nWARNING: No β-lactam hits found "
                      "(Element type == 'AMR' AND Class/Subclass contains 'beta-lactam')."
              )
              else:
              print(
                      "\n=== β-lactam subset ===\n"
                      "  kept %d of %d rows where Element type is 'AMR' and "
                      "Class/Subclass contains 'beta-lactam'\n"
                      % (len(df_beta), len(df))
              )
              make_beta_lactam_summaries(
                      df_beta=df_beta,
                      col_gene=col_gene,
                      col_seq=col_seq,
                      col_subcls=col_subcls,
                      col_ident=col_ident,
                      col_cov=col_cov,
                      col_iso=col_iso,
              )
      else:
              print(
              "\nWARNING: Cannot apply β-lactam filter because Element type and/or "
              "class/subclass columns were not found. Only amrfinder_* "
              "outputs were generated."
              )
  
      if __name__ == "__main__":
      main()

Title: Automated Kymograph Track Filtering & Lake File Generation (kymograph轨迹自动过滤与Lake文件生成流程)

Step 1 – Track Filtering with 1_filter_track.py

（用1_filter_track.py进行轨迹过滤） 运行命令：

python 1_filter_track.py  

核心思路：对每个原始*_blue.csv轨迹文件，根据位置和寿命（lifetime）进行过滤，将保留的轨迹和被剔除的轨迹分别存放于两个目录：

• filtered/ → 通过过滤条件保留下来的轨迹
• separated/ → 不满足过滤条件被剔除的轨迹 共有74个原始*_blue.csv文件。 确保对每个原始blue文件，针对每种过滤条件输出对应文件：
• 有轨迹通过过滤时，生成正常的filtered CSV（含数据行）
• 无轨迹通过过滤时，生成占位placeholder文件，格式正确，仅含header和注释，无数据 此设计确保后续2_update_lakes.py能正常读取，并判定该条件下无轨迹，保证流水线完整一致。

Step 2 – Organize filtered CSVs and Fix p940 Naming Bug

（整理过滤结果CSV，修正文件名命名错误） 创建文件夹：

mkdir filtered_blue_position filtered_blue_position_1s filtered_blue_position_5s filtered_blue_lifetime_5s_only

移动对应过滤文件：

1) 绑定位置2.2–3.8 µm

mv filtered/*_blue_position.csv filtered_blue_position

2) 绑定位置且寿命≥1s

mv filtered/*_blue_position_1s.csv filtered_blue_position_1s

3) 绑定位置且寿命≥5s

mv filtered/*_blue_position_5s.csv filtered_blue_position_5s

4) 寿命≥5s不限制位置

mv filtered/*_blue_lifetime_5s_only.csv filtered_blue_lifetime_5s_only

每个目录保留74个CSV文件（包含真实轨迹和header-only占位符）。 修正p940命名bug（文件名中p940与lake文件中940不匹配），统一去除多余的p：

find filtered_blue_position -type f -name 'p*_p[0-9][0-9][0-9]_*.csv' -exec rename 's/_p([0-9]{3})/_$1/' {} +
（同理在其它三个目录执行相同命令）

保证轨迹CSV名与lake文件中kymograph名称一一对应。

Step 3 – Write filtered tracks back to lake files

（把过滤后轨迹写回lake文件） 运行命令更新lake文件（每组过滤条件对应一组输出目录）：

python 2_update_lakes.py --merged_lake_folder lakes_raw --filtered_folder filtered_blue_position --output_folder lakes_blue_position_2.2-3.8 | tee blue_position_2.2-3.8.LOG

python 2_update_lakes.py --merged_lake_folder lakes_raw --filtered_folder filtered_blue_position_1s --output_folder lakes_blue_position_2.2-3.8_length_min_1s | tee blue_position_2.2-3.8_length_min_1s.LOG

python 2_update_lakes.py --merged_lake_folder lakes_raw --filtered_folder filtered_blue_position_5s --output_folder lakes_blue_position_2.2-3.8_length_min_5s | tee blue_position_2.2-3.8_length_min_5s.LOG

python 2_update_lakes.py --merged_lake_folder lakes_raw --filtered_folder filtered_blue_lifetime_5s_only --output_folder lakes_blue_length_min_5s | tee blue_length_min_5s.LOG

处理逻辑：

• 通过kymograph名称匹配filtered_*目录对应CSV
• 根据CSV内容重建blue轨迹文本，写回lake JSON
• 分类日志输出三种情况：

1. Updated：找到CSV且≥1条轨迹，更新保存轨迹
2. CSV存在但无轨迹或读取失败，移除kymograph及关联H5链接
3. 无匹配CSV，移除kymograph及H5链接
   • 日志统计统计各case数量、CSV总数、未使用“孤儿”CSV

最终实现每个replicate拥有多组更新的lake文件，各文件中kymographs、experiments[…].dataset、file_viewer的H5链接一致对应，确保完整性和可追踪性。

此流程自动化实现kymograph轨迹质量控制与lake文件二次生成，支持多样过滤条件，保证下游数据分析准确可靠。

FAU“身体活动与健康”硕士项目：申请指南与入学要求

如何申请 (How to apply)

“身体活动与健康”硕士项目（MA Programme Physical Activity and Health）只能在冬季学期开始（课程于2024年10月开课），针对2025/26冬季学期的申请将于2025年2月15日开始。申请截止日期为2025年5月31日。我们建议非欧盟公民最迟于2025年3月31日前提交申请，以便有充足时间办理签证手续。 所有所需申请材料必须通过线上系统[Campo(https://www.campo.fau.de/qisserver/pages/cs/sys/portal/hisinoneStartPage.faces)提交。（请不要邮寄任何申请材料到FAU，所有文件须通过Campo平台在线上传。）

所需材料 (Required Documents)

在申请“身体活动与健康”硕士项目时，需要提交以下文件：

• 个人简历
• 动机信（1–2页），说明你申请该项目的兴趣、动机与资质
• 德国高校毕业生：提交所有教育阶段的毕业证书及成绩单（如成绩单、Studienbuch等）复印件。
• 国际高校毕业生：提交所有教育阶段的经认证的毕业证书及成绩单复印件。
• 若你的学位为体育教育、心理学、社会学、政治学、人类学或医学：请提交与你本项目高度相关课程的清单，并附上至少一年全职的相关领域（运动科学/康复科学/治疗科学/公共卫生）的工作经验证明。
• 对于母语非英语且本科/硕士授课语言不是英语者：至少需提供B2级别英语能力证明。
• 对于母语非德语者（如有）：至少A1级别的德语语言证书。

动机信 (Cover letter)

动机信是你申请材料的重要部分。请说明为什么想加入本项目，以及你未来的职业规划。此外，应提及你先前在身体活动、物理治疗或公共卫生等主题领域的经验。篇幅应为1至2页。

个人简历 (Curriculum vitae/Resume)

简历应简要说明你的中学和大学学习经历，列出最近就读的所有学校或大学。包括与你申请项目相关的实习、兼职或全职工作经历。同时应注明出生日期与地点、国籍及现居地点。可以使用Europass简历模板（下载模板及说明，或访问Europass主页）。

经认证的复印件（仅限国际学位申请者）Certified copies (applicants with international degrees only)

需提交中学和大学期间所有学历及成绩单的经认证复印件。这些文件仅通过电子邮件提交（不接受邮寄或传真）。所有复印件必须经过正式认证。认证文件须：

• 含有认证机构印章；
• 由认证人员签字；
• 明确标注认证日期；
• 认证机构及人员具备认证资格。 通常，学校管理部门有权办理认证；如不确定，可咨询就近的德国大使馆或领事馆。

课程清单（适用于体育教育、心理学、社会学、政治学、人类学、医学等学位申请者）

项目对非运动科学、物理治疗、康复科学、健康教育等背景的学生开放。请列出与你本项目相关的所有课程，例如运动科学、体育教育、物理治疗、康复科学、老年学、公共卫生、流行病学、研究方法或统计学等。

Listing of courses/classes with high relevance to our programme (for applicants with degrees in Physical Education, Psychology, Sociology, Political Science, Anthropology, or Medicine only) The programme is open to students who do not have degrees in Sport Science, Kinesiology /Exercise Science, Physiotherapy, Rehabilitation Science, Health Education, Health Science/Public Health.

Such other degrees can be e.g. Physical Education, Psychology, Sociology, Political Science, Anthropology, or Medicine). This list should provide us with a brief summary of all classes or coursework that you have attended and that are relevant to the subject areas of physical activity and/or (public) health. Potential examples include courses/classes covering the topics of sport science, physical education, physical therapy, rehabilitation science, kinesiology, gerontology, public health, epidemiology, research methods, or statistics.

工作经验（适用于体育教育、心理学、社会学、政治学、人类学、医学等学位申请者）

具有上述专业背景的学生，需提供至少一年全职相关工作经验（运动科学、康复科学、治疗科学或公共卫生领域）证明，可由相关机构出具证明信。

Documentation of 1 year work experience in the fields of Sport Science/ Rehabilitation Science or Therapeutic Science/ Public Health (for applicants with degrees in Physical Education, Psychology, Sociology, Political Science, Anthropology, or Medicine only) The programme is open to students who do not have degrees in Sport Science, Kinesiology /Exercise Science, Physiotherapy, Rehabilitation Science, Therapeutic Science, Health Education, or Health Science/Public Health.

Such other degrees can be e.g. Physical Education, Psychology, Sociology, Political Science, Anthropology, or Medicine. Students with such degrees need to document at least 1 year of work experience (full-time) in the fields of Sport Science, Reahbilitation Science or Therapeutic Science, or Public Health in order to be eligible to apply to the programme. The documentation can be an attached letter from the institution/company.

英语语言证书（仅适用于母语非英语者）

本项目以英语授课，需要具备足够的听、说、读、写能力。若母语非英语且本科/硕士授课语言非英语，需提供语言证书证明达到我们要求的水平。最低要求为CEFR体系的B2级。详情见入学要求。

德语语言证书（仅适用于母语非德语者）

根据州级规定，所有母语非德语学生须在入学一年内至少达到A1级德语水平。若已有德语水平，请在申请材料中提供证明；若尚未具备，也可申请。大学提供免费德语课程，可在第一学年内学习。所有课程与考试均以英语进行。

申请提交地点

所有申请材料须通过Campo平台提交。（请不要邮寄任何文件至FAU，所有文件仅通过Campo上传。）

申请审核

DSS部门两位教师将依据以下标准评审申请：

• 运动科学、康复科学/治疗科学及公共卫生方面的先前知识；
• 相关领域（体育教育、心理学、社会学、政治学、人类学或医学）的知识背景；
• 研究方法（如统计学、质性研究）的知识；
• 在运动科学、康复科学/治疗科学及公共卫生领域的实践经验（如实习或工作经验）。 由于通常申请数量众多，请预留至少4周等待评审结果。

补充信息

如有关于硕士项目内容或申请流程的疑问，请联系项目顾问Karim Abu-Omar。 若你已通过Campo提交申请，请在联系时务必提供申请编号（application-ID），并在需要通过邮件发送的文件名称中注明该编号。]

“身体活动与健康”硕士项目的申请资格需通过以下条件证明：

具有以下学科之一的高等教育第一阶段学位（例如学士学位，或德国体系中的“Diplom”或“Staatsexamen”）：

• 运动科学（以健康为重点）
• 运动机能学/运动科学（以健康为重点）
• 康复科学/治疗科学
• 健康教育
• 健康科学/公共卫生

在特殊情况下，若申请人完成了以下相关领域的类似学位，也可被录取，例如体育教育、心理学、社会学、政治学、人类学或医学。申请人需提供证明，证明其已在运动科学/康复科学/治疗科学/公共卫生等领域修读了至少20个ECTS学分的课程，或在这些领域拥有至少1年的全职工作经验。

最低成绩要求：

• 对于采用百分制的评分体系：总成绩须达到75%或以上。
• 对于采用4分制GPA体系（如美国）：GPA须达到3.00或以上。
• 对于德国学生：成绩须达到2.5或以下。

目前仍在读本科的学生，在修完至少140个ECTS学分后即可申请。 正式录取前，必须提交最终成绩单及学士学位证书；被录取的申请者若尚未提交最终文件，其录取为有条件录取。

语言要求：英语

本硕士项目的所有课程均以英语授课。 所有母语非英语的申请者，须提供至少达到CEFR欧洲语言能力等级框架B2级的英语语言能力证明。

若你持有其他类型的语言证书，可参考以下证书等级对照表，以了解与CEFR B2等级约等的分数范围： 语言证书对照表 请注意：该对照表仅用于参考，不具法律效力。若提交的证书或成绩未标明CEFR等级，将由大学逐一评估是否符合要求。 未能提供CEFR B2水平英语证明的申请人，可能需要在入学前于大学语言中心参加英语水平测试。

语言要求：德语

入学时无须提供德语能力证明，但学生须在赴FAU就读的第一学年内学习德语，至少达到A1级。 建议申请者具备基本德语能力，特别是第二学年以项目研究为主的课程阶段。 大学语言中心为所有语言水平的学生提供免费的德语课程。

学费

本硕士项目不收取学费。