RNASeq数据分析

从原始的数据开始，进行reads回帖，到拼接转录本，计算表达量，分析差异表达，最后可视化分析结果。

TopHat是一个把reads回帖到基因组上的工具。首先用Bowtie把reads回帖到基因组上，然后通过拼接，我们就可以在基因组上看到一些reads堆叠起来的区域，称为consensus,这些consensus可能是一个真的外显子，也有可能是几个外显子拼在一起的，或者一些别的情况。我们知道，经典的剪切位点一般都有GT和AG这样的序列标志，在consensus的边界和内部，TopHat会去找这样的剪切位点，并且得到他们可能的组合。然后对于那些没有被Bowtie贴到基因组上的reads,TopHat会对他们建立索引，去和这些可能的剪切位点比对，这样就把跨越剪切位点的reads准确地贴到基因组上。

一些比较重要的命令行选项。

关于插入片段长度的选项：在RNA-Seq中，会把mRNA打断成小的片段，然后对片段长度进行iding筛选后拿去测序，如果选择的片段长度是300bp，两端各测序75bp的reads,中间的插入片段长度就应该设为150bp.

下面是设置插入片段长度的标准差，如果选择的片段长度比较集中，这个值可以设置的小一些，反之应该设置得大一些。

-G选项是提供哦呢一个已有的注释文件。如果你分析的基因组被注释得比较好了，最好能够提供这个文件， 这时TopHat就会先把reads往转录组上贴，没有贴到转录组上的再往基因组上贴，最后把结果合并起来。我们知道大多数的转录组都是比基因组小得多的，而且junction reads可以直接贴到转录本上，所以这样回帖的效力和准确度都可以

标准的Illumina平台是不分链的，我们无法知道配对的reads哪个方向和转录本一致，哪个和转录本反向互补。对于分链的数据，也有两种情况，在firststrand这种分链方法中，第二个read和转录本方向一致，第一个read和转录本反向互补，在另一种fr-secondstrand分链方法中，就刚好反过来了。所以在分析的时候一定要弄清楚自己的数据有没有分链，是怎么分链的。

下面是一个模拟的RNA-Seq数据集，双端测序，有两种处理，每种处理有3个重复，这里C代表处理，R代表重复，下面用C1R1进行演示

首先，要有参考序列fasta文件，也就通常说的基因组序列。 TopHat是利用Bowtie2回帖reads,我们首先需要建立Bowtie2的索引文件： bowtie2-build genome.fa(基因组文件) genome （注意程序和文件所在目录）

我们还需要reads的fastq文件，双端测序的数据，两个fastq文件分别以下划线1和2这样的形式结尾。在实际分析中，需要对拿到的数据进行质量评估和过滤等依稀类预处理工作，这些工作都是非常重要的。需要准备注释文件，当然它不是必须的。它可以是GTF或者GFF3格式的文件，对于注释得比较好的基因组，在UCSC可以下载。

准备好后就可以运行tophat了，-p是线程数，-G是注释文件，-o是输出文件夹，选项之后就是参考序列的索引，最后是两个reads的fastq文件。

看里面生成的文件，align-summary文件爱你，这个文件是reads回帖的一些统计信息。90%以上的回帖比例就非常好了，当然百分之70以上一个可以接受的范围。bam文件详细记录了reads回帖到基因组上的情况，由于这是一个二进制的文件，我们需要用samtools查看它。

Cufflinks是一套拼接转录本，计算表达量，计算差异表达的工具。尽可能拼接处最优可能的转录本的结构，并且估计它的表达量。

-G是提供一个注释文件，并且告诉Cufflinks不要去拼接新的转录本，只能用注释文件里提供的转录本。

-g 也是提供一个注释文件，但是Cufflinks会在这些已知转录本的指导下，拼接新的转录本。

-u / multi-read-count

-u是告诉Cufflinks用更准确的方法去处理贴到多个位点上的reads，如果没有-u,Cufflinks只会把这些reads简单地平均分配。比如一个read贴到了10个位置，那么每个位置分得十分之一。加-u后会先进行平均分配，然后按照这10个位置各自的表达量，计算read被分配到每个位置的概率。实际上Cufflinks会用EM算法进行迭代，计算在观察到当前数据的情况下，最优可能的reads分配。 library type和TopHat里面差不多。

这里的bam文件就是刚才TopHat运行的结果。

Cuffmerge

当我们使用Cufflinks处理多个数据之后，我们需要将其转录本数据整合为一个全面的转录本集合，Cuffmerge是一个将Cufflinks生成的gtf文件融合为一个更加全面的转录本注释结果的工具。如下图所示，图中的6个转录本被整合为一个转录本集合。同时我们可以利用基因组注释文件，获得更加准确可靠的结果。合并后的转录本集合为计算每个基因和转录本的表达量提供了一个统一的基础。

-g参数 指向参考GTF文件。

-p参数决定线程数

-s参考指向基因组DNA序列。如果是一个文件夹，，每个contig是一个fasta文件。如果是一个fasta文 件，则所有的contigs都需要在里面。

最后一项是一个列表，内容包括经过Cufflinks拼接的转录本的文件路径。

流程如下：

首先，我们需要用cat命令创建一个所有拼接的转录本的文件路径列表。

cat> assemblies.txt

复制上图的列表 然后运行cuffmerge,运行后的结果存储在merge_asm这个文件夹里面。其文件夹内包含一个Logs文件夹以及一个.gtf文件，也就是我们经过整合的转录本文件。

cuffdiff

当我们利用Cufflinks获得了拼接的转录本时，我们就可以计算不同样品中转录本的表达量。计算的简单原理在于测序深度和外显子长度一定时，Read的数量与对应的转录本数量成正比。通过对Reads进行计数计算转录本的表达量。同时cuffdiff可以计算不同条件下转录本表达水平的显著性差异。

-u /--multi-read-count # default:FALSE
-L /--labels

-u命令指cuffdiff对回帖的基因组中多个位置的read进行一个初步的估计，然后加权分配到各个基因组位置。而不是简单的平均分配，其功能与Cufflinks中的u命令相同。

-L为每个样品标上名称。

接下来Cuffmerge产生的gtf文件，Cuffdiff需要它提供的注释进行初始转录产物和可变剪切等定量分析。最后是TopHat产生的bam文件，如果一个样品中有多个实验重复，那么我们需要提供哦呢bam文件列表，文件名之间以逗号隔开。

运行之后，cuffdiff输出的文件在diff_out目录之下。其中包括一些按类别统计的表达水平结果，如果有相同的转录起始位点，或具有相同的编码区的转录本的表达水平，我么你可以利用他们进行下一步的分析。

CummeRbund

当我们对数据进行分析之后，我们希望将结果可视化，CummeRbund作为一个分析结果数据的软件包，在作图方面具有出色的能力。我们也可以在R环境下输入以下命令进行下载和安装。

常见的绘图方式有：密度分布图，闪点图，火山图，柱形图

这些命令的输入都是cuffdiff的输出文件，经过处理后的结果，我们会在之后进行详细的介绍。这是为不同条件的样本标记名称。

参考资料：

ps:听了这么久的课，总算讲了软件的用法，收获颇多，颇多。

山东大学 生物信息学课题组荣誉出品 http://www.crc.sdu.edu.cn/bioinfo 巩晶老师课件