19 tophat(二) 具体演示（重中之重来了）

A为对照组，B为处理组，A和B都采用illumina测序，reads读长分别为70bp，插入片段大小为200bp，reads之间的距离为200-70-70=60bp，将这个两个样品分别与人全基因组参考序列进行比对，也就是前面下载好的参考序列文件和GTF文件。

（因为数据量比较大，所以两个样本作者只截取500M数据进行展示）

后面会描述如何评估测序数据是否达到饱和。**

在做数据分析时，目录结构一定要规范，包括文件命名，否则样品多的时候或者文件大的时候就会非常混乱，要养成简洁高效的习惯

[if !supportLists]1. [endif]首先需要使用bowtie2_build对参考文件建立索引：

[if !supportLists]2. [endif]正式比对过程：cat tophat.sh     tophat -o A -G ~/Data/ref/hg38.gtf -r 60 --library-type frunstranded 2 -p 3 ~/Data/ref/hg38 ~/Data/reads/A.1.fq.gz （若有多个文件则用,号分割开）（优先使用-G，GTF文件；注意与参考序列要对应，里面文件需要一致。）

[if !supportLists]3. [endif]比对结果在输出目录中，tophat可以直接输出bam格式的结果，在数据结果的文件中会生成很多文件。

[if !supportLists]4. [endif]可以用samtools view 查看一下，通常只用1对1的比对情况，可以提取出来

samtools view accepted_hist.bam |grep “NH:1:1$”|le

Spliced比对:外显子的可变剪切问题，来自一个外显子的reads可以比对到多个转录本是正常的情况

[if !supportLists]5. [endif]除了Bam文件，tophat会产生三个bed文件

[if !supportLists]6. [endif]比对产生的结果，在align_summary.txt文件中，可以通过cat align_summary.txt进行查看，文件中会列出详细的比对信息，比如多少大于相应的比对次数等

[if !supportLists]7. [endif]比较bowtie与tophat那种能够比对到更多的结果，理论上tophat能够比对到更多的reads。因为前一种是unspliced（不剪切）的方法，后一种是进行剪切spliced的方法。王老师得出的结果----实际与理论相符合。

20 RNAseq数据评估

（用完tophat比对完成之后，其实已经完成RNAseq中最核心的共作）---质控评估下

接下来，基于比对结果进行更多操作。操作之前，还需要将比对完的结果进行评估，而非直接进行基因表达计算；差异表达分析以及差异基因功能注释（可能存在异常，那会浪费大量时间）

[if !supportLists]1. [endif]评估主要两个指标：

（1测序饱和度：可以用以推测测序所要的数据量，其实就是量足够了（比如说细胞RNA1ug量就够了）理想是所有转录本都能测到，（定性）；每个转录本的丰度信息都能测量到（定量）。若不饱和----测序饱和度图 横轴：reads条数  纵轴：检测到的基因表达数

如何评估饱和度呢？---统计tophat比对完的bam文件中reads比对到哪个基因上，根据GTF上的信息可以筛选到。然后依次取50万条reads，检测bam文件中reads比对到多少个基因上，然后再取100万进行统计，最终就可以获得一张饱和度的评估图，可以直观地找到数据饱和地临界点。

（2 测序随机性：用物理或化学的方法将转录本打断成小片段，上机测序。若打断随机性差，reads偏向到基因的特定区域，则会影响转录本的各项分析结果。理论上，希望有多长则测多长，但是当前二代测不了，三代测的了，但是丰度一般不够。二代测，采用小片段文库，测两端。随机性主要是看基因是否是随机覆盖，而非只测序两端。

利用reads在基因上的分布来评价随机性，但由于不同基因有不同长度，把reads在基因上的位置标准化到相对位置，即read在基因组上的位置与长度的比值，然后统计不同基因比对上的reads数。如果打断随机性好，reads在基因上各部分应该分布比较均匀。

21序列比对FAQ（分析RNAseq序列比对中的常见问题）

[if !supportLists]1. [endif]关于RNAseq测序数据量（一次测多少，才有意义）

在测序数据前，可以对物种转录组进行评估，预估一个数据量

（1有参考序列的物种，可以分析基因组信息 （2 无参：则查看相近物种的转录组大小

区别于DNA测序数据量，如人基因组数据是3G，要测序10倍以上，即30G数据才能完全覆盖。RNAseq测的是转录本，由于真核生物基因组上绝大部分是重复区，转录的区域并不多，因此RNAseq不需要测那么多，节约了测序成本。人1%为基因区则300M左右，但

一般真核生物推荐4-6G测序转录组数据量，4G够了但是不保险，6G基本完全够了，10倍去测的3G是肯定不够的。

[if !supportLists]2. [endif]RNAseq测序不饱和有哪些影响

（1从定性上来说，会导致一些低丰度转录本测序不到，不能反映真实情况

（2从定量上来说，由于没有饱和，每个基因的表达量统计会不准确，最终得到的差异表达也会有问题。

[if !supportLists]3. [endif]可以比对到基因组，但是比对不到基因基

（1参考序列亲缘关系比较远，导致比对率低（都低）

（2由于rRNA去除不干净,rRNA比对率高 ---mRNA少（多见于原核）

（3 存在DNA污染（会出现对全基因组覆盖度非常高）

[if !supportLists]4. [endif]基因比对率低会有哪些影响

（1如果核糖体去除不干净，分析结果不可靠

（2参考序列不近源的原因，很多基因无法进行定量

（3 unspliced比对方法的影响(一定要选择合适的比对策略)

22基因表达定量

[if !supportLists]1. [endif]表达量计算

根据比对结果，统计出每个基因（转录本）对应reads，在对其进行标准化处理

[if !supportLists]2. [endif]不同基因测序reads数（可以根据reads数测序基因表达量）

[if !supportLists]3. [endif]测序中的两种情况（1不同测序长度（同一个基因测的reads数一样，理论上短基因长度的基因表达量更高，因为其还有一大部分基因区域测到的reads数还没有进行统计）（2.不同测序深度，一个测序50倍，一个测序100倍得到的A基因reads数相同，测的少的A基因表达更高（测序量差异，如一个测了200万，一个只测了100万））即使reads数相同，但基因表达量不同

[if !supportLists]4. [endif]不能单独以reads表达量，确定基因表达量的差异。要计算基因的表达量与reads条数，基因长度以及基因深度有关系。综合上述三个因素，Mortazavi得出基因表达量计算公式RPKM（Reads per Kb per Million reads）法 （RPKM：每百万reads中每一千bp覆盖的reads数） 每百万reads固定基因测序量（测序深度），每一千bp:基因长度就固定了

RPKM=106C/(NL/103) 设RPKM(A)为基因A的表达量，则C为唯一比对到基因A的reads数，N为唯一比对到参考基因的总reads数，L为给予你A编码区的碱基数。RPKM法能消除基因长度和测序量差异对计算基因表达的影响。计算得到的基因表达量可用于不同样品间的基因表达差异。

如果一个基因存在多个转录本，则用该基因的最长转录本计算其测序覆盖度和表达量。

[if !supportLists]5. [endif]RPKM一般只适合于原核生物不适合真核生物。

因为真核发生了可变剪切，可变剪切reads无法比对回去,相应的定量就比较少

[if !supportLists]6. [endif]FPKM表示，每一百万个map的reads到外显子的每一千个碱基上的reads数，FPKM与RPKM计算方法基本一致。注意FPKM计算的是片段fragments，而RPKM计算的是reads。fragment比reads的含义更广，因此FPKM包含的意义也更广。

FPKM= total exon fragments/(mapped reads(Millions)*exon length(KB))

思考：片段与read的区别，为什么FPKM就能用来对具有可变剪切的真核生物进行定量呢？？

23 rpkmforgenes.py----计算rpkm值的程序（使用python开发的程序）

无需安装就能使用，但是注意需要有python相关的两个模块

利用公式计算rpkm

Python rpkmforgenes.py（看看其参数，很多都是输入输出相关的选项，也含有与基因模型相关的选项）

介绍与注释文件相关的几个选项（因为bam文件是与基因组比对生成的，所以就需要一个注释文件指明哪部分是基因，哪部分是转录本，这就需要一个注释文件；可以是GTF,BED,reference.txt）

均一化定量选项（选择外显子全或者全转录组啥的进行定量）

接着作者介绍了具体的操作命令与输出结果组成（每行是什么）

命令: python rpkmforgene.py -i~///_hist.bam-orpkm.txt-a ~///hg38.bed -fulltranscript

运行软件输入文件（tophat生成的bam文件）输出结果 参考基因组  可选参数输出结果每一列含义（第一行为样本名第二行为有多少reads比对上了；第三行为用于表达定量的reads数第四行为软件运行的选项参数）接下来每一列（与输出选项有关；常见reads序列；位置；RPKM值）

24基因差异表达筛选

[if !supportLists]1. [endif]计算每个基因的RPKM值（量化值）

[if !supportLists]2. [endif]相差多少才能说明具有差异呢？

需要考虑基因建库随机性的影响，并不能单纯比较RPKM与FPKM的倍数来判断基因表达差异是否显著，就像不能根据基因reads数量判断表达量一样，需要综合考虑所有的影响因素。在利用RNAseq数据进行比较分析时，不同样本之间同一个基因是否存在显著表达时，一般需要进行统计，选取两个标准，一个是FoldChange，另外一个是通过FDR（多重检验校正）校正。不能只用一个FoldChange,还需要通过FDR进行较正。那么如何较正呢？将p value校正得到一个q value值，p value越小越好 。

最终是FDR越小，FoldChange越大，则差异越大。具体q value设计越小越好，尽量不要超过0.05，否则没有意义。

主要是FoldChange q value FDR三个值筛选出差异表达基因

思考这个统计学的计算过程

25 cufflinks（一）-----后续主要讲新转录本的寻找与定量等

[if !supportLists]1. [endif]前面说明了如何计算基因表达量RPKM与FPKM与差异表达基因,这些工作看起来非常复杂，实际并不困难，因为目前有大量工具可以使用。直接Tophat+cufflinks（认真阅读这篇文章Nat .Protocol），作者介绍了整体过程

[if !supportLists]2. [endif]cufflinks主要包括4个东西

[if !supportLists]3. [endif]软件安装（建议下载编译好的版本，直接下载解压缩，否则需要先下载相应的许多依赖）

[if !supportLists]4. [endif]-G是用GTF文件，-g是利用GTF文件指导，进行转录本重构，将没有比对到参考基因组的reads进行重新拼接，拼接出新的转录本，并对新的转录本进行表达量计算。----寻找新的转录本按照 -g（小g）选项

26cufflinks（二）

[if !supportLists]1. [endif]输入格式是sam或者bam格式，默认为tophat处理好的sam或者bam格式。利用-G进行定量，-g进行转录本的重构，寻找新的转录本。另外还需要一个gtf文件

[if !supportLists]2. [endif]cufflinks -G ///（-G进行定量）

[if !supportLists]3. [endif]cufflinks -g///（-g进行转录本定量，由于-G是只需要进行表达定量统计表达数时间短，-g需要重新拼接时间需要较长）

[if !supportLists]4. [endif]采用已经运行的结果进行展示，每个结果都包含4个结构 le transcripts.gtf查看具体结果文件组成， :号，再输入q就能退出。（含有CUFF开头的，就是拼接出来的新转录本，其中拼接出来的转录本又称transfrag，transcript与fragment的合成）还有一个skipped的文件代表过滤掉的基因（skipped又取决于所设置的参数）

[if !supportLists]5. [endif]一定要区别内显子；外显子；isoform，基因等概念，一个基因对映多个isoform，FPKM_corf_hi上拼接，FPKM_corf_lo下拼接（用专门的计算公式）

27cuffmerge（将各个cufflink生成的transcript.ztf文件融合（打包）成一个更加全面的transcript注释文件，merge.ztf，以利于分析）---可以一下解决多个拼接的转录本。】

[if !supportLists]1. [endif]cuffmerge命令会直接得到其相应的选项参数信息。

[if !supportLists]2. [endif]cat assemlies.txt(将两个样本路劲写入)

[if !supportLists]3. [endif]cuffmerge -g ////...（具体命令待后续操作）

[if !supportLists]4. [endif]最终会生成一个merge.gtf文件。(既有原有转录本信息，又有重新拼接的)，内容更全，将这个gtf重新使用tophat工具重新进行一次短序列比对，得到新的bam文件，这样得到的文件将会更加准确。

28cuffcompare(使用cuffmerge的gtf结果，将其结果进行比较，主要用于新转录本的查找)

[if !supportLists]1. [endif]使用顺序，先用一个cufflinks对每个样本进行转录本重构，再使用cuffmerge合并成一个更全的转录本，非冗余的并集，得到merge.gtf。这时使用cuffcompare将merged.gtf与参考序列的gtf文件进行比较，就可以找到二者之间共有的部分与转录本部分。

[if !supportLists]2. [endif]介绍cuffcompare使用，同样输入cuffcompare就会显示其相应选项参数

[if !supportLists]3. [endif]王老师喜欢使用ll查看

[if !supportLists]4. [endif]介绍运行参数、运行命令与运行结果，最终也是生成gtf文件（cuffmerge更关注得到一个更全的转录本，而cuffcompare主要是为了得到与已知GTF有差异的东西）

29 cuffdiff（一）

（发现转录本表达，剪接，启动子所用的情况的明显变化。关注的差异表达就是通过cuffdiff工具直接得到的）

[if !supportLists]1. [endif]cuffdiff首先进行基因表达量FPKM值的计算，然后通过FPKM计算Foldchange，然后进行pvalue检验，再利用FDR校正，筛选出显著差异表达的基因，这些过程，cuffdiff一步就能完成，只需讲tophat比对的结果直接输入给cuffdiff即可。下面看软件参数，其与cufflinks存在很多相同之处。

[if !supportLists]2. [endif]cuffdiff显示相应参数。

[if !supportLists]3. [endif]--nodiff不进行差异分析

30cuffdiff（二）（）

具体案例展示（使用tohat输出的A与B的bam文件进行差异基因筛选）

介绍执行命令

分析执行结果（文件组成，关键结果文件的内容组成（行与列））

gene_exp.diff是最终需要的（介绍其格式，可以用excel文件打开）

Count_tracking文件评估来自哪个转录本

[if !supportLists]31. [endif]cummRbund(将差异结果进行可视化，cummRbund是针对cufflinks RNAseq输出结果分析与可视化发出的R包，目的简化分析，其会产生一个数据库，将这些数据建立索引，容易查询与使用，提供多种函数，满足一般需求需要)

Cufflinks与cummRbund是同一个作者，因此软件相互之间配合非常好，无需进行额外的格式调整。可以直接以cuffdiff的运行结果作为输入，只需使用几个简单的命令就可以输出文献能用的图，包括各种火山图，密度分布图，热图等。

属于bioconduct

sourcr(“http://bioconductor.org/bioclite.R”)

bioclite(“cummRbund”)安装

？？cummRbund（两个问号查看安装了没，安装了有帮助文档）

library (cummRbund)

八幅图的绘制以及其图的含义（分布图；密度分布图；箱线图；包含生物学重复的箱线图；四分位数盒型图；散点图（可以比较生物学重复的差异性与样品间的相关系数散点图，越接近于一表达相似性就越高）火山图（差异基因分布情况）；两两之间比较的火山图）

diffData()函数筛选差异表达基因

32 基因表达定量FAQ(问题)

cufflinks分析流程（8步）

1.为什么不直接从tophat到cuffdiff

Cufflinks分析策略一（3步，1.tophat比对 7.cuffdiff筛选差异基因 8.cummeRbund对cuffdiff结果进行可视化处理）

Cufflinks分析策略二（134678）多了一个cufflinks重构cuffmerge重新比对的过程（用到的这个gtf文件更加完整）

两种策略的选择（如果测序样本量多，测序深度不深merge之后误差大，人基因小鼠的注释文件信息比较全了，都比较适合策略一，直接用cuffdiff）

[if !supportLists]3. [endif]基因表达量为零，如何进行差异表达计算（将0设置为非常小的值如0.001，就能进行差异分析计算）

[if !supportLists]4. [endif]如何提取差异表达的基因----gene_exp.diff；可以通过grep文件抓取出来

[if !supportLists]5. [endif]关于生物学重复实验---1.证明实验操作可重复 2.证明差异分析所需的差异基因是所需要的（相关性检查---散点图）

[if !supportLists]6. [endif]cuffllinks软件使用中的注意事项5点（1.cufflinks有算法，自身可以提供相关参数2.cufflinks使用加权分布更加精确而非平均分布；3.细菌不推荐使用cufflinks；4.cuffdiff所有基因和转录本FPKM=0，可能GTF中的染色体名字和BAM里的名字不匹配；5）

33 新转录本分析

[if !supportLists]1. [endif]新转录本定义（相较于原有的转录本，也就是已知的转录本）

[if !supportLists]2. [endif]为什么会有新转录本（1原有基因预测不准确 2.可变剪切 3.非基因区转录，如lnRNA）许多新转录本是肿瘤发生的重要标志

[if !supportLists]3. [endif]新转录本分析（具体步骤3步-----，也就是cufflinks中所使用的）

[if !supportLists]4. [endif]转录本重构策略比较（1.每个样品分别构建，merge合并 2.所有样品数据合并，然后构建。一般建议第一种，第二种太浪费资源与内存了）

34lncRNA分析（lncRNA>200nt）

1.lncRNA(序列上保守，表达量低，组织特异性强，通常与蛋白编码基因协同表达，共同参与生物过程)

2. 8种lncRNA的作用机制

[if !supportLists]4. [endif]lncRNA的结构特征

[if !supportLists]5. [endif]lncRNA发展史

[if !supportLists]6. [endif]编码潜能分析

[if !supportLists]7. [endif]lncRNA预测原理（预测其可能发挥的作用，预测具有假阳性）

[if !supportLists]8. [endif]iseeRNA支持人和小鼠的使用。