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