3 min read

原核转录组分析的原理和流程

本文是我在华中农业大学JC Bioinformatics 2018暑期生物信息学培训班上的讲义

原核转录组

原核转录组课程讲授利用RNA-seq方法研究细菌基因组中基因表达差异的原理和方法。课程中不仅介绍了RNA-seq生物信息学分析的内容,还对RNA-seq的原理、文库构建等湿实验进行解读。

为了降低学习门槛,课程中的实例取自公开发表的细菌基因组,但是这不意味着相关原理和方法不能应用于其它微生物基因组。

另外,本课程涉及到的文库构建原理、方法、生物信息学工具软件在其它物种相关研究(各种基于高通量测序的研究方法)中也具有普适性。需要学生自己体会和进一步学习。

在开始本课程之前

你需要:

  • 在你的电脑上已经安装了Linux双系统或虚拟机(推荐)
    • Linux发行版推荐Ubuntu
    • 虚拟机软件推荐VirtualBox
  • 对Linux系统下的工作环境有一定认识
    • 使用Bash或类似的命令行工具
    • 了解命令参数的使用
    • 会使用程序内嵌的帮助系统(help或man)
    • 能够使用系统的包管理系统安装软件(如Ubuntu下的apt-get)
    • 能够安装预编译好的生物信息学工具(包括设置PATH环境变量等)
  • 在你的电脑上已经安装了最新版本的R
    • 能够安装CRAN里面的包
    • 能够安装Bioconductor里面的包
    • 会使用R的命令行
    • 会操作R的变量和数据
  • 对ggplot作图有初步的认识
    • 图形语法
    • 能够使用ggplot2绘制简单的图形
  • 在你的电脑上已经安装了最新版本的RStudio
    • 已经初步认识了RStudio的界面
    • 能够在RStudio中使用命令行
    • 能够在RStudio中查看变量和输出

课程主要内容

  • 转录组的原理:通过测序mRNA等的丰度得到基因表达情况
  • RNA-seq的实验和分析流程
    • 文库的构建:ssRNA-seq的介绍
    • 短序列与基因组的比对
    • 基因组注释信息在分析中的应用
    • 差异表达基因的鉴定(差异基因数目上、表达量上的差异)
    • 差异表达基因的富集分析
    • 如何针对数据下一个结论
  • 商业公司RNA-seq结果解析
    • 文件夹的结构
    • 结题报告(有用、无用、错误和正确的信息)
  • 从头开始RNA-seq分析的流程
    • Alignment -> Read Count
    • 靠谱的基因表达差异工具:DESeq2
    • 靠谱的基因富集分析工具:ClusterProfiler
  • 案例研究:重现论文中的数据图表
    • 研究目的:BS、PP对VP基因表达的影响
    • 实验设计:共培养和单独培养的RNA-seq
    • RNA-seq数据分析
    • 结论

转录组的原理:通过测序mRNA的丰度得到基因表达情况

前一段时间,有一个非常火的电视剧:三国机密之潜龙在渊。讲述了汉室与曹操之间纠葛的关系。 里面有曹丕与曹植争宠。让我想到了他们还有个弟弟:曹冲。历史总是很有意思的。 曹冲十三岁就死了,让人遗憾。但是,曹冲称象的故事广为流传。

曹冲称象图

古时候没有那么大量程的秤,采用了化整为零的方法。所以有一头大象 = 一堆石头的换算。

那如果来了一群大象,就可以分别用不同的一堆石头来定量,得出大象之间的重量差异

  • 大象A = 石头堆A
  • 大象B = 石头堆B
  • 大象C = 石头堆C
  • ……

这其实跟RNA-seq检测基因表达差异的情况非常类似。

我们知道中心法则:DNA编码的基因通过转录变成了RNA(其中大量是rRNA)。RNA-seq就是通过比较RNA的丰度来定量基因表达的差异。

中心法则

问题: - 基因组DNA的长度范围是什么? - RNA的长度范围是什么? - 测序的长度是什么?

因为测序的读长(人类和机器的能力)远低于DNA或RNA的长度(自然的能力),所以我们不得不采用曹冲称象的方法来还原基因组和转录组。

RNA-seq 分析流程

建测序文库是将RNA打断成小片段的过程。测序之后,得到一堆一堆的小石头(称为Reads,长度100-150bp)。

生物信息学分析是将小石头对应到基因表达上的过程,是建库的逆过程。通过比较各个基因石头堆中Reads的数目,可以得出基因表达差异。

RNA-seq的实验和分析流程

RNA-seq实验一般至少有两组样品,分别是对照组和实验组,每组一般设3个重复。

对照组 实验组 实验设计
wild type mutant of a transcriptional regulator 转录因子调控的靶基因
non-treated treated 处理引起(诱导)的基因表达变化情况
0 h 6h, 12h, etc 基因表达随时间的变化情况

实验设置重复对于发表论文是非常重要的。没有重复的结果很难被杂志和审稿人接受。

一个成功的RNA-seq实验,合理的实验设计是基本需要。获得实验中所需的样品可能要花很多精力。另外,由于RNA不稳定的特性,样本制备和保存也要特别注意。

测序一般由商业公司完成。剩下的主要工作就是生物信息学分析。

注意:这一部分先讲基本原理,下一部分讲生物信息学操作。

文库的构建:ssRNA-seq的介绍

目前主流的RNA-seq文库构建方法分为链特异性和非特异性文库。链特异性文库包含了RNA所在DNA链的信息,成本只比后者略高,因此建议构建链特异性文库(Strand-specific RNA-seq,ssRNA-seq)。

我们知道,DNA是双链,RNA是单链。所以同样一段DNA可以有两种序列不同的RNA产生。ssRNA-seq在建库的过程中对序列标记,保留了链特异性信息。

常见的链特异性文库构建方法是dUTP方法

dUTP

在合成第二链的时候,使用dUTP替代dTTP,合成之后通过UDG酶降解第二链,从而得到特异性的文库。

参考文献:

  • Zhong S et al. (2011). High-Throughput Illumina Strand-Specific RNA Sequencing Library Preparation. Cold Spring Harb Protoc 2011: pdb.prot5652.
  • Levin JZ et al. (2010). Comprehensive comparative analysis of strand-specific RNA sequencing methods. Nat Meth 7: 709–715.
  • Zhang Z et al. (2012). Strand-specific libraries for high throughput RNA sequencing (RNA-Seq) prepared without poly(A) selection. Silence 3: 9.

短序列与基因组的比对

我们前面提到,每个基因(大象)测序(称量)之后都对应一堆Reads(石头)。但问题是,这些石头是一大堆。如何将它们分开呢?

这个过程就是短序列比对

短序列比对需要将百万量级的Reads快速、准确的比对到基因组上。采用了不同于BLAST的算法、结果呈现方式等。

基因组注释和差异表达基因的鉴定

在上一步将reads比对到基因组上后,可以根据基因组注释信息确定对应基因上序列的个数。

Reads count

简单说来,会有下面的结果。

Gene Control Treated Fold Change
Gene1 3 1 13
Gene2 2 3 32
Gene3 1 1 11
Gene4 2 1 12

当然,事实上的情况比这要复杂。上图中,Control一共有7个Reads,而Treated一共有5个Reads;Gene4的长度也比其它要长;Gene3和Gene4还有一部分重叠。这些情况在真正的分析时都要考虑,会使用相应的模型和统计算法来完成。

好消息是,生物信息学软件将为你做好了这一切。所以,你不要太担心。

差异表达基因的富集分析

在上一步之后,我们将会得到差异表达的基因。这些基因可能会有很多,表达量变化幅度可能也不一样。比如,我们分析发现:200 个基因表达上调,150 个基因表达下调。这有什么生物学意义呢?

这就是富集分析要解决的问题。富集分析可以告诉你,基因表达变化主要发生在哪些代谢途径?哪个细胞结构(细胞壁、细胞膜、细胞器)?

假如基因组中一共有2000个基因,其中200个是与氨基酸代谢相关的(酶、调控蛋白等),占全部基因的十分之一。而恰好200个上调表达的基因里面,100个都是氨基酸代谢相关的基因。那意味着什么?

如果基因表达差异是随机的,那我随机选择200个基因,按照十分之一的几率,其中最有可能是有20个左右的氨基酸代谢相关基因。但是在200个上调表达的基因中,却有一半的基因是氨基酸代谢相关基因。说明氨基酸代谢途径是显著受影响的代谢途径,在差异表达基因基因中富集。这就是富集分析。

富集程度以超几何分布概率公式计算。其中N是基因总数,m是具某一属性(如氨基酸代谢相关)的基因总数(m < N)。n是总的差异表达基因,k是具同一属性(如氨基酸代谢相关)基因的个数。

$$f(k;n,m,N) = {{{m \choose k} {{N-m} \choose {n-k}}}\over {N \choose n}}$$

在上述的例子中,计算phyper(100,1800,200,200)的值为$3.0 \times 10^{-56}$。如果出现这样的结果,肯定不是一个巧合。

如何针对数据下一个结论

小概率事件不可能发生的,发生了小概率事件,那只能说明处理组中氨基酸代谢相关途径显著上调了。意味着细胞开始大量代谢氨基酸。是合成?还是降解?都可以用类似的富集分析方法去分析。

这就是RNA-seq的结论部分。

商业公司RNA-seq结果解析

测序公司都会出具一个结题报告,并将分析结果打包给你。

基本上,上面描述的分析,测序公司都会包含在内(只多不少)。

但是测序公司的结题报告,存在一些不足之处。

  • 千篇一律,有用的信息夹杂在无用的信息之间;
  • 做得好测序,却做不好数据分析和解读;基本上,越到结题报告的后面,越让人看不下去。
  • 公司的分析数据库常常缺乏更新;
  • 你不能设定分析的参数;
  • 甚至包含着错误的信息。

结题报告(有用、无用、错误和正确的信息)

下面是某大基因公司的结题报告。

下面是某诺公司的结题报告。

一般来讲,共有类似的文件夹结构。

文件夹的结构

--Results
 +- Rawdata    | CleanDataAssess
 +- Genome     
 +- Map        | MappingAssess
 +- Expression | QuantitativeAnalysis
 +- DESeq
 +- Structure  | StructureAnalysis
report.pdf
...

你需要什么内容,应该可以按图索骥了。

从头开始 RNA-seq 分析

将 raw reads 比对到参考基因组

使用的生物信息学工具/软件:bowtie2或hisat2或其它短序列比对软件。

Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. HISAT2 is a fast and sensitive alignment program for mapping next-generation sequencing reads (whole-genome, transcriptome, and exome sequencing data) against the general human population (as well as against a single reference genome).

测序结果是以 fastq格式分两个Read保存的,通常会压缩,扩展名推荐为 *.fastq.gz

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

现在的测序结果一般是双端测序,会有两个分开的文件。

短序列比对虽然与常用的BLAST比对差别较大,但是二者实际上做的是同样的工作:都是把一段较短的序列比对到数据库中较长的序列上。

首先需要建立一个基因组数据库索引,然后才能把测序得出的Paired-reads(fastq格式文件)比对到基因组上,比对结果输出为*.sam格式文件。

关于fastq格式,参考:FASTQ格式

关于*.sam格式,参考:SAM格式说明

## 从基因组序列 genome.fasta 创建名为 genome 的 hisat2 专用索引
hisat2-build genome.fasta genome

## 运行下列命令将 forward.fq.gz 和 reverse.fq.gz 中的测序结果比对到 genome 索引,文件输出为 file.sam
hisat2 -1 forward.fq.gz -2 reverse.fq.gz -x genome -S file.sam

file.samTAB分隔的文本文件,可以用less file.sam查看。

计算落到每个基因上面的 Reads 个数 - Reads count

使用的生物信息学软件有:samtools,htseq,NCBI数据库。

Samtools is a set of utilities that manipulate alignments in the BAM format. It imports from and exports to the SAM (Sequence Alignment/Map) format, does sorting, merging and indexing, and allows to retrieve reads in any regions swiftly.

HTSeq is a Python package that provides infrastructure to process data from high-throughput sequencing assays.

下图是一个计数模型。

count-model

## 格式转换 SAM -> BAM
samtools view -bS file.sam -o file.bam

## 将比对结果按照 read 的名称排序
samtools sort -n file.bam #sort bam by name

## 计算落到各个基因上面的 Reads 个数
htseq-count -f bam --stranded=yes --type=CDS --idattr=id file.bam gene.gff

靠谱的基因表达差异工具:DESeq2

能用于RNA-seq分析的工具很多。不同工具采用的算法和模型各异,导致结果存在差异。

2017年,Nature Communication上发表的一篇论文详细比较了RNA-seq各种工具的性能,最终得出:DESeq2是最靠谱的分析工具之一。

生信人公众号解读:必看|史上最全的39个RNA-seq分析工具

推荐使用!

靠谱的基因富集分析工具:ClusterProfiler

同样的,基因富集分析的工具也比较多,这里推荐ClusterProfiler,原因有三:

  • 支持人、动物、植物和微生物等各物种
  • 软件仍然还在维护中
  • 分析数据库实时更新
  • 内嵌了丰富的可视化工具,随便修改一下就是发表级

推荐使用!

作者的微信公众号:biobabble ClusterProfiler应用文章目录点这里