前一段时间跟着孟浩巍大神的视频学习，在自己的小破笔记本上还是跑完了整个RNAseq差异表达的分析流程（ tophat2 + cufflink + cuffdiff ）虽然这个流程比较老了，现在做分析一般使用的都是 HTseq + DESeq2 等其他的流程，但是作为入门的知识还是比较容易理解，这篇文章先更一下流程，后面会再更一篇 安装 Linux 子系统（已更新），安装 anconda 和一些分析软件（已更新）的流程，凑一个真正完整的入门生信的操作流程。
火山图、热图在 R 中可视化部分已更新
我的电脑配置，真的是战五渣。
但还是在一天内跑完了整个流程
运行环境python2.7
原始数据如下：
包括四个文件：
- bowtie2_hg19 index 文件（这里已经提前使用bowtie2建立好了index文件可以直接使用）
- raw_data illumina 双端测序原始文件（对照组和实验组各两个，就是八条测序文件）
- rRNA rRNA index 序列文件（用于去除 rRNA 的影响）
-
RefSeq_genes_hg19.gtf 基因组注释文件
链接：https://pan.baidu.com/s/1Q4SwosKaG16-aYA74SaASQ
提取码：gss3
分析流程
- RNA-seq的原始数据（raw data）的质量评估
- raw data的过滤和清除不可信数据（clean reads）
- reads回帖基因组和转录组（alignment）
- 计数（count ）
- 基因差异分析（Gene DE）
- 数据的下游分析（这次先不做这个，下次会单独写）
下面开始正式分析
1、fastqc质控
mkdir fastqc_result.raw #(建立输出文件夹)
fastqc -q -t 3 -o ../fastqc_result.raw/ *.fq.gz & #（使用fastqc质控软件，对所有raw_data进行质控检测）
2、multiqc整合质控文件（因为得到很多的质控检测结果，需要整合）
multiqc . #(这一步就是对 fastqc_reault 文件夹下所有文件进行整合)
3、根据质控结果，使用fastx_tolltik去除不良序列
zcat hela_ctrl_rep1_R1.fq.gz | fastx_trimmer -f 11 -l 140 -z -o out_rep1_R1.fq.gz &
zcat hela_ctrl_rep2_R1.fq.gz | fastx_trimmer -f 11 -l 140 -z -o out_rep2_R1.fq.gz &
zcat hela_ctrl_rep2_R2.fq.gz | fastx_trimmer -f 11 -l 140 -z -o out_rep2_R2.fq.gz &
zcat hela_ctrl_rep1_R2.fq.gz | fastx_trimmer -f 11 -l 140 -z -o out_rep1_R2.fq.gz &
zcat hela_treat_rep1_R1.fq.gz | fastx_trimmer -f 11 -l 140 -z -o out_t_rep1_R1.fq.gz &
zcat hela_treat_rep1_R2.fq.gz | fastx_trimmer -f 11 -l 140 -z -o out_t_rep1_R2.fq.gz &
zcat hela_treat_rep2_R1.fq.gz | fastx_trimmer -f 11 -l 140 -z -o out_t_rep2_R1.fq.gz &
zcat hela_treat_rep2_R2.fq.gz | fastx_trimmer -f 11 -l 140 -z -o out_t_rep2_R2.fq.gz &
因为当时我还不会写 shell 脚本，正则表达式也不懂，就一个一个写了，但是 shell 才是生产力啊啊啊啊,这一步也可以放后台的，当时为了看结果就一个一个跑了
zcat解压缩，文件名，fastx_trimmer -f x -l y 保留从x-y的序列 -z压缩命令 -o输出结果 &后台运行
trimmer，可以使所有序列长度一致
4、cutadaptor去掉read两端的adaptor
nohup cutadapt --times 1 -e 0.1 -o 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o cut_out_c_rep1_R1.fq.gz -p cut_out_c_rep1_R2.fq.gz out_c_rep1_R1.fq.gz out_c_rep1_R2.fq.gz &
nohup cutadapt --times 1 -e 0.1 -o 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o cut_out_c_rep2_R1.fq.gz -p cut_out_c_rep2_R2.fq.gz out_c_rep2_R1.fq.gz out_c_rep2_R2.fq.gz &
nohup cutadapt --times 1 -e 0.1 -o 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o cut_out_t_rep1_R1.fq.gz -p cut_out_t_rep1_R2.fq.gz out_t_rep1_R1.fq.gz out_t_rep1_R2.fq.gz &
nohup cutadapt --times 1 -e 0.1 -o 3 --quality-cutoff 6 -m 50 -a AGATCGGAAGAGC -A AGATCGGAAGAGC -o cut_out_t_rep2_R1.fq.gz -p cut_out_t_rep2_R2.fq.gz out_t_rep2_R1.fq.gz out_t_rep2_R2.fq.gz &
--times 1一条序列只去一次Adaptor；
-e 0.1在匹配时可以有10%的错误率；
-o 3 adaptor序列必须和测序序列有3个碱基以上的overlap才可以；
--quality-cutoff常用6；
-m 50如果处理之后低于50的话就扔掉序列，短序列测序质量可能不是很好；
-a和-A是 Illumina 常用的通用引物，之所以输入两个，是因为是一个双端测序的结果，需要对两个文件内容进行分别去除，-a对应Reads1，-A对应reads2
-o 上一步输出fastq1 -p 上一步输出fastq2
> 最后是写入log文件
//其中nohup：不挂断地运行命令。
& 就是放后台
//2>1 是传递给shell脚本的第一个参数；
5、利用bowtie2将reads比对到rRNA上，除去rRNA的影响
nohup bowtie2 -x $rRNA_index -1 cut_out_c_rep1_R1.fq.gz -2 cut_out_c_rep1_R2.fq.gz -S sam_out_c_rep1_bowtie -p 2 --un-conc-gz fastq_unmap_c_rep1_R1 &
nohup bowtie2 -x $rRNA_index -1 cut_out_c_rep2_R1.fq.gz -2 cut_out_c_rep2_R2.fq.gz -S sam_out_c_rep2_bowtie -p 2 --un-conc-gz fastq_unmap_c_rep2_R1 &
nohup bowtie2 -x $rRNA_index -1 cut_out_t_rep2_R1.fq.gz -2 cut_out_t_rep2_R2.fq.gz -S sam_out_t_rep2_bowtie -p 2 --un-conc-gz fastq_unmap_t_rep2_R1 &
nohup bowtie2 -x $rRNA_index -1 cut_out_t_rep1_R1.fq.gz -2 cut_out_t_rep1_R2.fq.gz -S sam_out_t_rep1_bowtie -p 2 --un-conc-gz fastq_unmap_t_rep1_R1 &
（这里要先为rRNA进行index建库，如果有别人建好的库可以直接下载使用）
rRNA_index=/mnt/c/Users/wt/Desktop/test_data/rnaseq_test_date/rRNA/bowtie2_rRNA_human
一般通过map到rRNA中的比例来衡量建库的质量。一般的要求rRNA的比例不超过10%。
-x 对应rRNA的索引序列；
-1，-2 是刚输出的reads1和reads2；
-S 是比对结果的输出文件；
-p 2 使用2个核心去运算（p就是processor吧！）；
--un-conc-gz 输出比对不上的结果；(比对不上的才是需要的)
输出结果如下
其实还应当将log文件输出，用于查看运行过程，如果报错容易查看进行处理，但是我第一次做的时候输出的都是空白，也不知道哪里有问题，。这里sam_out文件有点问题，虽然没有用，本应当输出的是比对上的比例，是一个log文件，后面会再看看这里怎么回事。
到这里才算质控做完！
6、使用 tophat2 将过滤掉 rRNA 的 reads 比对到 ref（参考）基因组上
（如果mRNA直接比对到人的DNA上，可能会出现问题，有可能跨越了一个内含子，tophat2考虑了这个问题，它将reads根据注释文件分开成短序列，重新比对；）
需要先建库（这里别人建好了）
hg19_index=/mnt/c/Users/wt/Desktop/test_data/rnaseq_test_date/bowtie2_hg19/hg19_only_chromosome
nohup tophat2 -p 2 -o top_out1 $hg19_index fastq_unmap_c_rep1_R1.fq.1.1.gz fastq_unmap_c_rep1_R1.fq.1.2.gz &
nohup tophat2 -p 2 -o top_out2 $hg19_index fastq_unmap_c_rep2_R1.fq.2.1.gz fastq_unmap_c_rep2_R1.fq.2.2.gz &
nohup tophat2 -p 1 -o top_out3 $hg19_index fastq_unmap_t_rep1_R1.fq.1.1.gz fastq_unmap_t_rep1_R1.fq.1.2.gz &
nohup tophat2 -p 1 -o top_out4 $hg19_index fastq_unmap_t_rep2_R1.fq.2.1.gz fastq_unmap_t_rep2_R1.fq.2.2.gz &
-o top_out1， 结果输出到这个文件夹中
hg19 是人的第19个基因组版本，常用的包括hg19和hg38，hg38会取代hg38，hg19包含的信息比较丰富
所使用序列是上一步未比对上的序列unmap(也就是我们所需要的质控后的结果)
输出结果如下
.bam 是最终的比对结果;
.txt 是比对中的总结情况；
3个bed文件，软件检测出来的 deletions 缺失， insertions 插入， junctions 区域
.info 有一些reads没有直接 map 到连续的 DNA 基因组上，需要切一些reads,加工 reads 的过程文件保存在 info 里；
unmapped.bam 是没有 map 上去，一层层都没有比对到的，可能是基因组上未注释过的、测序问题等。
7、使用cuffliks对表达量进行评估(这一步在正常情况下没什么用)
cufflinks -o cufflink_out1 -p 4 -G ../RefSeq_genes_hg19.gtf top_out1/accepted_hits.bam
cufflinks -o cufflink_out2 -p 4 -G ../RefSeq_genes_hg19.gtf top_out2/accepted_hits.bam
cufflinks -o cufflink_out3 -p 4 -G ../RefSeq_genes_hg19.gtf top_out3/accepted_hits.bam
cufflinks -o cufflink_out4 -p 4 -G ../RefSeq_genes_hg19.gtf top_out4/accepted_hits.bam
-G 转录组的参考文件
cufflinks输出结果如下：
gene.fpkm gene的 fpkm 计算结果（基因表达量）
isoforms.fpkm isoforms (可以理解为 gene 的各个外显子)的 fpkm 计算结果（转录本表达量）
skipped.gtf 跳过的基因的转录本信息
transcripts.gtf 转录组的gtf,该文件包含Cufflinks的组装结果isoforms
fpkm 是衡量基因表达量的数值，一个基因有不同的内含子和外显子，不同的外显子之间可以形成不同的转录本，每一个转录本可以翻译成不同的蛋白，这些蛋白互相之间就是 isoforms（亚型），对于不同的转录本来说基因有一个表达量，这就是基因的 fpkm 和 isoform 的 fpkm 。
8、cuffdiff计算基因表达差异
ctrl_bam=top_out1/accepted_hits.bam,top_out2/accepted_hits.bam
treat_bam=top_out3/accepted_hits.bam,top_out4/accepted_hits.bam
label=hela_ctrl,hela_treat
cuffdiff -o diff_out1 -p 7 --labels $label --min-reps-for-js-test 2 ../RefSeq_genes_hg19.gtf $ctrl_bam $treat_bam
--lables 是文件的输入次序，如上 label=hela.ctrl,hela_treat；
--min 每个 treat 里有几个 repeat ，上边 ctrl_bam 是两个，要和 treat_bam 数量一致且>=2（最小重复）
然后比对到参考转录本上
运行结果如下：
主要用到genes_exp.diff文件后续分析
以上文件含义查看 cuffdiff 输出文件的笔记内容（有时间我补充上来）
至此 rnaseq 分析结束一部分，可视化会另外做
