转录组分析(五)--王通老师讲解

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支持人和小鼠的使用。






 

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 206,723评论 6 481
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 88,485评论 2 382
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 152,998评论 0 344
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 55,323评论 1 279
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 64,355评论 5 374
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,079评论 1 285
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,389评论 3 400
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,019评论 0 259
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 43,519评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,971评论 2 325
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,100评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,738评论 4 324
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,293评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,289评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,517评论 1 262
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,547评论 2 354
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,834评论 2 345

推荐阅读更多精彩内容