RNAseq定量分析方案

RNAseq定量分析方案
RNAseq定量分析方案 1
一、实验目的： 2
二、实验大致流程 2
三、实验前的准备活动 2
3.1准备数据 2
3.2确定涉及软件是否安装完毕，及其输入输出文件。

3
四、实验过程 4
4.1.利用bowtie2-bulid命令根据提供的参考基因组序列建立对应的索引文件。

5
4.2.利用tophat命令将分别将代比较的reads maping到参考基因组序列上 6
4.3利用cufflinks软件分别将待测样品的转录组reads拼接起来，并同时计算每个样品各
个基因的rpkm值 7
4.4.利用coffmerge和cuffdiff软件计算每个样品各个基因的fpkm值。

8
五．利用R软件查看结果文件。

10
一、实验目的：
利用已有的水稻基因组数据对来自两棵不同的水稻进行转录组水平差异的研究。

二、实验大致流程
1. 利用bowtie2-bulid命令根据提供的参考基因组序列建立对应的索
引文件
2. 利用tophat命令将分别将代比较的reads maping到参考基因组序
列上
3. 利用cufflinks软件分别将待测样品的转录组reads拼接起来，计算
每个样品各个基因的fpkm值
4. 利用coffmerge和cuffdiff软件计算每个样品各个基因的fpkm值。

5. 利用R软件查看比较结果。

三、实验前的准备活动
3.1准备数据
像大多数生物实验一样，做生物信息学实验之前也是需要事先准备好“药品”和“仪器”，不然到了关键时刻也是一样会手忙脚乱的。

现在我们先来谈一谈生物信息实验需要准备的“药品”—数据。

对于RNAseq实验而言，这里至少需要准备一下几个文件：
1.参考基因组序列文件如 refrence.fa
参考基因组数据：
2.参考基因组注释文件如 refrence.gtf
a参考基因组序列设为refrence.fa b参考基因注释设为refrence.gtf
F1_sample_R1.fastq
待比较样品F1：
(为RNA测序数据) F1_sample_R1.fastq
P1_sample_R1.fastq
待比较样品P2：
(为RNA测序数据) P1_sample_L007_R2.fastq
3.2确定涉及软件是否安装完毕，及其输入输出文件。

通过上面的大致流程，我们知道，我们至少需要用到以下几个软件，查看你的系统是否安装了这些软件只需要直接输入这些命令就可以了：
1.bowtie2-build
/local/doc/bowtie2.html
Usage:
bowtie2-build [options]* <reference_in> <bt2_index_base>
reference_in comma-separated list of files with ref sequences bt2_index_base write .bt2 data to files with this dir/basename 2.tophat
/tutorial.shtml#toph
TopHat maps short sequences from spliced transcripts to whole genomes. Usage:
tophat [options] <bowtie_index> <reads1[,reads2,...]> [reads1[,reads2,...]] \
[quals1,[quals2,...]] [quals1[,quals2,...]]
3.cufflinks
/manual.html
Usage:
cufflinks [options] <hits.sam>
4.cuffmerge
/manual.html
Usage:
cuffmerge [Options] <assembly_GTF_list.txt>
5.cuffdiff
/manual.html
Usage:
cuffdiff [options] <transcripts.gtf> <sample1_hits.sam> <sample2_hits.sam> [... sampleN_hits.sam]
Supply replicate SAMs as comma separated lists for each condition: sample1_rep1.sam,sample1_rep2.sam,...sample1_repM.sam 由于这些命令的参数都不少，所以不能一一在这里讲解，想要了解这些软件的参数作用最好的方法就是直接读这些软件的manual。

四、实验过程
本实验数据存储在/share/Public/BioinfoTraim下面的Train_RNAseq目录里。

如果你想做练习的话，只要在/share/Public/BioinfoTraim将数据重新拷贝一份，然后再拷贝的数据下进行操作就可以了。

方法如下所示：
现在正式进入RNAseq的实验：
1.进入到Train_RNAseq这个目录里面(你则进入到自己拷贝的目录里面)，查看一下该目录下都有些什么文件。

我们可以看到，在该目录下，有我们需要的所有数据文件，说明数据已经准备就绪。

4.1.利用bowtie2-bulid命令根据提供的参考基因组序列建立对应的索引文件。

我们知道参考基因组的序列文件时Osativa.fa文件，所以这就是bowite2-bulid命令的输入文件，而bowtie2-build的输出文件时什么呢？根据bowtie2-build的Usage:
bowtie2-build [options]* <reference_in> <bt2_index_base>
reference_in comma-separated list of files with ref sequences
bt2_index_base write .bt2 data to files with this dir/basename
可知，bowtie2-build的输出文件名可以由我们制定的，实际上是输出文件不止一个，我们制定的只是这些输出文件的文件名的前缀而已。

照此说来，我们这里可以用这句命令来实现我们的目标：
# bowtie2-build Osativa.fa indexname
但是你会看到输出的软件，把我的实验环境搞的软七八糟的，所以我不想这样做，所以我决定先建立一个文件夹，来储存这些输出文件。

所以我该用了下面的语句：
#mkdir refrence_index
#bowtie2-build Osativa.fa refrence_index/indexname
查看一下refrence_index文件下面都有些什么文件。

我们会在refrence_index这个目录下面看到很多个文件，主要为两种类型，一种是bt2，还有一种是rev.x.bt2。

但这些文件都有一个共同特点，文件名的前缀都是indexname，这时你应该能理解bowtie2-bulid这句命令的作用了吧。

4.2.利用tophat命令将分别将代比较的reads maping到参考基因组序列上
建立完参考基因组序列的索引文件之后，我们就该讲样品的数据reads拿去maping了。

也许你会问为什么要建立索引？实际上这是为了加快mapping的速度。

tophat最简单的格式是:
#tophat --max-intron-length 20000 索引文件样品R1文件样品R2文件
结合本实验环境我们可以写这样的命令:
#tophat --max-intron-length 20000 refrence_index/indexname
F1_sample_R1.fastq\ F1_sample_R2.fastq
(备注命令中的“\”是用来实现连接下一行的，第一二行本来应该为一个句子)
还是那句话，我很讨厌输出的一大堆文件弄的我的实验环境乱七八糟的，所以我还是决定为这些输出文件找个好归宿，所以我会用下面这几命令：
#tophat -o F1 --max-intron-length 20000 refrence_index/indexname\
F1_sample_R1.fastq F1_sample_R2.fastq
（备注：F1一般就令其为样品名）
考虑到我们还有参考基因组的注释文件，所以我们可以一起把这个文件加进去。

$tophat -o F1 -p 8 -G Osativa.gtf --max-intron-length 20000\
refrence_index/indexname F1_sample_R1.fastq F1_sample_R2.fastq
可见这句命令使用到了三类文件：
参考基因组序列的索引文件
参考基因组的注释文件
样品的数据文件，包括R1和R2两个文件。

这里用到了tophat的三个参数：
-o 指定输出文件的存放目录。

-G 指定参考基因组的注释文件（可有也可没有，有的话，最后
对于分析结果有帮助）。

-p 指定运行该命令使用多少个cpu ,没有指定情况下默认为一
个。

-p 8 则指定使用8个。

--max-intron-length 最大的intron长度。

Tophat会忽略长度大于该
值的donor/acceptor pairs，除非有long read支持。

根据实际情况
来制定。

最后得到的数据存放在了F1里面，有不少类型的文件，但最主要的是accepted_hits.bam文件。

对F1样品进行maping操作之后，同样需要对P1样本进行同样的操作，输出文件制定为P1。

$tophat -o P1 -p 8 -G Osativa.gtf --max-intron-length 20000\
refrence_index/indexname P1_sample_R1.fastq P1_sample_R2.fastq
4.3利用cufflinks软件分别将待测样品的转录组reads拼接起来，并同时计算每个样品各个基因的rpkm值
由cufflinks的命令格式Usage: cufflinks [options] <hits.sam> 可知，该命令只有一类输入文件，虽然文件类型要求是hits.sam，实际上数据hits.bam 文件也是可以的，这两个文件的关系，可通过samtools命令来实现转换。

(/s/blog_6c9eaa150100xi81.html）
对于cufflinks不想讲太多的东西，因为这个命令输入文件即我们上一步由tophat产生的F1或P1里的accepted_hits.bam。

所以对于样品F1的操作命令可以写成：
$cufflinks -p 8 -o F1_assemble_transcripts F1/accepted_hits.bam
cufflinks参说明：
-o F1_assemble_transcripts 指定输出文件的存放目录为
F1_assemble_transcripts 。

-p 8 则指定使用8个cpu跑这个命令。

你可以看看F1_assemble_transcripts下都有哪些文件，顺便看看最重要的这个genes.fpkm_tracking文件里面的内容。

自己把下面这个图片拉大，查看文件内容。

对样品P2进行同样的操作
$cd /share/Public/BioinfoTraim/Train_RNAseq
$cufflinks -p 20 -o P1_assemble_transcripts P1/accepted_hits.bam
4.4.利用coffmerge和cuffdiff软件计算每个样品各个基因的fpkm值。

这一个步骤中有一些地方会让一些初学者感到奇怪，比如
1. cuffmerge
/manual.html
Usage:
cuffmerge [Options] <assembly_GTF_list.txt>
这个命令中的assembly_GTF_list.txt 文件到底指的是什么?实际上，命令帮助里是有说到的“cuffmerge takes two or more Cufflinks GTF files and merges them into a single unified transcript catalog. Optionally, you can provide the script with a reference GTF, and the script will use it to attach gene names and other metadata to the merged catalog”这句话的意思是cuffmerge的作用是将两个或者更多cufflinks产生的trainscripts.gtf文件合成一个统一的转录组目录。

既然涉及到两个或者更多的将两个或者更多cufflinks产生的trainscripts.gtf文件，那么怎么告诉cuffmerge这个软件，我这些文件都在哪里呢？我们可以将这些文件所在的地方都存放在一个文本文件里，如惯例使用assembly_GTF_list.txt这个文件来存储这些信
息。

我们可以用下面这个语句，来填写assembly_GTF_list.txt这个文件。

$ find -name transcripts.gtf > assemblly_GTF_list.txt
2. assemblly_GTF_list.txt文件准备好之后，我们就开始运行cuffmerge
命令吧。

$cuffmerge -g Osativa.gtf -s Osativa.fa -p 8 assemblly_GTF_list.txt Cuffmerge参数说明
-o 指定输出文件存放的目录，如果你没指定，它会默认存到
merged_asm这个文件里。

-g 可选项，主要是指定参考基因组的注释文件。

-s 指定参考基因组的序列文件。

现在让我们来看一下，运行完之后的变化吧，多产生了一个merged_asm目录，这个文件里有logs目录和一个merged.gtf文件。

这时候让我给实验来上最后沉重的一击吧,在
/share/Public/BioinfoTraim/Train_RNAseq下运行：
$ cuffdiff -o F1-vs-P1 -b Osativa.fa -p 20 -L F1,P1 -u
merged_asm/merged.gtf\
F1/accepted_hits.bam P1/accepted_hits.bam
Usage:
cuffdiff [options] <transcripts.gtf> <sample1_hits.sam> <sample2_hits.sam> [...
sampleN_hits.sam]
由该命令即可知，cuffdiff需要一下几个文件
cuffmerge阶段产生的merged.gtf文件
tophat产生的 F1和P1的xx_hits.bam文件
cuffdiff参数说明
-o 指定输出文件存放的目录
-b 使用偏定校正，后面接参考基因组的序列，实际作用还不是
很清楚。

-p 指定运行该命令使用多少个cpu
-L 条件标签逗号分隔的列表
-u use 'rescue method' for multi-reads
五．利用R软件查看结果文件。

1.下载安装R软件
/share/link?shareid=2674820005&uk=1010159686
2.启动R软件，安装需要的软件包CummeRbund package.通过输入以下命令就可以了
> source('/biocLite.R')
> biocLite('cummeRbund')
Use WinSCP to copy folder F1-vs-P1 (/export/xxxx/rnaseq/cufflinks/F1-vs-P1) to your computer. Download it to C:\ so that you can remember the path.
> library(cummeRbund)
> cuff_data = readCufflinks('C:/F1-vs-P1')
> csDensity(genes(cuff_data))
> csScatter(genes(cuff_data), 'F1', 'P1')
> csVolcano(genes(cuff_data), 'F1', 'P1')
> gene_diff = diffData(genes(cuff_data))
> sig_gene_diff = subset(gene_diff, (significant == 'yes'))
> nrow(sig_gene_diff)
[1] 2
> sig_gene_diff
gene_id sample_1 sample_2 status value_1 value_2 log2_fold_change
43 XLOC_000043 F1 P1 OK 689.648 2386.220 1.79079
151 XLOC_000151 F1 P1 OK 0.000 734.691 Inf
test_stat p_value q_value significant
43 0.823764 5e-04 0.0265 yes
151 NA 5e-05 0.0053 yes
> mygene = getGene(cuff_data, 'LOC_Os01g01870')
> expressionBarplot(mygene)。