Summary
这里介绍的方案描述了一个完整的管道，用于分析从原始读数到功能分析的RNA测序转录组数据，包括质量控制和预处理步骤到高级统计分析方法。
Abstract
病原体可引起多种传染病。宿主对感染的反应诱导的生物过程决定了疾病的严重程度。为了研究这些过程，研究人员可以使用高通量测序技术（RNA-seq），以测量宿主转录组在感染的不同阶段，临床结果或疾病严重程度的动态变化。这项调查可以更好地了解疾病，并发现潜在的药物靶点和治疗方法。这里介绍的方案描述了一个完整的管道，用于分析从原始读取到功能分析的RNA测序数据。管道分为五个步骤：（1）数据的质量控制;（2）基因的作图和注释;（3）统计分析，鉴定差异表达基因和共表达基因;（4）测定样品扰动的分子程度;（5）功能分析。步骤 1 删除了可能影响下游分析质量的技术工件。在第2步中，根据标准文库协议绘制基因图谱并进行注释。步骤3中的统计分析可识别感染样本中差异表达或共表达的基因，与未感染的基因进行比较。在步骤4中使用分子扰动度方法验证样品变异性和潜在生物异常值的存在。最后，步骤5中的功能分析揭示了与疾病表型相关的途径。所提出的管道旨在通过来自宿主 - 病原体相互作用研究的RNA-seq数据分析来支持研究人员，并推动未来的体外 或 体内 实验，这对于了解感染的分子机制至关重要。
Introduction
登革热、黄热病、基孔肯雅热和寨卡等虫媒病毒与几次地方性疫情广泛相关，并已成为过去几十年中导致人类感染的主要病原体之一1，2。感染基孔肯雅病毒（CHIKV）的个体经常出现发热、头痛、皮疹、多关节痛和关节炎3，4，5。病毒可以破坏细胞的基因表达并影响各种宿主信号通路。最近，血液转录组研究利用RNA-seq鉴定与急性CHIKV感染相关的差异表达基因（DEGs）与康复6 或健康对照组7进行比较。CHIKV感染的儿童具有上调的基因，这些基因参与先天免疫，例如与病毒RNA，JAK / STAT信号传导和Toll样受体信号通路的细胞传感器相关的基因6。急性感染CHIKV的成人还显示出与先天免疫相关的基因的诱导，例如与单核细胞和树突状细胞活化相关的基因，以及与抗病毒反应相关的基因7。富含下调基因的信号通路包括与适应性免疫相关的信号通路，例如T细胞活化以及T细胞和B细胞的分化和富集7。
可以使用几种方法来分析宿主和病原体基因的转录组数据。通常，RNA-seq文库的制备从成熟poly-A转录本的富集开始。该步骤去除大部分核糖体RNA（rRNA），在某些情况下去除病毒/细菌RNA。然而，当生物学问题涉及病原体转录本检测并且RNA独立于先前的选择进行测序时，可以通过测序检测许多其他不同的转录本。例如，亚基因组mRNA已被证明是验证疾病严重程度的重要因素8。此外，对于某些病毒，如CHIKV和SARS-CoV-2，即使是富含poly-A的文库也会生成病毒读数，可用于下游分析9，10。当专注于宿主转录组的分析时，研究人员可以研究样品之间的生物扰动，鉴定差异表达的基因和富集的途径，并生成共表达模块7，11，12。该协议突出了使用不同生物信息学方法对CHIKV感染患者和健康个体的转录组分析（图1A）。来自先前发表的一项研究7 的数据包括20名健康和39名CHIKV急性感染个体，用于产生具有代表性的结果。
Subscription Required. Please recommend JoVE to your librarian.
Protocol
该协议中使用的样品由圣保罗大学生物医学科学研究所微生物学系和塞尔希培联邦大学伦理委员会批准（协议分别为54937216.5.0000.5467和54835916.2.0000.5546）。
1. Docker 桌面安装
注意：准备 Docker 环境的步骤因操作系统 （OS） 而异。因此，Mac 用户必须遵循列为 1.1 的步骤，Linux 用户必须遵循列为 1.2 的步骤，Windows 用户必须遵循列为 1.3 的步骤。
- 在 MacOS 上安装。
- 访问 Get Docker 网站（材料表），单击 Docker Desktop for Mac ，然后单击" 从 Docker Hub 下载 "链接。
- 通过单击" 获取 Docker "按钮下载安装文件。
- 执行 Docker.dmg文件以打开安装程序，然后将图标拖到 "应用程序" 文件夹中。本地化并执行" 应用程序" 文件夹中的 Docker.app 以启动程序。
注：顶部状态栏中的软件特定菜单指示软件正在运行，并且可从终端访问该软件。
- 在 Linux 操作系统上安装容器程序。
- 访问 Get Docker Linux 网站（材料表），并按照使用 Docker Linux 存储库链接上提供的存储库部分进行安装的说明进行操作。
- 使用命令行更新所有 Linux 软件包：
sudo apt-get update - 将所需的软件包安装到 Docker：
sudo apt-get install apt-transport-https ca-certificates curl gnupg lsb-release - 创建软件归档密钥环文件：
curl -fsSL https://download.docker.com/linux/ubuntu/gpg |sudo gpg --dearmor -o /usr/share/keyrings/docker-archive-keyring.gpg - 在 source.list 文件中添加 Docker deb 信息：
echo "deb [arch=amd64 signed-by=/usr/share/keyrings/docker-archive-keyring.gpg] https://download.docker.com/linux/ubuntu $（lsb_release -cs） stable" |sudo tee /etc/apt/sources.list.d/docker.list > /dev/null - 再次更新所有软件包，包括最近添加的软件包：
sudo apt-get update - 安装桌面版本：
sudo apt-get install docker-ce docker-ce-cli containerd.io - 选择地理区域和时区以完成安装过程。
- 在 Windows 操作系统上安装容器程序。
- 访问 Get Docker 网站（材料表），然后单击" 入门"。查找 Docker Desktop for Windows 的安装程序。下载文件并将其本地安装在计算机上。
- 下载后，启动安装文件（.exe）并保留默认参数。确保标记了"安装 WSL 2 所需的 Windows 组件"和"将快捷方式添加到桌面"这两个选项。
注意：在某些情况下，当此软件尝试启动服务时，它会显示错误：WSL 安装不完整。要找出此错误，请访问网站 WSL2-Kernel（材料表）。 - 下载并安装最新的 WSL2 Linux 内核。
- 以管理员身份访问 PowerShell 终端并执行以下命令：
dism.exe /online /enable-feature /featurename：Microsoft-Windows-Subsystem-Linux /all /norestart - 确保软件 Docker 桌面已成功安装。
- 从 Docker 中心（材料表）上的 CSBL 存储库下载映像。
- 打开 Docker 桌面，并验证工具栏左下角的状态是否为"正在运行"。
- 转到 Windows PowerShell 终端命令行。从 Docker 中心上的 CSBL 存储库下载此协议的 Linux 容器映像。执行以下命令下载镜像：
Docker pull csblusp/transcriptome
注意：下载映像后，可以在 Docker 桌面中看到该文件。若要创建容器，Windows 用户必须遵循步骤 1.5，而 Linux 用户必须遵循步骤 1.6。
- 初始化 Windows 操作系统上的服务器容器。
- 从工具栏在桌面应用管理器中查看 Docker 映像文件，并访问"映像"页。
注意：如果管道映像已成功下载，则将有一个 csblusp/transcriptome 映像可用。 - 通过单击" 运行 "按钮从 csblusp/transcriptome 图像启动容器。展开 "可选设置" 以配置容器。
- 定义容器名称（例如，服务器）。
- 将本地计算机中的文件夹与 Docker 中的文件夹相关联。为此，请确定主机路径。在本地计算机中设置一个文件夹，以存储将在结束时下载的已处理数据。设置容器路径。定义 csblusp/transcriptome 容器文件夹并将其链接到本地计算机路径（使用名称"/opt/transferdata"作为容器路径）。
- 之后，单击" 运行 "以创建 csblusp/transcriptome 容器。
- 要从 csblusp/transcriptome 容器访问 Linux 终端，请单击 CLI 按钮。
- 键入 bash 终端以获得更好的体验。为此，请执行以下命令：
砰 - 执行 bash 命令后，确保终端显示 （root@
：/#）：
root@ac12c583b731：/#
- 从工具栏在桌面应用管理器中查看 Docker 映像文件，并访问"映像"页。
- 初始化 Linux 操作系统的服务器容器。
- 执行以下命令以基于映像创建 Docker 容器：
docker run -d -it --rm --name server -v：/opt/transferdata csblusp/transcriptome
注意：<主机路径>：定义本地文件夹计算机的路径。 - 执行以下命令以访问 Docker 容器的命令终端：
Docker exec -it server bash - 确保 Linux 终端的可用性，以便使用命令行执行任何程序/脚本。
- 执行 bash 命令后，确保终端显示 （root@
：/#）：
root@ac12c583b731：/#
注意：默认情况下，root 密码为"转录组"。如果需要，可以通过执行以下命令来更改 root 密码：
路卫 - 首先，执行 source 命令以 addpath.sh 以确保所有工具都可用。执行命令：
源 /选择/添加路径.sh
- 执行以下命令以基于映像创建 Docker 容器：
- 检查RNA测序文件夹的结构。
- 访问转录组管道脚本文件夹，并确保来自RNA测序的所有数据都存储在该文件夹中：/home/transcriptome-pipeline/data。
- 确保从分析中获得的所有结果都存储在路径/home/transcriptome-pipeline/results的文件夹中。
- 确保基因组和注释参考文件存储在路径/home/transcriptome-pipeline/datasets的文件夹中。这些文件将有助于支持所有分析。
- 确保所有脚本都存储在路径/home/transcriptome-pipeline/scripts的文件夹中，并按每个步骤分隔，如下所述。
- 下载注释和人类基因组。
- 访问脚本文件夹：
cd /home/transcriptome-pipeline/scripts - 执行以下命令以下载参考人类基因组：
巴什·downloadGenome.sh - 要下载注释，请执行以下命令：
bash downloadAnnotation.sh
- 访问脚本文件夹：
- 更改参考基因组的注释或版本。
- 打开 downloadAnnotation.sh 并 downloadGenome.sh 以更改每个文件的 URL。
- 将 downloadAnnotation.sh 和 downloadGenome.sh 文件复制到传输区域，然后在本地操作系统中进行编辑。
cd /home/transcriptome-pipeline/scripts
cp downloadAnnotation.sh downloadGenome.sh /opt/transferdata - 打开" 主机路径" 文件夹，该文件夹在步骤 1.5.4 中选择用于在主机和 Docker 容器之间进行链接。
- 使用首选编辑器软件编辑文件并保存。最后，将修改后的文件放入脚本文件夹中。执行命令：
cd /opt/transferdata
cp downloadAnnotation.sh downloadGenome.sh /home/transcriptome-pipeline/scripts
注意：这些文件可以使用 vim 或 nano Linux 编辑器直接编辑。
- 接下来，使用命令行配置 fastq-dump 工具：
vdb-config --interactive
注意：这允许从示例数据下载排序文件。- 使用 Tab 键导航 "工具" 页面，然后选择当前文件夹选项。导航到 "保存 "选项，然后单击" 确定"。然后， 退出 快速转储工具。
- 开始下载以前发表的论文中的读数7。每个样本的 SRA 加入编号是必需的。从SRA NCBI网站（材料表）获取SRA编号。
注意：要分析公共数据库中可用的RNA-Seq数据，请遵循步骤1.12。要分析私有RNA-seq数据，请遵循步骤1.13。 - 分析特定的公共数据。
- 访问国家生物技术信息中心（NCBI）网站并查找特定主题的关键字。
- 单击"基因组"部分中的"生物项目的结果"链接。
- 选择并单击特定研究。单击 "SRA 实验"。将打开一个新页面，其中显示了可用于此研究的所有样本。
- 单击加入号上方的 "发送到：" 。在 "选择目标" 选项中，选择" 文件 和 格式" 选项，然后选择" RunInfo"。单击 "创建文件" 以导出所有库信息。
- 将 SraRunInfo.csv 文件保存在 1.5.4 步骤中定义的主机路径中，然后执行下载脚本：
cp /opt/transferdata/SraRunInfo.csv /home/transcriptome-pipeline/data
cd /home/transcriptome-pipeline/scripts
bash downloadAllLibraries.sh
- 分析私有和未发布的测序数据。
- 将排序数据组织到名为 Reads 的文件夹中。
注意：在 "读取" 文件夹中，为每个示例创建一个文件夹。对于每个示例，这些文件夹必须具有相同的名称。将每个示例的数据添加到其目录中。如果它是一个配对端的RNA-Seq，每个样本目录应该包含两个FASTQ文件，这些文件必须分别显示根据模式{sample}_1.fastq.gz和{sample}_2.fastq.gz，正向和反向序列结束的名称。例如，名为"Healthy_control"的示例必须具有同名的目录和名为 Healthy_control_1.fastq.gz 和 Healthy_control_2.fastq.gz 的 FASTQ 文件。但是，如果文库测序是单端策略，则只应保存一个读取文件以供下游分析。例如，同一示例"正常对照"必须具有名为 Healthy_control.fastq.gz 的唯一 FASTQ 文件。 - 创建包含所有样本名称的表型文件：将第一列命名为"样本"，将第二列命名为"类"。用样本名称填充"样本"列，样本名称必须与样本目录的名称相同，并在"类"列中填充每个样本的表型组（例如，对照组或感染）。最后，保存一个名为"metadata.tsv"的文件，并将其发送到 /home/transcriptome-pipeline/data/ 目录。查看现有的 metadata.tsv 以了解表型文件的格式。
cp /opt/transferdata/metadata.tsv
/home/transcriptome-pipeline/data/metadata.tsv - 访问步骤 1.5.4 中定义的 主机路径 目录，并复制新的结构化目录示例。最后，将示例从 /opt/transferdata 移动到管道数据目录。
cp -rf /opt/transferdata/reads/*
/首页/转录组-管道/数据/读取/
- 将排序数据组织到名为 Reads 的文件夹中。
- 观察所有读取都存储在文件夹 /home/transcriptome-pipeline/data/reads 中。
2. 数据的质量控制
注：以图形方式评估序列读取中出错的概率。删除所有技术序列，例如适配器。
- 使用 FastQC 工具访问文库的测序质量。
- 要生成质量图表，请运行 fastqc 程序。执行命令：
巴什·FastQC.sh
注意：结果将保存在/home/transcriptome-pipeline/results/FastQC文件夹中。由于序列适配器用于文库制备和测序，因此在某些情况下，适配器序列的片段可能会干扰映射过程。
- 要生成质量图表，请运行 fastqc 程序。执行命令：
- 删除适配器序列和低质量读数。访问 "脚本" 文件夹并执行 Trimmomatic 工具的命令：
cd /home/transcriptome-pipeline/scripts
bash trimmomatic.sh
注：用于测序滤波的参数为：去除前导低质量或3个碱基（低于质量3）（领先：3）;删除尾随的低质量或3个基数（低于质量3）（尾随：3）;使用4底座宽的滑动窗口扫描读数，当每个底座的平均质量降至20以下时进行切割（滑动窗口：4：20）;Drop的读数低于36个碱基长度（MINLEN：36）。可以通过编辑 Trimmomatic 脚本文件来更改这些参数。- 确保结果保存在以下文件夹中：/home/transcriptome-pipeline/results/trimreads。执行命令：
ls /home/transcriptome-pipeline/results/trimreads
- 确保结果保存在以下文件夹中：/home/transcriptome-pipeline/results/trimreads。执行命令：
3. 样本的映射和注释
注意：在获得高质量的读数后，需要将这些读数映射到参考基因组。对于此步骤，使用 STAR 映射器来映射示例示例。STAR 映射器工具需要 32 GB 的 RAM 内存来加载和执行读取和基因组图谱。对于没有 32 GB RAM 内存的用户，可以使用已映射的读取。在这种情况下，请跳转到步骤 3.3 或使用 Bowtie2 映射器。此部分包含 STAR（结果显示在所有图中）和 Bowtie2（低内存所需的映射器）的脚本。
- 首先为参考基因组编制索引以进行映射过程：
- 使用命令行访问 "脚本" 文件夹：
cd /home/transcriptome-pipeline/scripts - 对于 STAR 映射器，请执行：
bash indexGenome.sh - 对于 Bowtie 映射器，请执行：
巴什 indexGenomeBowtie2.sh
- 使用命令行访问 "脚本" 文件夹：
- 执行以下命令将过滤后的读数（从步骤2中获得）映射到参考基因组（GRCh38版本）。STAR 和 Bowtie2 映射器都是使用默认参数执行的。
- 对于 STAR 映射器，请执行：
bash mapSTAR.sh - 对于 Bowtie2 映射器，执行：
bash mapBowtie2.sh
注意：最终结果是存储在 /home/transcriptome-pipeline/results/mapreads 中的每个样本的二进制对齐映射 （BAM） 文件。
- 对于 STAR 映射器，请执行：
- 使用 FeatureCounts 工具对映射的读数进行注释，以获取每个基因的原始计数。运行为读取内容添加批注的脚本。
注意：特征计数工具负责将映射的测序读取分配给基因组特征。在生物学问题之后可以更改的基因组注释的最重要方面包括，检测同种型，多个映射读取和外显子 - 外显子连接，对应于参数，GTF.attrType="gene_name"基因或不指定元特征水平的参数，允许MultiOverlap=TRUE和juncCounts=TRUE。- 使用命令行访问脚本文件夹：
cd /home/transcriptome-pipeline/scripts - 要注释映射的读数以获取每个基因的原始计数，请执行命令行：
Rscript 注释。R
注：用于注释过程的参数为：返回基因短名称（GTF.attrType="gene_name"）;允许多个重叠（允许多重叠 = TRUE）;并指示库已配对端 （isPairedEnd=TRUE）。对于单端策略，请使用参数 isPairedEnd=FALSE。结果将保存在/home/transcriptome-pipeline/countreads文件夹中。
- 使用命令行访问脚本文件夹：
- 使基因表达正常化。
注意：归一化基因表达对于比较结果（例如，健康和感染样本）至关重要。还需要归一化来执行微扰分析的共表达和分子程度。- 使用命令行访问 "脚本" 文件夹：
cd /home/transcriptome-pipeline/scripts - 使基因表达正常化。为此，请执行命令行：
Rscript 规范化样本。R
注意：在此实验中，原始计数表达式使用 M 值的修剪平均值 （TMM） 和每百万计数 （CPM） 方法进行归一化。本步骤旨在通过进行文库大小归一化来消除由于技术影响而导致的基因表达差异。结果将保存在/home/transcriptome-pipeline/countreads文件夹中。
- 使用命令行访问 "脚本" 文件夹：
4. 差异表达基因和共表达基因
- 使用开源 EdgeR 包鉴定差异表达的基因。这涉及找到与对照相比表达更高或更低的基因。
- 使用命令行访问 "脚本" 文件夹：
cd /home/transcriptome-pipeline/scripts - 要识别差异表达的基因，请使用命令行执行DEG_edgeR R 脚本：
Rscript DEG_edgeR.R
注意：包含差异表达基因的结果将保存在/home/transcriptome-pipeline/results/degs文件夹中。数据可以传输到个人计算机。
- 使用命令行访问 "脚本" 文件夹：
- 从 csblusp/transcriptome 容器下载数据。
- 将已处理的数据从 /home/transcriptome 管道传输到 /opt/transferdata 文件夹（本地计算机）。
- 通过执行命令行将所有文件复制到本地计算机：
cp -rf /home/transcriptome-pipeline/results /opt/transferdata/pipeline
cp -rf /home/transcriptome-pipeline/data /opt/transferdata/pipeline
注意：现在，转到本地计算机，以确保所有结果、数据集和数据都可以在主机路径中下载。
- 识别共表达模块。
- 访问共表达模块识别工具 （CEMiTool） 网站 （表格
材料）。此工具从用户提供的表达式数据集中识别共表达式模块。在主页上，单击右上角的 "运行 "。这将打开一个新页面来上传表达式文件。 - 单击"表达文件"部分下方的"选择文件"，然后从主机路径上传归一化的基因表达矩阵"tmm_expression.tsv"。
注意：步骤 4.4.是非强制性的。
- 访问共表达模块识别工具 （CEMiTool） 网站 （表格
- 探索共表达模块的生物学意义。
- 单击"示例表型"部分中的"选择文件"，然后从"下载数据"步骤 4.2.2 metadata_cemitool.tsv 上传包含示例表型的文件。执行基因集富集分析 （GSEA）。
- 按"基因相互作用"部分中的"选择文件"以上传具有基因相互作用的文件（cemitool-interactions.tsv）。可以使用webCEMiTool提供的基因相互作用文件作为示例。相互作用可以是蛋白质 - 蛋白质相互作用，转录因子及其转录基因或代谢途径。此步骤为每个共表达模块生成一个交互网络。
- 单击"基因集"部分中的"选择文件"，以上传基因基质转置（GMT）格式文件中功能相关的基因列表。基因集文件使工具能够对每个共表达模块执行富集分析，即过度表示分析（ORA）。
注意：此基因列表可以包含通路、GO 项或 miRNA 靶基因。研究人员可以使用血液转录模块（BTM）作为该分析的基因集。BTM 文件 （BTM_for_GSEA.gmt）。
- 设置用于执行共表达分析的参数并获得其结果。
- 接下来，通过单击加号展开" 参数 "部分以显示默认参数。如有必要，请更改它们。选中 应用 VST 复选框。
- 在" 电子邮件 "部分中编写电子邮件，以电子邮件形式接收结果。此步骤是可选的。
- 按 "运行 CEMiTool" 按钮。
- 点击右上角的下载完整报告， 下载完整的分析报告 。它将cemitool_results.zip下载压缩文件。
- 使用 WinRAR 提取cemitool_results.zip的内容。
注意：包含提取内容的文件夹包含多个文件，其中包含所有分析结果及其建立的参数。
5. 样品分子扰动程度的测定
- 分子扰动度 （MDP） 网络版。
- 要运行 MDP，请访问 MDP 网站（材料表）。MDP计算每个样品与参考的分子距离。单击" 运行 "按钮。
- 在 "选择文件" 链接上，将表达式文件tmm_expression.tsv。然后，从下载数据步骤 4.2.2 上传表型数据文件 metadata.tsv。也可以提交GMT格式的通路注释文件，以计算与疾病相关的通路的扰动评分。
- 上传数据后，定义包含 MDP 使用的表型信息的"类"列。然后，通过选择与控件类对应的标签来定义控件类。
注意：有一些可选参数会影响示例分数的计算方式。如有必要，用户可以更改统计平均方法，标准偏差和扰动基因的最高百分比。 - 之后，按 "运行MDP "按钮，将显示MDP结果。用户可以通过单击每个 图中的"下载图" 以及"下载 MDP 分数 文件"按钮上的 MDP 分数来下载 数字。
注意：如果对如何提交文件或 MDP 的工作原理有疑问，只需浏览教程和关于网页即可。
6. 功能丰富分析
- 创建一个下调 DG 列表和另一个上调 DEX 列表。基因名称必须根据Entrez基因符号。列表中的每个基因必须放在一行上。
- 以 txt 或 tsv 格式保存基因列表。
- 访问 Enrichr 网站（材料表）以执行功能分析。
- 通过单击" 选择文件"选择基因列表。选择其中一个 DEG 列表，然后按 "提交 "按钮。
- 单击网页顶部的 Pathways ，使用 ORA 方法执行功能富集分析。
- 选择路径数据库。"Reactome 2016"通路数据库被广泛用于获取人类数据的生物学意义。
- 再次单击路径数据库的名称。选择 条形图 并检查是否按 p 值排名排序。如果没有，请单击条形图，直到按 p 值排序。此条形图包括根据 p 值的前 10 个路径。
- 按 配置 按钮，选择红色进行上调基因分析，选择蓝色进行下调基因分析。通过单击 svg、 png 和 jpg，以多种格式保存条形图。
- 选择 表 ，然后单击条形图左下角 的将条目导出到表中 ，以获取txt文件中的功能丰富分析结果。
注意：此功能富集结果文件在每行中包含一个通路的名称、提交的DEG列表和通路之间的重叠基因数、p值、调整后的p值、优势比、组合评分以及DEG列表中存在的参与通路的基因的基因符号。 - 对其他 DEG 列表重复相同的步骤。
注意：使用下调DEG的分析为下调基因提供了富集的途径，使用上调基因的分析为上调基因提供了富集的途径。
Subscription Required. Please recommend JoVE to your librarian.
Representative Results
转录组分析的计算环境是在 Docker 平台上创建和配置的。这种方法允许初学者Linux用户在没有先验管理知识的情况下使用Linux终端系统。Docker 平台使用主机操作系统的资源来创建包含特定用户工具的服务容器（图 1B）。创建了一个基于Linux OS Ubuntu 20.04发行版的容器，并完全配置为转录组分析， 可通过 命令行终端访问。在此容器中，有一个用于数据集和脚本的预定义文件夹结构，这是所有管道分析所必需的（图 1C）。我们的研究小组7 发表的一项研究用于分析，其中包括来自健康个体的20个样本和来自CHIKV急性感染个体的39个样本（图1D）。
总RNA测序过程可能会产生读取错误，这可能是由具有两个或多个转录本的簇或试剂耗尽引起的。测序平台返回一组"FASTQ"文件，其中包含每个核苷酸碱基的序列（读取）和相关质量（图2A）。Phred质量量表表示每个碱基读数不正确的概率（图2B）。低质量的读数可能会产生偏倚或不正确的基因表达，从而触发下游分析的连续错误。开发了Trimomatic等工具来识别和删除样品中的低质量读数，并增加映射读取的可能性（图2C，D）。
该绘图模块预先配置了STAR对准器和GRCh38人类宿主作为参考基因组。在此步骤中，从上一步中恢复的高质量读数用作输入，以与人类参考基因组对齐（图3A）。STAR 对准器以 BAM 格式文件输出映射到参考基因组的对齐。基于这种对齐方式，FeatureCounts工具使用GTF文件格式的人类宿主的参考注释来执行这些对齐读取的特征（基因）的注释（图3B）。最后，生成将每个基因名称作为一行，并将每个样品作为一列的表达矩阵（图3C）。还需要提供包含样本名称和相应样本组的附加元数据文件，以便进一步进行下游分析。基因表达矩阵表示映射到样本中每个基因的计数数量，可用作EdgeR输入以识别DEG。此外，使用TMM和CPM对该基因表达基质进行归一化，以消除技术变异性，并通过考虑样品中表达基因在总文库大小中的比例来纠正RNA-seq测量。该矩阵进一步用作共表达和MDP分析的输入。
CEMiTool识别并分析共表达模块12。位于同一模块中的基因是共表达的，这意味着它们在数据集的样本中表现出相似的表达模式。该工具还允许探索每个已识别模块的生物学意义。为此，它提供了三个可选分析 - GSEA 的功能富集分析、过度表示分析 （ORA） 的功能丰富分析和网络分析。GSEA的功能富集分析提供了有关每个模块在每种表型上的基因表达的信息（图4A）。根据这一点，它能够识别在每种表型上被抑制或诱导的模块。ORA分析显示了按调整后的p值排序的每个模块的前10个显着富集的生物学功能。可以将GSEA和ORA结果结合起来，以确定受损的生物过程，以及它们是否被感兴趣的表型抑制或诱导。网络分析提供了每个模块的交互组（图4A）。它能够可视化每个模块的基因如何相互作用。除此之外，网络分析还提供了有关连接最紧密的基因（中心）的信息，这些基因由它们在网络中的名称来识别。节点的大小表示连接程度。
为了识别 DEG，我们开发了一个内部脚本，以单向和简洁的命令行运行端到端差异分析。该脚本执行执行 DEG 分析所需的所有步骤，比较用户在元数据文件中提供的不同示例组。此外，DEG结果存储在下调和上调基因的单独列表中，然后使用Bioconductor的增强型Volcano R包编译成可发表的图（图4B）。
通过 MDP 工具对分子扰动程度进行分析，我们可以识别来自健康和受感染个体的扰动样本11。考虑每个CHIKV感染样本的所有表达基因并考虑健康样本作为参考组，计算扰动评分（图5A）。MDP还仅使用这些样本中最受干扰的基因的前25%进行分析（图5B）。鉴于遗传背景，年龄，性别或其他先前疾病，样本可以呈现出很大的变异性。这些因素可以改变转录组的概况。基于此，MDP建议哪些样本是潜在的生物异常值，以去除它们并改善下游结果（图5A，B）。
ORA的功能富集分析可以使用富集剂进行，以确定DEGs的生物学意义。基于下调基因列表提供的结果表明所研究表型中抑制的生物过程，而基于上调基因列表提供的结果则呈现了在感兴趣的表型中诱导的生物过程。Enrichr生成的条形图中显示的生物过程是基于p值排名的前10个富集基因集（图6）。
图 1：环境 Docker 和示例研究。 （A） Docker 平台使用 OS Host 资源为 Linux 系统创建"容器"，其中包含用于转录组分析的工具。（B） Docker 容器模拟 Linux 系统来执行管道脚本。（C） 创建并组织了转录组管道文件夹结构，以存储用于分析的数据集和脚本。（D）我们小组的研究被用作转录组分析的一个例子。 请点击此处查看此图的放大版本。
图2：测序的质量控制。 （A）FASTQ格式文件用于表示序列和核苷酸碱基质量。（B） Phred 分数方程，其中每 10 个增加一个误读基数的对数概率。（C）和（D）Boxplot分别表示Trimomatic执行之前和之后每个核苷酸碱基的质量分布。 请点击此处查看此图的放大版本。
图3：从序列到基因计数表达的映射和注释过程。 （A）映射包括对齐转录本中的序列和基因组中的序列，以识别基因组定位。（B）映射到参考基因组的读数是根据其重叠的基因组定位进行注释的。（C）基于特征计数等映射文件工具，对基因表达进行汇总。 请点击此处查看此图的放大版本。
图4：共表达基因网络和DEG的统计分析（A）基于模块基因表达的共表达模块和蛋白质 - 蛋白质相互作用网络。（B）CHIKV急性感染和健康个体的统计分析，以及红色（p值和log2FC标准），紫色（仅p值），绿色（仅log2FC）和灰色（无显着性）的差异基因表达。请点击此处查看此图的放大版本。
图5：CHIKV急性感染和健康个体的分子扰动度（MDP ）。（A）使用转录组所有表达基因的每种样品的MDP评分。（B）每个样本的MDP评分仅使用最扰动基因的前25%。 请点击此处查看此图的放大版本。
图6：DEGs的功能分析（A）上调和（B）下调基因被提交给 Enrichr 网站工具，以评估生物学途径或代表性基因集。计算每种途径的P值，图中仅显示显着差异。请点击此处查看此图的放大版本。
Subscription Required. Please recommend JoVE to your librarian.
Discussion
测序文库的制备是以最佳方式回答生物学问题的关键一步。该研究感兴趣的转录本类型将指导选择哪种类型的测序文库并推动生物信息学分析。例如，根据测序的类型，从病原体和宿主相互作用的测序中，可以识别来自两者的序列或仅来自宿主转录本的序列。
下一代测序设备，例如Illumina平台，测量测序质量评分，这代表了碱基被错误调用的概率。下游分析对低质量序列非常敏感，并导致基因表达不足或误读。执行正确分析和解释的另一个障碍是适配器序列。适配器序列有助于文库的制备和测序，在大多数情况下，适配器也是测序的。最近的研究发现，制图工具对最终结果的影响很小13。然而，在病原体 - 宿主研究中，当测试不同的阈值以最小化多映射位点序列问题时，映射过程可以产生稍微更好的结果。
差异基因表达结果应谨慎解释，特别是当每组样本数量非常少且样本来自不同的测定并受到批次效应干扰DEG结果时。这些结果对几个因素很敏感：（一）所应用的数据过滤，例如去除低表达基因和要维持的样本数量;（ii）研究设计，仅比较样本组或每个感染患者与所有对照患者，如CHIKV研究7所示;（iii）用于识别DEGs的统计方法。在这里，我们用Edger举一个基本示例，以识别假设阈值p值为0.05的DEGs。文献中还知道，与其他基准测试方法相比，EdgeR在识别DEGs14方面可以具有很大的变异性。人们可以考虑这些不同方法之间的权衡，并考虑可用的重复次数和实验设计的复杂性14。
CEMiTool执行共表达模块分析12。该工具可通过Bioconductor存储库上的R包获得，也可以通过webCEMiTool以用户友好的版本获得;后者是当前协议中使用的版本。与WGCNA15 相关的替代软件16相比，它具有几个优点，包括它更人性化17。此外，该工具具有过滤基因的自动方法，而在WGCNA中，用户必须在使用WGCNA之前过滤基因。此外，此工具已建立默认参数，而在 WGCNA 中，用户必须手动选择参数分析。手动参数选择会损害再现性;因此，自动参数选择保证了更高的再现性。
在某些情况下，CEMiTool无法找到合适的软阈值，也称为β值。在这种情况下，用户应检查RNA-seq数据是否具有较强的均值-方差依赖性。如果均值与方差（考虑所有基因）表现出很强的线性关系，则用户必须重新运行分析，检查"应用VST"参数以消除转录组学数据的均值-方差依赖性。检查数据中是否存在较强的均值-方差依赖关系并在存在时将其删除始终至关重要。
CEMiTool已被广泛用于识别和探索共表达模块的生物学意义。一项 CHIKV 急性感染研究显示，在出现症状 2 至 4 天后，患者中具有较高活性的模块7。ORA对该模块的功能富集表现出单核细胞和嗜中性粒细胞的增加7。一项流感疫苗接种研究使用从基线到疫苗接种后第7天的血液转录组，提出了功能丰富的共表达模块，用于与T，B和自然杀伤细胞，单核细胞，嗜中性粒细胞，干扰素反应和血小板活化相关的生物过程18。
考虑到转录组数据集的变异性，识别和量化数据异质性可能是一个挑战，因为许多变量都会影响基因表达谱7，11。MDP提供了一种通过以下步骤来识别和量化来自健康和受感染受试者的扰动样本的方法：（i）计算对照样本的中心性方法（中位数或平均值）和标准偏差;（ii）使用获得的值计算所有基因的z评分;（iii）设置阈值z得分绝对大于2，表明与对照样品的代表性偏差;（iv）使用为每个样本过滤的分数计算基因值的平均值。尽管 scRNA-seq 分析存在一些局限性，但该工具在确定微阵列和 RNA-seq 数据的扰动评分方面具有功能11。此外，之前的一项研究已经使用该工具证明了结核病和糖尿病患者血液转录组的分子扰动程度升高19。在这项工作中，已经显示了使用健康个体作为参考组的对照组和CHIKV急性感染样本的扰动。
Enrichr执行的功能富集分析是ORA20，21。ORA 是一种类型的功能扩充分析，其中用户必须向工具提供 DEG 列表。DEG列表通常在下调的DEG列表和上调的DEG列表中分开。还有其他工具可以执行ORA，其中包括gProfiler，它以用户友好的Web版本22 提供，goseq23 在Bioconductor上作为R包提供。另一种类型的功能富集分析是GSEA。要执行GSEA，用户必须提供排名列表中的所有基因。该列表通常根据折叠变化中的基因表达进行排名。
富集器始终根据条形图结果中的 p 值提供富集的前 10 个基因集。因此，用户在解读结果时必须警惕，如果富集基因集少于10个，条形图也会显示非富集的生物过程。为了避免这种错误，用户必须为p值建立一个截止值，并在假设条形图的所有基因集都被富集之前观察通路的p值。此外，用户必须意识到，条形图中显示的10个基因组的顺序是根据p值，而不是调整后的p值。如果用户想要在条形图中显示所有丰富的路径，甚至根据调整后的 p 值重新排序，建议用户使用下载的表格创建自己的条形图。用户可以使用Excel甚至R软件制作新的条形图。
Subscription Required. Please recommend JoVE to your librarian.
Disclosures
作者没有什么可透露的。
Acknowledgments
HN由FAPESP（资助编号：#2017/50137-3，2012/19278-6，2018/14933-2，2018/21934-5和2013/08216-2）和CNPq（313662/2017-7）资助。
我们特别感谢为研究员提供的以下赠款：ANAG（FAPESP Process 2019/13880-5），VEM（FAPESP Process 2019/16418-0），IMSC（FAPESP Process 2020/05284-0），APV（FAPESP Process 2019/27146-1）和RLTO（CNPq Process 134204/2019-0）。
Materials
| Name | Company | Catalog Number | Comments |
| CEMiTool | Computational Systems Biology Laboratory | 1.12.2 | Discovery and the analysis of co-expression gene modules in a fully automatic manner, while providing a user-friendly HTML report with high-quality graphs. |
| EdgeR | Bioconductor (Maintainer: Yunshun Chen [yuchen at wehi.edu.au]) | 3.30.3 | Differential expression analysis of RNA-seq expression profiles with biological replication |
| EnhancedVolcano | Bioconductor (Maintainer: Kevin Blighe [kevin at clinicalbioinformatics.co.uk]) | 1.6.0 | Publication-ready volcano plots with enhanced colouring and labeling |
| FastQC | Babraham Bioinformatics | 0.11.9 | Aims to provide a simple way to do some quality control checks on raw sequence data coming from high throughput sequencing |
| FeatureCounts | Bioinformatics Division, The Walter and Eliza Hall Institute of Medical Research | 2.0.0 | Assign mapped sequencing reads to specified genomic features |
| MDP | Computational Systems Biology Laboratory | 1.8.0 | Molecular Degree of Perturbation calculates scores for transcriptome data samples based on their perturbation from controls |
| R | R Core Group | 4.0.3 | Programming language and free software environment for statistical computing and graphics |
| STAR | Bioinformatics Division, The Walter and Eliza Hall Institute of Medical Research | 2.7.6a | Aligner designed to specifically address many of the challenges of RNA-seq data mapping using a strategy to account for spliced alignments |
| Bowtie2 | Johns Hopkins University | 2.4.2 | Ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences |
| Trimmomatic | THE USADEL LAB | 0.39 | Trimming adapter sequence tasks for Illumina paired-end and single-ended data |
| Get Docker | Docker | 20.10.2 | Create a bioinformatic environment reproducible and predictable (https://docs.docker.com/get-docker/) |
| WSL2-Kernel | Windows | NA | https://docs.microsoft.com/en-us/windows/wsl/wsl2-kernel |
| Get Docker Linux | Docker | NA | https://docs.docker.com/engine/install/ubuntu/ |
| Docker Linux Repository | Docker | NA | https://docs.docker.com/engine/install/ubuntu/#install-using-the-repository |
| MDP Website | Computational Systems Biology Laboratory | NA | https://mdp.sysbio.tools |
| Enrichr Website | MaayanLab | NA | https://maayanlab.cloud/Enrichr/ |
| webCEMiTool | Computational Systems Biology Laboratory | NA | https://cemitool.sysbio.tools/ |
| gProfiler | Bioinformatics, Algorithmics and Data Mining Group | NA | https://biit.cs.ut.ee/gprofiler/gost |
| goseq | Bioconductor (Maintainer: Matthew Young [my4 at sanger.ac.uk]) | NA | http://bioconductor.org/packages/release/bioc/html/goseq.html |
| SRA NCBI study | NCBI | NA | https://www.ncbi.nlm.nih.gov/bioproject/PRJNA507472/ |
References
- Weaver, S. C., Charlier, C., Vasilakis, N., Lecuit, M. Zika, Chikungunya, and Other Emerging Vector-Borne Viral Diseases. Annual Review of Medicine. 69, 395-408 (2018).
- Burt, F. J., et al. Chikungunya virus: an update on the biology and pathogenesis of this emerging pathogen. The Lancet. Infectious Diseases. 17 (4), 107-117 (2017).
- Hua, C., Combe, B.
- Suhrbier, A., Jaffar-Bandjee, M. -C., Gasque, P.
- Nakaya, H. I., et al. Gene profiling of chikungunya virus arthritis in a mouse model reveals significant overlap with rheumatoid arthritis. Arthritis and Rheumatism. 64 (11), 3553-3563 (2012).
- Michlmayr, D., et al. Comprehensive innate immune profiling of chikungunya virus infection in pediatric cases. Molecular Systems Biology. 14 (8), 7862 (2018).
- Soares-Schanoski, A., et al. Systems analysis of subjects acutely infected with the Chikungunya virus. PLOS Pathogens. 15 (6), 1007880 (2019).
- Alexandersen, S., Chamings, A., Bhatta, T. R. SARS-CoV-2 genomic and subgenomic RNAs in diagnostic samples are not an indicator of active replication. Nature Communications. 11 (1), 6059 (2020).
- Wang, D., et al. The SARS-CoV-2 subgenome landscape and its novel regulatory features. Molecular Cell. 81 (10), 2135-2147 (2021).
- Wilson, J. A. C., et al. RNA-Seq analysis of chikungunya virus infection and identification of granzyme A as a major promoter of arthritic inflammation. PLOS Pathogens. 13 (2), 1006155 (2017).
- Gonçalves, A. N. A., et al. Assessing the impact of sample heterogeneity on transcriptome analysis of human diseases using MDP webtool. Frontiers in Genetics. 10, 971 (2019).
- Russo, P. S. T., et al. CEMiTool: a Bioconductor package for performing comprehensive modular co-expression analyses. BMC Bioinformatics. 19 (1), 56 (2018).
- Costa-Silva, J., Domingues, D., Lopes, F. M. RNA-Seq differential expression analysis: An extended review and a software tool. PloS One. 12 (12), 0190152 (2017).
- Seyednasrollah, F., Laiho, A., Elo, L. L. Comparison of software packages for detecting differential expression in RNA-seq studies. Briefings in Bioinformatics. 16 (1), 59-70 (2015).
- Zhang, B., Horvath, S. A general framework for weighted gene co-expression network analysis. Statistical Applications in Genetics and Molecular Biology. 4, Article17 (2005).
- Cheng, C. W., Beech, D. J., Wheatcroft, S. B. Advantages of CEMiTool for gene co-expression analysis of RNA-seq data. Computers in Biology and Medicine. 125, 103975 (2020).
- Cardozo, L. E., et al. webCEMiTool: Co-expression modular analysis made easy. Frontiers in Genetics. 10, 146 (2019).
- de Lima, D. S., et al. Long noncoding RNAs are involved in multiple immunological pathways in response to vaccination. Proceedings of the National Academy of Sciences of the United States of America. 116 (34), 17121-17126 (2019).
- Prada-Medina, C. A., et al. Systems immunology of diabetes-tuberculosis comorbidity reveals signatures of disease complications. Scientific Reports. 7 (1), 1999 (2017).
- Chen, E. Y., et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 14, 128 (2013).
- Kuleshov, M. V., et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Research. 44, 90-97 (2016).
- Raudvere, U., et al. g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Research. 47, 191-198 (2019).
- Young, M. D., Wakefield, M. J., Smyth, G. K., Oshlack, A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology. 11 (2), 14 (2010).
