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 上。
下载地址

推荐写一个 sh 脚本，批量下载，即列出全部的 连接。

而后使用 sra-tools 把 sra 格式序列转换成 fq 格式序列

脚本以下：

#!/bin/bash

cd /home/user/download/myfile/

wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871481/SRR1871481.sra

wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871482/SRR1871482.sra

wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871483/SRR1871483.sra

wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871484/SRR1871484.sra

wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871485/SRR1871485.sra

wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP%2FSRP055%2FSRP055992/SRR1871486/SRR1871486.sra

# use sra-tools to transform

export PATH=$PATH:/home/user/sratoolkit.2.5.2-ubuntu64/bin

fastq-dump *.sra

这样就把原始序列 fq 文件获得了。

为了后面分析方便，把相应的序列文件名改为相应的组

mv SRR1871481.fastq WT_Rep1.fastq
mv SRR1871482.fastq WT_Rep2.fastq
mv SRR1871483.fastq WT_Rep3.fastq

mv SRR1871484.fastq athb1_Rep1.fastq
mv SRR1871485.fastq athb1_Rep2.fastq
mv SRR1871486.fastq athb1_Rep3.fastq

Pre-Alignment QC

使用fastqc 软件来对原始序列fastq 文件进行质量检测

export PATH=/home/maque/FastQC/:$PATH
fastqc *.fastq

这样每一个 fastq 文件都能生成一个 html 报告文件，很详细

使用 tophat 和 bowtie 进行比对

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 数据分析的简单过程，不少细节没有提，并且还有不少其余步骤能够进行扩展，这些还须要再进一步深刻理解。