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 數據分析的簡單過程，不少細節沒有提，並且還有不少其餘步驟能夠進行擴展，這些還須要再進一步深刻理解。