http://www.fungenomics.com/article/30 【专题】基因组学技术专题（二）—— 为什么说FPKM/RPKM是错的
下载数据
wget是linux下一个从网络上自动下载文件的常用自由工具。它支持HTTP，HTTPS和FTP协议，可以使用HTTP代理。一般的使用方法是: wget + 空格 + 参数 + 要下载文件的url路径，例如：
1wget http://www.linuxsense.org/xxxx/xxx.tar.gz
Wget常用参数
-b：后台下载，Wget默认的是把文件下载到当前目录。
-O：将文件下载到指定的目录中。
-P：保存文件之前先创建指定名称的目录。
-t：尝试连接次数，当Wget无法与服务器建立连接时，尝试连接多少次。
-c：断点续传，如果下载中断，那么连接恢复时会从上次断点开始下载。
-r：使用递归下载
格式转换
// use sra-tools to transform > fastq-dump *.sra
为了后面分析方便，把相应的序列文件名改成相应的组
mv SRR1871481.fastq WT_Rep1.fastq
质量控制
Pre-Alignment QC
使用fastqc 软件来对原始序列fastq 文件进行质量检测 export PATH=/home/maque/FastQC/:$PATH fastqc *.fastq 这样每个 fastq 文件都能生成一个 html 报告文件，很详细
序列比对
使用 tophat 和 bowtie 进行比对
待添加
tophat2 -p 8 --bowtie1 -G genes.gtf -o WT_Rep1_output ../genome WT_Rep1.fastq
其他5个 fastq 文件与上面一致
-p 8 使用8线程
--bowtie1 使用bowtie1 , 默认是bowtie2
-G 后面跟注释文件 gtf
-o 后跟输出文件夹
最后面跟 原始序列 fastq 文件
经完成比对，生成了一个 accepted_hits.bam 文件， 这个就是 samtools 生成的，后面主要也是利用这个文件继续分析。 建议提前利用 IGV 软件查看一下 该 bam 文件，可以知道mapping 的情况。
查看bam文件
如果想查看，需要先对 该bam文件进行 index ,比如： samtools index WT_Rep1_output/accepted_hits.bam Use Cufflinks to generate expression estimates from the SAM/BAM files
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 数据分析的简单过程，很多细节没有提，而且还有很多其他步骤可以进行扩展，这些还需要再进一步深入理解。
