《生物信息学：导论与方法》--新一代测序NGS：转录组分析RNA-Seq--听课笔记（十四）

第八章 新一代测序NGS：转录组分析RNA-Seq

8.1 转录组介绍

• A transcriptome is a collection of all the transcripts present in a given cell.
• 所谓转录组（transcriptome），是指特定细胞中全体转录本（transcript）的集合。
• 也可以说是细胞特定时刻基因表达谱的一个快照（snapshot of expression profile）
• 在转录组中，既包括早已熟悉的、编码蛋白的信使RNA（mRNA），也包括近年来新发现的不编码蛋白的miRNA、long non-coding RNA(lncRNA)等非编码RNA。
• 转录组中的RNA转录本彼此协同作用，共同来调控细胞生长、发育、凋亡等一系列重要的生理过程。
• 对转录组的研究包括定性（qualitative）和定量（quantitative）。
• qualitative(鉴定出所有表达的转录本)：identify all transcrits, i.e. expressed genes as well as their isoforms(亚型)
• quantitative(确定这些转录本各自的表达量)：estimate the expression level of these transcripts, i.e. the transcript abundance of expressed genes/isoforms
• 实时荧光定量PCR（Real Time quantitative Reverse Transcription PCR, Real-Time qRT-PCR）：通过对经典PCR扩增反应中每一个循环产物荧光信号的实时检测，可以实现对起始模板的定量分析。
• 通过正确设计引物（primer）和探针（probe），qRT-PCR技术可以在很大范围内定量的检测目标转录本的拷贝数，也即表达水平，因此常被作为转录组分析中的金标准（Gold Standard）。
• qRT-PCR技术一次只能测定一个转录本的表达水平，同时，这个技术需要事先知道待检测转录本的序列，难以用来发现未知的转录本。Low-throughput; Prior knowledge of transcript sequences needed.
• 在深度测序获得大规模应用前，Microarray是主要的高通量转录组表达分析技术。
• 微阵列（microarray），也称基因芯片（gene chip），是通过将几十万个不等的探针（probe）分子固定在约1cm见方的固体片基上制成的。利用核苷酸分子在形成双链时碱基互补原理，microarray可以一次性检测出样本中所有与探针互补的核苷酸片段，从而快速得到样本中基因的表达谱（expression profile）。
• 与qRT-PCR技术相比，microarray虽然在通量上有了显著的提高，但仍然需要事先确定待检测转录本的序列。
• 表达序列签（Expressed Sequence Tag, EST）技术通过对一个随机选择的cDNA克隆进行单次测序来获得cDNA的部分序列。
• 与microarray不同，EST是基于测序的，并不需要事先知道待检测转录本的序列，因此可以被用来发现新的转录本。
• 由于当时测序技术通量的限制，一次EST通常只能得到几千个转录本的序列，远远无法进行全转录组水平的profiling。
• RNA-Seq技术，即深度测序技术的出现，使得研究人员首次可以在全转录组水平利用测序技术同时进行定量与定性分析。
• RNA-Seq技术首先对生物样品中的RNA反转录为cDNA，而后，将这些cDNA打碎为较小片段，上机进行测序。
• RNA-Seq技术使得研究人员可以快速确定转录组，进而鉴定存在的可变剪切体（Alternative splicing isoform）；另一方面，通过对基因组特定位点上reads深度的计数，可以对表达量进行估计。
• RNA-Seq技术本质上是对转录本序列的随机抽样，因此，其检测效力（power）和灵敏度（sensitivity）高度依赖于测序的深度。如果深度不够，就会难以检出低拷贝的基因。原则上，只有在饱和曲线（saturation curve）达到平台期（plateau）后，才能认为深度足够。
• 在随机抽样的情况下，可以认为map到特定转录本上的read数目正比于起表达量（transcript abundance），因此我们可以利用落在某个转录本上的reads的总数来估计其表达量。另一方面，落在一个转录本上的reads的数目，也与其长度和总测序深度成正比。
• 针对上述问题，在实际分析中，通常会将原始的reads数目利用线性放缩，转换为RPKM值来进行归一化处理。
•     

1. C：the number of mapped reads for specified transcript (贴到这段转录本上的reads总数目)
2. N：the number of total mapped reads（这次实验中总的reads数目，即测序深度）
3. L：the length of the specified transcript（这段序列的长度）

• 在假定不同样本中RNA的总体分布一致的前提下，RPKM可以正确处理由于转录本长度和测序深度引起的artifact，从而使得来自不同基因、不同sequencing run乃至不同样本之间的表达数据彼此之间可以比较。
• RPKM并非是唯一的归一化方法。
• 通过考虑不同的误差因素（bias effectors）、引入不同的生物学假设，可以构造不同的归一化方法。
• 目前已有研究表明，相比于后续提出的TMM、DESeq等方法，RPKM方法在样本间差异基因表达检验等分析中的效果并不是最理想。
• 另一个需要在RNA-Seq技术中引入注意的地方是链特异性（strand-specific）。DNA的两条链都可以转录，形成不同的转录本。然而常用的illumina RNA-Seq kit是不分链的，也就是说我们无法知道配对的reads哪个方向和转录本一致，哪个和转录本反向互补。
• 在分析之前一定要弄清楚数据有没有分链，是怎样分链的。

8.2 RNA测序数据回帖与组装

• 同样是基于深度测序技术，因此在reads长度、数量、质量等方面，RNA-Seq产生的reads与之前基因组重测序产生的DNA reads具有相似的性质：长度短，数量大，质量参差不齐，错误率较高。
• 在数据上，由于RNA-Seq测序数据来源于RNA转录本，特别的，在DNA转录成mRNA的过程中，内含子被切掉，外显子会在剪切位点连接到一起，对于这些跨过剪切位点的reads，也就是所谓的junction reads，如果不把它们从中间断开，就无法准确贴到基因组上。
• junction reads是确定剪切位点的直接证据，对于正确重构转录本结构至关重要。
• 算法在mapping时，需要“意识到”junction site和intron的存在，从而正确处理这些junction reads。
• Handle Junction Reads："Join exon" strategy

1. 基于已知转录本中所有exon来构建所有可能的junctions的library，而这个library未必是已知的，而是所有可能的组合。
2. 采取之前DNA reads类似的unspliced方式基因组比对，将非junction reads map到基因组，对于无法直接map的junction reads，则再与第一步中构建的junction library进行比对。

• Join exon策略实质上相当于对之前DNA reads mapping算法的一个“补丁” （patch），通过构建所有可能的junction库，该策略可以发现新的splicing isoform。然而对于以前未知的exon，这个策略就无能为力了。？
• 为了克服上述问题，可以使用“Split reads” strategy。
• 采取之前DNA reads类似的unspliced方式基因组比对，将非junction reads map到基因组，对于无法直接map的junction reads，则参照之前BLAST的方法，切分成若干个长度为k的种子（seed），而后再利用这些种子，重试mapping，即在更小的粒度上尝试寻找junction site；最后再将临近的mapped seeds重新组合起来，来得到最终全read的alignment。
• Split reads策略相对Join exon策略速度更慢；但这个策略不依赖于先验的exon注释，可以发现新的exon乃至新的基因。
• 目前常用的RNA-Seq工具通常会组合两个策略，以平衡检测的灵敏度与速度。
• TopHat2工具：首先利用join exon策略快速检测已知的junction site，然后再利用split reads策略来发现新的junctions，针对不同的阶段采用了不同的索引（index），可以进一步提高mapping的速度。
• 完成mapping后，还需进行：

1. Assembly: reconstruct full-length transcript sequences from the (mapped) reads.
2. Quantification: estimate the expression abundance for each transcripts

• 在正确mapping后，可以将转录本组装问题描述为一个对有向图（directed graph）的遍历问题。通过对图中的边给予不同的权值作为约束，可以利用图论（graph theory）中的寻路算法（path finding algorithm）找到一条或多条符合约束的最优路径及其对应的转录本序列。
• Cufflinks工具：一个针对RNA-Seq数据进行转录本组装和表达分析的工具软件。
• 实际中，为了更为准确的估计表达量，常常会采用EM等方法来反复迭代的考虑转录本组装与表达量估计。

8.3 RNA-Seq数据分析（TopHat 和 cufflinks工具介绍）

• TopHat、Cufflinks、Cuffmerge、Cuffdiff、CummeRbund（一个R包）
• RNA-Seq分析流程：从raw data开始，进行reads回帖，到拼接转录本，计算表达量，分析差异表达，最后可视化分析结果。
• TopHat: 把reads mapping到基因组的工具，A spliced read mapper for RNA-Seq
• 经典的剪切位点一般都有GT和AG这样的序列标志。
• 
• TopHat----Important options：

1. -r/ --mate-inner-dist <int> : 插入片段长度的选项。
2. --mate-std-dev <int>: 设置插入片段长度的标准差，如果选择的片段长度比较集中，这个值可以设置得小一些，反之应该设置得大一些。
3. -G/ --GTF <GTF/GFF3 file>：提供一个已有的注释文件。如果分析的基因组被注释得比较好了，最好能够提供这个文件。
4. 设置了-G，TopHat就会先把reads往转录组上贴，剩下的再往基因组上贴，最后把结果合并起来。大多数转录组都是比基因组小得多的，而且junction reads可以直接贴到转录本上，所以这样mapping的效率和准确度都会得到提高。
5. --library-type

• 演示数据集：GSE32038
• 实际分析过程中，需要对拿到的数据进行质量评估和过滤等一系列预处理工作，这些工作非常重要。
• 注释得比较好的基因组的注释文件，可以在UCSC下载。
• tophat -p 16 (线程数) -G genes.gtf  -o C1_R1_thout genome  GSM794486_C1_R1_1.fq GSM794486_C1_R1_2.fq
• 真实数据中，90%以上的回帖比例结果就非常好了，60%~70%也是一个可以接受的范围。
• 二进制文件用samtools查看。
• Cufflinks是一套列拼接转录本、计算表达量、计算差异表达的工具。
• Cufflinks就是尽可能地拼接处最有可能的转录本结构，并且估计它的表达量。
• Cufflinks----Important options：

1. -G/--GTF <reference_annotation.(gtf/gff)>：Use the supplied annotation and not assemble novel transcripts, 告诉Cufflinks不要去拼接新的转录本，只能用注释文件里提供的转录本。
2. -g/--GTF-guide <reference_annotation.(gtf/gff)>：Use the supplied annotation to guide the assembly of novel transcripts. 提供一个注释文件，但是Cufflinks会在这些已知转录本的指导下，拼接新的转录本。
3. -u/--multi-read-correct：告诉Cufflinks用更准确的方法去处理贴到多个位点上的reads，如果没有-u，Cufflinks只会把这些reads简单地平均分配。比如一个read贴到了10个位置，那么每个位置分得十分之一。有-u，则比如说会先进行平均分配，然后按照这10个位置各自的表达量，计算read被分配到每个位置的概率。会用EM算法进行迭代，计算在观察当前数据的情况下，最有可能的reads分配。
4. --library-type

• cufflinks -p 16 (线程数) -g genes.gtf  -o C1_R1_clout  C1_R1_thout/accepted_hits.bam(Mapped Bam/Sam file)
• 当我们使用Cufflinks处理多个数据之后，我们需要将转录本数据整合为一个全面的转录本集合。
• Cuffmerge是一个将Cufflinks生成的gtf文件融合成一个更加全面的转录本注释结果的工具。同时，还可以利用基因组注释文件，获得更加准确可靠的结果。
• 合并后的转录本集合为计算每个基因和转录本的表达量提供了一个统一的基础。
• Cuffmerge----Important options：

1. -g/ --ref-gtf                #optional
2. -p/--num-threads      #default:1
3. -s/--ref-sequence<seq_dir>/<seq_fastq>

• cuffmerge ‐g genes.gtf ‐s genome.fa ‐p 16 ‐o cuffmerge_out assemblies.txt
• 当利用Cufflinks获得了拼接的转录本时，我们可以计算不同样品中转录本的表达量。计算的简单原理在于测序深度和外显子长度一定时，read的数量与对应的转录本数量成正比。通过对reads进行计数计算转录本的表达量。同时Cuffdiff可以计算不同条件下转录本表达水平的显著性差异。
• Cuffdiff----Important options：

1. -u/--multi-read-correct     #default:FALSE 指Cuffdiff对回帖到基因组中多个位置的read进行一个初步估计，然后加权分配到各个基因组位置，而不是简单的平均分配。与Cufflinks中的u命令相同。
2. -L/--labels                        #为每个样本标上名称

• 
• CummeRbund：A R package to aid and simplify the task of analyzing Cufflinks RNA-Seq output.
• CummeRbund----Important options：

1. > csDensity(genes(cuff_data))
2. > csScatter(genes(cuff_data), 'C1', 'C2')
3. > csVolcano(genes(cuff_data), 'C1', 'C2')
4. > expressionBarplot(mygene)
5. > expressionBarplot(isoforms(mygene))

• 密度分布图可以表示不同表达水平的转录本的密度分布。
• 散点图可表示两种条件下转录表达水平的情况，位于直线两侧的点对应的转录本其表达水平具有不同条件偏好。
• 火山图可表示不同条件下基因表达是否具有显著性差异，横坐标指不同条件下基因表达水平之比的对数值，纵坐标为T检验中p值的对数值。高于一定阈值的点对应的基因存在显著差异性表达。然后可以用getGene命令对特定的基因进行代入分析。