利用RNA-seq研究分子进化与基因表达的生物信息学分析流程

摘要

本实验方案的目的是利用RNA测序数据研究候选基因的进化与表达。

摘要

整理和报告大规模数据集（如全基因组或转录组数据）通常是一项艰巨的任务。一种分解结果的方法是聚焦于一个或多个对研究生物体和研究主题具有重要意义的基因家族。在本实验方案中，我们概述了用于构建系统发育树并量化目标基因表达水平的生物信息学步骤。系统发育树有助于揭示基因在种内和种间的进化模式，并可推断直系同源关系。结合RNA-seq数据，可进一步比较这些基因在不同个体或组织中的表达情况，从而增强分析结果。分子进化与表达研究能够揭示基因在不同物种间的进化模式及其功能保守性。对一个基因家族的表征可为后续研究提供基础，并可在新发表的基因组或转录组研究论文中突出一个重要基因家族的存在与意义。

引言

测序技术的进步促进了非模式生物基因组和转录组的测序。除了从多种生物中测序DNA和RNA的可行性不断提高外，大量公开可用的数据也为研究目标基因提供了便利。本实验方案的目的是提供一系列生物信息学步骤，用于探究可能在目标生物中发挥重要作用的基因的分子进化与表达特征。

研究一个基因或基因家族的进化过程，可以为生物系统的演化提供重要见解。基因家族的成员通常通过识别保守基序或同源基因序列来确定。以往对基因家族进化的研究主要利用远缘关系的模式生物基因组进行¹。这种方法的一个局限性在于，尚不清楚这些基因家族在近缘物种中的演化过程，以及不同环境选择压力所起的作用。在本实验方案中，我们纳入了在近缘物种中搜索同源基因的步骤。通过构建门水平的系统发育树，我们可以观察到基因家族进化中的趋势，例如保守基因的存在或谱系特异性的基因重复事件。在此层次上，我们还可以分析基因是直系同源（orthologs）还是旁系同源（paralogs）。尽管许多同源基因可能具有相似的功能，但这并非绝对²。在这些研究中引入系统发育树对于判断同源基因是否为直系同源至关重要。在真核生物中，许多直系同源基因在细胞内保持相似功能，这一点已得到证实，例如哺乳动物蛋白能够恢复酵母中直系同源蛋白的功能³。然而，也存在非直系同源基因执行特定功能的情况⁴。

系统发育树开始揭示基因与物种之间的关系，但仅基于遗传关系无法确定基因功能。基因表达研究结合功能注释和富集分析 可为基因功能提供有力支持。在能够对个体之间或组织类型之间的基因表达进行定量和比较的情况下，更能揭示其潜在功能。以下实验方案参考了研究Hydra vulgaris⁷中视蛋白基因所采用的方法，但这些方法可应用于任何物种和任何基因家族。此类研究结果为在非模式生物中进一步探究基因功能和基因调控网络奠定了基础。例如，通过分析视蛋白（opsins）的系统发育关系——这类蛋白可启动光转导级联反应——有助于理解眼睛和感光功能的演化过程^8,9,10,11。在此背景下，非模式生物尤其是刺胞动物或栉板动物等基础动物类群，有助于阐明光转导级联和视觉系统在不同分支中的保守性或演化变化^12,13,14。类似地，确定其他基因家族的系统发育关系、表达模式及其调控网络，将有助于我们理解适应性状背后的分子机制。

方案

本方案遵循加州大学欧文分校的动物护理指南。

1. RNA-seq 文库构建

1. 使用以下方法分离RNA。
   1. 收集样品。若需在后续提取RNA，应将样品迅速冷冻或置于RNA保存液中¹⁵（材料表）。
   2. 处死并解剖生物体，分离目标组织。
   3. 使用RNA提取试剂盒提取总RNA，并使用RNA纯化试剂盒进行纯化（材料表）
      注意：针对不同物种和组织类型，可能存在更适用的实验方案和试剂盒^16,17。我们已从蝴蝶的不同体组织¹⁸以及胶状水母Hydra¹⁹中成功提取RNA（见讨论部分）。
   4. 测定各样品RNA的浓度和质量（材料表）。构建cDNA文库时应选用RNA完整性数值（RIN）高于8、理想情况下接近9的样品²⁰。
2. 按以下步骤构建cDNA文库并进行测序。
   1. 根据文库构建说明书操作手册构建cDNA文库（见讨论部分）。
   2. 测定cDNA的浓度和质量（材料表）。
   3. 对文库进行多重化（multiplex）并进行测序。

2. 访问计算机集群

注意：RNA-seq 分析需要处理大型文件，最好在计算机集群上进行（材料表）。

1. 在终端（Mac）或 PuTTY（Windows）应用程序窗口中，使用命令 ssh username@clusterlocation 登录计算机集群账户。

3. 获取RNA-seq读段

1. 从测序机构获取 RNA-seq 测序读段，或对于已发表研究中的数据，从其存放的数据存储库中获取（3.2 或 3.3）。
2. 从 ArrayExpress 等存储库下载数据时，请执行以下操作：
   1. 使用登录编号在网站上进行搜索。
   2. 找到下载数据的链接，单击鼠标左键并选择 复制链接。
   3. 在终端窗口中输入 wget，然后选择 粘贴链接，将数据复制到用于分析的目录中。
3. 下载 NCBI 短读段存档（SRA）数据时，请遵循以下替代步骤：
   1. 在终端中使用 wget 下载 SRA Toolkit v. 2.8.1。
      注意：向计算机集群下载和安装程序可能需要 root 权限，若安装失败，请联系您的计算机集群管理员。
   2. 通过输入 tar -xvf $TARGZFILE 完成程序的安装。
   3. 在 NCBI 中搜索您要下载的样本对应的 SRA 登录编号，编号格式应为 SRRXXXXXX。
   4. 在终端窗口中输入 [sratoolkit 位置]/bin/prefetch SRRXXXXXX 获取 RNA-seq 数据。
   5. 对于双端测序文件，输入 [sratoolkit 位置]/bin/fastq-dump --split-files SRRXXXXXX 以获得两个 fastq 文件（SRRXXXXXX_1.FASTQ 和 SRRXXXXXX_2.FASTQ）。
      注意：若要进行 Trinity de novo 组装，请使用命令 [sratoolkit 位置]/bin/fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files SRRXXXXXX

4. 去除接头序列和低质量读段（可选）

1. 在计算集群上安装或加载 Trimmomatic²¹ v. 0.35。
2. 在存放 RNA-seq 数据文件的目录中，输入一条命令，其中包含 trimmomatic jar 文件的位置、输入 FASTQ 文件、输出 FASTQ 文件，以及读长和质量等可选参数。
   注意：具体命令会因原始数据及期望的读长和质量要求而有所不同。对于带有 Nextera 接头的 Illumina 43 bp 读长数据，我们使用的命令为：java -jar /data/apps/trimmomatic/0.35/trimmomatic-0.35.jar PE $READ1.FASTQ $READ2.FASTQ paired_READ1.FASTQ unpaired_READ1.FASTQ paired_READ2.FASTQ unpaired_READ2.FASTQ ILLUMINACLIP:adapters.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:17 MINLEN:30。

5. 获取参考基因组序列

1. 在 Google、EnsemblGenomes 以及 NCBI Genomes 和 Nucleotide TSA（转录组鸟枪法组装）数据库中搜索目标物种的参考基因组或组装完成的转录组（图1）。
   注意：如果参考基因组或转录组不可用或质量较低，请转至第6步，进行从头（de novo）组装。
2. 若存在参考基因组或组装完成的转录组，请将其以 fasta 文件格式下载至分析执行的位置，具体操作步骤如下：
   1. 找到下载基因组的链接，单击鼠标左键并复制链接。
   2. 在终端窗口中输入wget命令，然后粘贴链接地址。如有提供，同时复制该参考基因组的 GTF 文件和蛋白质 FASTA 文件。

6. 进行从头组装（de novo assembly）（替代步骤5）

1. 在终端窗口中输入 cat *READ1.FASTQ > $all_READ1.FASTQ 和 cat *READ2.FASTQ > all_READ2.FASTQ，将所有样本的 RNA-seq READ1 和 READ2 fastq 文件合并。
2. 在计算集群上安装或加载 Trinity²² v.2.8.5。
3. 在终端中输入以下命令进行生成与组装：Trinity --seqType fq --max_memory 20G --left $all_READ1.FASTQ --right $all_READ2.FASTQ。

7. 将序列比对到基因组 (7.1) 或de novo转录组 (7.2)

1. 使用 STAR²³ v. 2.6.0c 和 RSEM²⁴ v. 1.3.0 将测序读段比对到参考基因组。
   1. 在计算集群上安装或加载 STAR v. 2.6.0c 和 RSEM v. 1.3.0。
   2. 通过输入命令 rsem-prepare-reference --gtf $GENOME.GTF --star -p 16 $GENOME.FASTA $OUTPUT 对基因组进行索引构建。
   3. 通过输入命令 rsem-calculate-expression -p 16 --star --paired-end $READ1.FASTQ $READ2.FASTQ $INDEX $OUTPUT 对每个样本的读段进行比对并计算基因表达量。
   4. 使用命令 mv RSEM.genes.results $sample.genes.results 将结果文件重命名为具有描述性的名称。
   5. 通过输入命令 rsem-generate-data-matrix *[genes/isoforms.results] > $OUTPUT 生成所有计数的矩阵。
2. 使用 RSEM 和 bowtie 将 RNA-seq 数据比对到 Trinity de novo 组装结果。
   1. 安装或加载 Trinity²² v.2.8.5、Bowtie²⁵ v. 1.0.0 以及 RSEM v. 1.3.0。
   2. 通过输入命令 [trinity_location]/align_and_estimate_abundance.pl --prep-reference --transcripts $TRINITY.FASTA --seqType fq --left $READ1.FASTQ --right $READ2.FASTQ --est_method RSEM --aln_method bowtie --trinity_mode --output_dir $OUTPUT 对读段进行比对并计算每个样本的表达量。
   3. 使用命令 mv RSEM.genes.results $sample.genes.results 将结果文件重命名为具有描述性的名称。
   4. 通过输入命令 [trinity_location]/abundance_estimates_to_matrix.pl --est_method RSEM *[genes/isoforms].results 生成所有计数的矩阵。

8. 鉴定目的基因

注意：以下步骤可使用核苷酸或蛋白质 FASTA 文件进行，但使用蛋白质序列时效果最佳且更为简便。在不同物种之间进行搜索时，采用蛋白质对蛋白质的 BLAST 比对更有可能获得结果。

1. 对于参考基因组，使用步骤 5.2.2 中的蛋白质 FASTA 文件或参见补充材料以生成自定义基因特征 GTF 文件。
2. 对于从头转录组（de novo transcriptome），使用 TransDecoder 生成蛋白质 FASTA 文件。
   1. 在计算机集群上安装或加载 TransDecoder v. 5.5.0。
   2. 通过输入命令 [Transdecoder 位置]/TransDecoder.LongOrfs -t $TRINITY.FASTA 来查找最长开放阅读框及预测的肽序列。
3. 在 NCBI GenBank 中搜索近缘物种中的同源基因。
   1. 打开一个网页浏览器窗口，访问 https://www.ncbi.nlm.nih.gov/genbank/。
   2. 在搜索栏中输入目标基因名称以及已测序的近缘物种名称、属名或门名；在搜索栏左侧选择“protein”，然后点击“search”。
   3. 点击 Send to，然后选择 File。在“Format”下选择 FASTA，然后点击 Create File 以导出序列。
   4. 通过在本地终端窗口中输入 scp $FASTA username@clusterlocation:/$DIR 将同源基因的 FASTA 文件传输至计算机集群，或使用 FileZilla 在本地计算机与集群之间传输文件。
4. 使用 BLAST+²⁶ 搜索候选基因。
   1. 在计算机集群上安装或加载 BLAST+ v. 2.8.1。
   2. 在计算机集群上，通过输入命令 [BLAST+ 位置]/makeblastdb -in $PEP.FASTA -dbtype prot -out $OUTPUT，基于基因组或转录组翻译得到的蛋白质 FASTA 文件创建 BLAST 数据库。
   3. 通过输入命令 [BLAST+ 位置]/blastp -db $DATABASE -query $FASTA -evalue 1e-10 -outfmt 6 -max_target_seqs 1 -out $OUTPUT，将来自 NCBI 的同源基因序列与目标物种数据库进行比对。
   4. 使用 more 命令查看输出文件。复制目标物种中唯一的基因 ID 至一个新的文本文件。
   5. 通过输入命令 perl -ne 'if(/^>(\S+)/){$c=$i{$1}}$c?print:chomp;$i{$_}=1 if @ARGV' $gene_id.txt $PEP.FASTA > $OUTPUT 提取候选基因的序列。
5. 通过反向 BLAST 验证基因注释。
   1. 在网页浏览器中访问 https://blast.ncbi.nlm.nih.gov/Blast.cgi。
   2. 选择 tblastn，粘贴候选序列，选择“Non-redundant protein sequence”数据库，然后点击 BLAST。
6. 通过使用基因本体（GO）术语对基因组或转录组中的所有基因进行注释，识别更多候选基因（参见讨论部分）。
   1. 将蛋白质 FASTA 文件传输至本地计算机。
   2. 下载并安装 Blast2GO^27,28,29 v. 5.2 至本地计算机。
   3. 打开 Blast2GO，点击 File，选择 Load，再选择 Load Sequences，点击 Load Fasta File (fasta)。选择 FASTA 文件后点击 Load。
   4. 点击 Blast，选择 NCBI Blast，然后点击 Next。可编辑参数，或直接点击 Next，再次编辑参数后点击 Run，以查找最相似的基因描述。
   5. 点击“mapping”，然后点击 Run，以搜索相似蛋白的基因本体（Gene Ontology）注释。
   6. 接着点击 interpro，选择 EMBL-EBI InterPro，点击 Next。可编辑参数，或点击 Next 后点击 Run，以搜索已知基因家族和结构域的特征序列。
   7. 点击 File，选择 Export，点击 Export Table。点击 Browse，命名文件，点击 Save，再点击 Export，以导出注释结果。
   8. 在注释表格中搜索感兴趣的 GO 术语，以识别额外的候选基因。从 FASTA 文件中提取相应序列（见步骤 8.4.5）。

9. 系统发育树

1. 下载并安装 MEGA³⁰ v. 7.0.26 到您的本地计算机。
2. 打开 MEGA，点击 对齐，点击 编辑/构建比对，选择 创建新的比对 点击 好的，选择 蛋白质.
3. 当对齐窗口打开时，点击 编辑，点击 从文件插入序列 选择候选基因及可能同源基因的蛋白质序列FASTA文件。
4. 选择所有序列。找到臂状符号并将其悬停。此时应显示“使用 MUSCLE 进行序列比对”³¹ 算法。单击臂状符号，然后单击 比对蛋白质 以比对序列。编辑参数或单击 好的 使用默认参数进行比对。
5. 目视检查并对齐窗口进行必要的手动修改，然后保存并关闭对齐窗口。
6. 在主 MEGA 窗口中，单击 模型，点击 寻找最佳DNA/蛋白质模型（机器学习），选择比对文件并设置相应参数，例如： 分析：模型选择（ML），所用树：自动（邻接法建树），统计方法：最大似然法，替换类型：氨基酸， gaps/缺失数据处理：使用所有位点，分支位点过滤：无。
7. 确定最佳数据模型后，转到 MEGA 主窗口。单击 系统发育 并点击 构建/检验最大似然树 然后根据需要选择对齐方式。为系统发育树选择适当的参数： 统计方法：最大似然法，系统发育检验：采用100次重复的自举法，替代类型：氨基酸，模型：LG加频率（+F），位点间速率：伽马分布（G），分为5个离散伽马类别，缺失数据/空缺处理：使用所有位点，最大似然启发式方法：最近邻交换法（NNI）。

10. 使用 TPM 可视化基因表达

1. 对于 Trinity，在计算机集群中进入运行 abundance_estimates_to_matrix.pl 的目录，其中一个输出文件应为 matrix.TPM.not_cross_norm。将此文件传输至本地计算机。
   注意：交叉样本标准化的详细信息请参见补充材料。
2. 对于基因组分析得到的 TPM 值，请按以下步骤操作。
   1. 在计算机集群中，进入 RSEM 安装目录。通过输入命令 scp rsem-generate-data-matrix rsem-generate-TPM-matrix 复制 rsem-generate-data-matrix 文件。使用 nano 编辑新文件，将“my $offsite = 4”中的数值从 4 修改为 5（适用于 TPM），修改后应为“my $offsite = 5”。
3. 进入包含 RSEM 输出文件 .genes.results 的目录，运行 rsem-generate-TPM-matrix *[genes/isoforms.results] > $OUTPUT 以生成 TPM 矩阵。将结果传输至本地计算机。
4. 使用 ggplot2 可视化结果。
   1. 在本地计算机上下载 R v. 4.0.0 和 RStudio v. 1.2.1335。
   2. 打开 RStudio，在屏幕右侧进入 Packages 选项卡并点击 Install。输入 ggplot2 并点击 install。
   3. 在 R 脚本窗口中，输入以下命令读取 TPM 表格：data<-read.table("$tpm.txt",header = T)
   4. 绘制类似于 图 4 的柱状图，可输入类似以下命令：p<- ggplot() + geom_bar(aes(y=TPM, x=Symbol, fill=Tissue), data=data, stat="identity")
      fill<-c("#d7191c","#fdae61", "#ffffbf", "#abd9e9", "#2c7bb6")
      p<-p+scale_fill_manual(values=fill)
      p + theme(axis.text.x = element_text(angle = 90))

结果

上述方法总结于图1，并应用于Hydra vulgaris组织的数据集。H. vulgaris是一种淡水无脊椎动物，属于刺胞动物门（Cnidaria），该门还包括珊瑚、水母和海葵。H. vulgaris可通过出芽方式进行无性繁殖，并在被切成两段时再生其头部和足部。在本研究中，我们旨在探讨Hydra中视蛋白基因的进化与表达⁷。尽管Hydra没有眼睛，但其表现出光依赖性行为³²。视蛋白基因编码的蛋白质在视觉中起重要作用，可检测不同波长的光并启动光转导级联反应。在基础物种中研究该基因家族的分子进化与表达，有助于理解动物眼睛及光感知能力的演化过程。

我们使用 Hydra 2.0 生成了引导式组装³³ 参考基因组和公开可用的RNA-seq数据（GEO编号GSE127279）

讨论

本实验方案的目的是概述利用RNA-seq数据对基因家族进行表征的各个步骤。这些方法已被证明适用于多种物种和数据集^4,34,35。本文建立的分析流程已尽可能简化，即使是生物信息学初学者也能轻松遵循。本方案的重要意义在于，它详尽列出了完成一项可发表水平分析所需的全部步骤和必要软件工具。方案中的一个关键步骤是必须获得正确组装的全长转录本，这依赖于高质量的基因组或转录组数据。为获得合适的转录本，需要高质量的RNA和/或DNA，以及下文将讨论的良好注释信息。

在进行RNA-seq文库构建时，我们列出了适用于水螅（Hydra）¹⁹和蝴蝶¹⁸等小型生物组织的试剂盒（材料表）。需要注意的是，对于低起始量RNA样本，我们采用了改良的实验流程

材料

本文使用的材料清单

姓名	公司	目录编号	评论
生物分析仪-DNA试剂盒	Agilent	5067-4626	湿实验材料
生物分析仪-RNA试剂盒	Agilent	5067-1513	湿实验材料
BLAST+ v. 2.8.1			在计算集群上 https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
Blast2GO（在您的个人电脑上）			在本地计算机上 https://www.blast2go.com/b2g-register-basic
boost v. 1.57.0			在计算集群上
Bowtie v. 1.0.0			在计算集群上 https://sourceforge.net/projects/bowtie-bio/files/bowtie/1.3.0/
计算集群（强烈推荐）			注意：基因组数据分析最好在高性能计算集群上进行，因为文件非常大。
Cufflinks v. 2.2.1			在计算集群上
edgeR v. 3.26.8（在R中）			在Rstudio中 https://bioconductor.org/packages/release/bioc/html/edgeR.html
gcc v. 6.4.0			在计算集群上
Java v. 11.0.2			在计算集群上
MEGA7（在您的个人电脑上）			在本地计算机上 https://www.megasoftware.net
MEGAX v. 0.1			在本地计算机上 https://www.megasoftware.net
NucleoSpin RNA II试剂盒	Macherey-Nagel	740955.5	湿实验材料
perl 5.30.3			在计算集群上
python			在计算集群上
Qubit 2.0 荧光仪	ThermoFisher	Q32866	湿实验材料
R v.4.0.0			在计算集群上 https://cran.r-project.org/src/base/R-4/
RNA保存液（RNAlater）	ThermoFisher	AM7021	湿实验材料
RNeasy试剂盒	Qiagen	74104	湿实验材料
RSEM v. 1.3.0			计算机软件 https://deweylab.github.io/RSEM/
RStudio v. 1.2.1335			在本地计算机上 https://rstudio.com/products/rstudio/download/#download
Samtools v. 1.3			计算机软件
SRA Toolkit v. 2.8.1			在计算集群上 https://github.com/ncbi/sra-tools/wiki/01.-Downloading-SRA-Toolkit
STAR v. 2.6.0c			在计算集群上 https://github.com/alexdobin/STAR
StringTie v. 1.3.4d			在计算集群上 https://ccb.jhu.edu/software/stringtie/
Transdecoder v. 5.5.0			在计算集群上 https://github.com/TransDecoder/TransDecoder/releases
Trimmomatic v. 0.35			在计算集群上 http://www.usadellab.org/cms/?page=trimmomatic
Trinity v.2.8.5			在计算集群上 https://github.com/trinityrnaseq/trinityrnaseq/releases
TRIzol	ThermoFisher	15596018	湿实验材料
TruSeq RNA文库制备试剂盒v2	Illumina	RS-122-2001	湿实验材料
TURBO DNA-free试剂盒	ThermoFisher	AM1907	湿实验材料
*在计算集群上下载和安装可能需要root权限。请联系您的网络管理员。