主要介紹如何分析RNA-seq 數據
參考文檔
paper: RNA-Seq: a revolutionary tool for transcriptomics
TopHat: discovering splice junctions with RNA-Seq
推薦：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 進行比對
在開始之前，需要下載擬南芥的基因組序列，注釋文件以及 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
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
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 數據分析的簡單過程，很多細節沒有提，而且還有很多其他步驟可以進行擴展，這些還需要再進一步深入理解。