E-mail: yzeu@qq.com分享 http://blog.sciencenet.cn/u/Jasion (Personal website:http://yz2.net);(School email: 2015212003@student.cup.edu.cn)

博文

RNA测序数据处理

已有 2944 次阅读 2018-4-12 21:16 |个人分类:所见所闻|系统分类:科研笔记

 

深度测序数据研究转录组 explore transcriptome using NGS

An Overview

 

转录组:细胞全体转录本的集合,也是特定时刻基因表达谱的快照。

包含编码蛋白的信使RNA  mRna ,非编码RNA :micrRNA lcnRNA

彼此协作共同调控细胞发育。

对转录组的研究包括定性和定量两方面。

EST expressed sequence tag表达序列标签  技术

随机选择的CDNAD进行克隆,单次测序获得此CDNA部分序列。60—500bp 1%错误率

全转录组水平的深度测序 RNA-seq

过程: 生物sample中的RNA反转录成reverse transcription +RCR CDNA,然后将这些CDNA打碎为较小片段之后上机测序。

1.      快速确定转录组中所有可变剪切体。2.通过对基因组上特定位点的reads计数,估计表达量的水平。对转录组同时进行定性、定量研究。

RNA-seq本质是对转录本的随机采样。The detection power is highly dependent on sequenceing depth. 随机抽样的情况下,认为map到特定转录本上reads的数目正比于它的表达量,也正比于转录本长度和测序深度。

从原始计数到表达水平:

RPKM方法 

C:贴到特定转录本上reads总数

N:这次实验总的reads数,也就是测序深度

L:特定转录本序列的长度

基因表达量算法类似。

在假定不同样品中RNA总体分布一致的情况下,可以用RPKM比较不同基因、不同样品之间的表达数据比较。

后续提出的TMM、DESeq效果更好。

测序数据按照是否分链进行分类,包含三类:

链特异性,DNA两条链都可以转录,形成不同转录本。常用iilumina rna-seq是不分链的,无法知道配对的reads哪个方向和转录本一致,哪个方向是反向互补的,standard iilumina fr-unstranded

对于分链数据存在两种情形,dutp标记的方法,第二个read方向和转录本一致。第一个read和转录本反向互补。Iilumina Fr-firststrand .

在另一种solid平台就刚好相反,fr-secondstrand,在分析数据之前一定弄清楚数据有没有分链,如何分链。

数据归一化参考文献:

A comprehensive  evaluation of normalization methods for iiumina high-throughput RNA sequencing data analysis. Briefings in bioinformatics.

 

测序数据长度短:36~80bp per read

数据量大、质量参差不齐、错误率高

 

DNA转录成Mrna过程中内含子被切掉,外显子在剪切位点连接到一起,这些跨过剪切位点的reads叫做junction reads 。如果无法把junctionreads正确从中断开就不能正确贴到基因组上。Junction reads是确定剪切位点的直接证据,对正确重构转录本至关重要

算法在mapping时要能够识别junction和intron的存在,目前有两种策略:

1.join exon已知转录本的所有exon构建可能存在的junctions组合。

2,split reads。无法直接map到基因组的reads切分成长度为k的种子,在更小的粒度上尝试mapping。最后将临近的reads组合起来得到全域的alignment。此方法虽然速度慢,但不依赖已知的exon注释,可用来发现新的exon或新的基因。

目前的工具组合了两种策略,平衡灵敏度和速度,例如tophat2

完成mapping只是数据分析的第一步。

还需要将这些reads组装成转录本,并对每个转录本估计相应的表达量。

包括junction reads在内的所有reads都正确map到基因组之后,转录本的组装问题变为一个有向图的遍历问题,找到一条或多条符合约束的最优路径。以及相对应的转录本序列。

常用工具cufflinks,针对RNA-seq的转录本组装和表达分析的软件。

原则上在转录本组装正确完成,且每一个exon表达量均归一化处理之后,转录本表达量可以直接从各个exon表达量推出来。

实际reads分配中cufflinks还要考虑reads长度分布,反复迭代考虑转录本组装和表达量估计。

 

Tophat、 cufflinks(包含cuffmerge、 cuffdiff) 是macOS和Linux下的工具

CummeRbund是R的工具包

 

Tophat将RNA-seq回帖到基因组上,原理:先用bowtie把reads贴到外显子上,形成堆叠的consensus。经典的剪切位点一般都有GT和AG这样的序列标志,tophat会去consensus中找这样的位点,并得到可能的组合。把没有被bowtie贴到基因组的reads,建立索引和可能的剪切位点比对。

Tophat中的命令行选项:

插入片段长度的选项  -r  /--mate-inner-dist <int>

Mrna打断后筛选的片段长度300bp 两端测序各75bp, 则插入长度为150bp


--mate-std-dev <int> 设置插入片段长度的标准差

如果筛选用的片段长度分布较为集中,这个值应当设置小一些,反之大一些。

 

-G/--GTF <GTF/GFF3 FILE> 提供已有的基因组注释文件,最好能提供此文件。Tophat会先把reads贴到转录组上,一些junction reads可以直接贴上去,剩下的往基因组上贴。最后合并结果。

 

--library-type 文库分链类型

Standard iilumina不分链 fr-unstranded 无法知道配对的reads哪个和转录本一致,哪个和转录本反向互补。

 

分链数据也有两种情况:

fr-firststrand :Dutp, NSR ,NNSR 第二个read和转录本方向一致,第一个与转录本反向互补。

fr-secondstrand : Ligation, Standard Solid ,与上述相反。

操作演示

-p是线程数      -G是注释文件    -o是输出文件夹     genome是基因组索引index

后面是两个reads的fastq文件

需要:

1. 基因组序列fasta文件 genome.fa

2. 建立bowtie2的索引index文件  bowtie2-build genome.fa genome

得到.bt2的文件就是需要的index文件

3.测序得到的.fastq文件

4.基因组注释文件 .gtf  .gff3 (本文件可选)给tophat提供剪切位点的信息

        1       3        4          5         7                     9         

染色体  exon   起始位置    结束位置    +在正链上;/-     第几个外显子

  

后面是基因、转录本的id和名字

 

.log是tophat运行的整个过程记录

 

align_summary.txt是reads回帖的统计信息,90%以上回帖率就非常好了,70%也可以接受。

 

Accepted_hits.bam 详细记录了reads回帖到基因组上的情况。是一个二进制文件,无法直接查看,需要samtools,  samtools view *hist.bam | less

Bam的每一行是read回帖到基因组情况的描述

1             2                3                      4

read名字   回帖情况描述   贴到参考序列的名字      read最左端在参考序列上位置

 

   5             6 :/75M/66M76N9M                          7

  回帖质量   贴上的碱基数目,可能贴到剪切位点     和它配对的read在参考序列的位置

 

    8

 片段长度

 

XS:A:+ 正链

还有和它配对的read是否贴上

运行结束后生成5个文件:

Accepted_hits.bam  是reads mapping信息

deletions.bed     insertions.bed   splice_junctions.bed插入、缺失、剪切位点

  align_summary.txt 对mapping信息统计的文件

 

Cufflinks是一套拼接转录本,计算表达量,计算差异表达的工具

通过所有read map结果拼接转录本,并估计表达量。

命令行选项:

-G/--GTF <reference_annotation.gtf/gff> 提供注释文件,不允许拼接新的转录本

 

-g/--GTF-guide <reference_annotation.gtf/gff> 提供注释文件,允许新的转录本

 

-u/--multi-read-correct  使用更加准确的方法(EM算法迭代)处理同时贴到多个位点的

reads,没有-u则将reads简单分配到各个位点。

 

 

--library-type  和tophat一样 the same as tophat 分链

 

演示:

输出:

transcripts.gtf是cufflinks拼接的转录本,包含转录本位置、RPKM信息

会生成三个文件:

基因表达文件  gene expression

转录本表达文件  transcript expression

组装好的转录本  assembled transcripts

 

使用cufflinks处理多个样本之后,cuffmerge将转录本数据.gtf整合为一个全面的转录本集合。

上图将6个转录本整合为一个转录本,还可以同时利用基因组注释文件获得更加全面可靠的结果。合并后转录本集合为计算每个基因、转录本的表达量提供了统一的基础。

 

命令行选项:

-g/--ref-gtf  参考注释gtf文件

 

-p/--num-threads 线程数

 

-s/--ref-sequence <seq_dir>/<seq_fastq> 指向基因组的DNA序列

 

演示

assemblies.txt 是cufflinks拼接的所有转录本文件路径  .gtf

 

 merge_asm文件夹储存了最后merge的结果,包含.gtf整合的转录本文件

 

cuffdiff

利用cufflinks获得了拼接的转录本时,可以计算不同样品中转录本的表达量

原理:测序深度和外显子长度一定时,reads数量与转录本表达水平成正比。

对reads计数可以计算转录本的表达量,同时cuffdiff可以计算不同条件下转录本表达水平的显著性差异。

 

命令参数:

-u/--multi-read-correct  默认关闭 ,对回帖到基因组中多个位置的reads进行估计后加权分配到各个基因组位置。功能用法与cufflinks中-u相似

 

-L/--labels  为每个样本标上名称

 

演示

 


 

merge.gtf  是cuffmerge产生的.gtf文件,此注释用于初始转录产物、可变剪切定量分析。

 

  .bam 是tophat产生的bam文件,如果一个样品有多个实验重复,需要提供.bam文件列表,文件之间逗号隔开。

根据需要设置参数,如 FDR(0.05)、 normalization方法(几何平均数)、min alignment count

输出文件在diff_out目录,包括按类别统计的表达水平结果,如:相同的转录起始位点、相同编码区的转录本表达水平...... 可以利用它们进行后续分析

Cuffdiff生成文件共有21个:

第一类是 FPKM, 保存基因、转录本、CDS区、TSS的表达量

第二类是count文件,保存每个FPKM对应的read cout

..................................................

差异表达检验文件,包含差异表达基因的相关注释信息以及统计检验相关参数,significant列是yes,表明此基因在两组中表达差异显著。

CDS FPKM差异表达分析需要基因的所有转录本都有p_id属性。 Isoforms of gene

 

数据分析后可视化 cummeRbund 用于分析cufflinks套件RNA-seq处理后的文件。

常用绘图命令:


密度分布图、散点图、火山图、柱形图

csBoxplot csHeatmap

R 进入R环境

 

Fold change 不同条件下基因表达水平比值

 

可对特定基因在不同条件下表达水平分析,特定基因不同转录本表达差异

 

 

可能的影响因素 possible effect factors

Reads quality

Sequencing depth

Mapping parameters

Normalization methods of cufflinks/cuffdiff

Biology duplicate  实验重复的次数,越多越可靠。




http://blog.sciencenet.cn/blog-734913-1108768.html

上一篇:笔记
下一篇:matlab 调用马夸特算法

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2019-11-18 09:06

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部