RNA-seq 数据的简单分析

主要介绍如何分析RNA-seq 数据

参考文档

Wikipedia-RNA-seq

paper: RNA-Seq: a revolutionary tool for transcriptomics

TopHat

Cufflinks

CummeRbund

TopHat: discovering splice junctions with RNA-Seq

TopHat2 paper

Nature Protocol

推荐：griffithlab/rnaseq_tutorial
以下的文档就按上面的这个教程来组织

软件安装

需要安装的软件：
sra-tools，samtools, bam-readcount, bowtie, bowtie2, tophat, star, cufflinks, htseq-count, R, cummeRbund, fastqc, picard-tools, and samstat

其中已经安装的：

sra-tools
samtools
bowtie
bowtie2
tophat
cufflinks
R
fastqc
samstat

案例

以下以一个例子来说明如何进行一般的 rna-seq的分析
GEO number : GSE66666
GSE66666

从中了解实验是如何设计的，想解决什么问题，多少样本，该研究所发表的文章等相关信息。

下载原始序列

原始序列一般以 SRA 的格式保存在 NCBI 上。
ref

在开始之前，需要下载拟南芥的基因组序列，注释文件以及 INDEX文件
iGenomes
选择 Ensembl tigr10 版本, 解压

cd /media/文档/rna-seq-arabi/ 
#原始序列与 tigr10 文件夹都放在该文件夹下

export PATH=/home/maque/samtools-1.2/bin:$PATH
export PATH=/home/maque/tophat-2.1.0.Linux_x86_64/:$PATH
export PATH=$PATH:/home/maque/bowtie-1.1.2

tophat2 -p 8 --bowtie1 -G Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Annotation/Genes/genes.gtf -o WT_Rep1_output  Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Sequence/BowtieIndex/genome WT_Rep1.fastq

# 其他5个 fastq 文件与上面一致
# -p 8    使用8线程
# --bowtie1   使用bowtie1 , 默认是bowtie2
# -G    后面跟注释文件 gtf
# -o    后跟输出文件夹
# 最后面跟 原始序列 fastq 文件

这些过程完成后，说明已经完成比对，在每个新建的文件夹里面，应该有一些信息，最主要的是生成了一个 accepted_hits.bam 文件， 这个就是 samtools 生成的，后面主要也是利用这个文件继续分析。

建议提前利用 IGV 软件查看一下 该 bam 文件，可以知道mapping 的情况。
如果想查看，需要先对 该bam文件进行 index ,比如：

samtools index  WT_Rep1_output/accepted_hits.bam

Use Cufflinks to generate expression estimates from the SAM/BAM files

ref

export PATH=/home/maque/cufflinks-2.2.1.Linux_x86_64/:$PATH
cufflinks -p 8 -o WT_Rep1_cuffout WT_Rep1_output/accepted_hits.bam
# 其他5个与之类似
# -p 8   使用8线程
# -o WT_Rep1_cuffout    输出目录
#  最后面跟相应的 bam 文件

该过程完成后，会生成相应的文件夹，里面有相应的文件，后面会使用 transcripts.gtf 文件

Differential Expression

ref

ls -1 *cuffout/transcripts.gtf > assembly_GTF_list.txt
cuffmerge -p 8 -o merged -g Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Annotation/Genes/genes.gtf -s Arabidopsis_thaliana_Ensembl_TAIR10/Arabidopsis_thaliana/Ensembl/TAIR10/Sequence/WholeGenomeFasta/genome.fa assembly_GTF_list.txt

# -p 8   使用8线程
# -o merged  后跟目录
# -g  后跟参考基因的gtf文件
# -s  后跟基因组序列文件
#  最后跟上一步新建的 assembly_GTF_list.txt

接下来使用 cuffdiff

cuffdiff -o rna_de/diff_out -p 8 -L WT,athb merged/merged.gtf WT_Rep1_output/accepted_hits.bam,WT_Rep2_output/accepted_hits.bam,WT_Rep3_output/accepted_hits.bam athb_Rep1_output/accepted_hits.bam,athb_Rep2_output/accepted_hits.bam,athb_Rep3_output/accepted_hits.bam

# -o  后跟输出文件目录 
# -p 8 使用8线程
# -L WT,athb    '-L' tells cuffdiff the labels to use for samples
# 后跟 上一步由 cuffmerge 生成的 merged.gtf 文件
# 最后跟6个bam 文件， 组内由逗号隔开，组间由空格隔开。

该过程会新建一个diff_out 文件夹，里面包含了很多信息
这些信息可以使用 R 包 cummeRbund 很方便的进行分析

使用cummeRbund

推荐流程

可以按照推荐流程文档中的步骤去分析

下面主要说一些注意点：

安装

source("http://bioconductor.org/biocLite.R")
biocLite("cummeRbund")

读入数据

首先先 cd 到上一个步骤生成的 diff_out 文件夹

library(cummeRbund)
cuff <- readCufflinks()

这样即完成读入数据。

各种分析图

可以按照推荐流程中的去分析，也可以参考 Nature Protocol

寻找差异表达基因

大部分进行RNA-Seq 实验的目的主要还是寻找实验组与对照组之间的差异表达基因。

一种方法是:

mySigGeneIds<-getSig(cuff,alpha=0.05,level='genes')
mySigGeneIds
length(mySigGeneIds)

mySigGenes<-getGenes(cuff,mySigGeneIds)
mySigGenes

diffData(mySigGenes)
featureNames(mySigGenes)

另一种方法是：

gene_diff_data <- diffData(genes(cuff))
sig_gene_data <- subset(gene_diff_data, (significant == 'yes'))
sig_gene_data

这些方法列出的 gene_id 是像 XLOC_000268 这样的格式， 它所对应的通用的gene_id 比如AT1G06080 , 这种一一对应关系文件可以在合并的 merged.gtf 文件中寻找
而AT1G06080 这种gene_id 所对应的基因类型，基因名称等信息可以在 tair10 文件夹中的 gene.gtf 文件中找到
AT1G06080 这种gene_id 所对应的基因名称也可以在 同一文件夹中的 refFlat.txt 文件中找到。

也可以先把上一步输出的 gene_id 放到一个 list.txt 中，注意不要有空行，最好使用 vim , 然后利用 grep 即可：

    grep -f  list.txt  merged.gtf  | less - S

以上就是rna-seq 数据分析的简单过程，很多细节没有提，而且还有很多其他步骤可以进行扩展，这些还需要再进一步深入理解。