本文是我在华中农业大学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](/post/2019-07-17-RNA-seq_files/RNA-seq-reads count.png)

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

Gene Control Treated Fold Change
Gene1 3 1 1/3
Gene2 2 3 3/2
Gene3 1 1 1/1
Gene4 2 1 1/2

当然,事实上的情况比这要复杂。上图中,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应用文章目录点这里

作者简介

Chun-Hui Gao is a Research Associate at Huazhong Agricultural University.

重复使用

Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The source code is licensed under MIT. The full source is available at https://github.com/yihui/hugo-prose.

欢迎修订

如果您发现本文里含有任何错误(包括错别字和标点符号),欢迎在本站的 GitHub 项目里提交修订意见。

引用本文

如果您使用了本文的内容,请按照以下方式引用:

gaoch (2019). 原核转录组分析的原理和流程. BIO-SPRING. https://bio-spring.top/post/2019/07/17/rna-seq/

BibTeX citation

@misc{
  title = "原核转录组分析的原理和流程",
  author = "gaoch",
  year = "2019",
  journal = "BIO-SPRING",
  note = "https://bio-spring.top/post/2019/07/17/rna-seq/"
}