前一段时间跟着孟浩巍大神的视频学习，在自己的小破笔记本上还是跑完了整个RNAseq差异表达的分析流程( tophat2 + cufflink + cuffdiff )虽然这个流程比较老了，现在做分析一般使用的都是 HTseq + DESeq2 等其他的流程，但是作为入门的知识还是比较容易理解，这篇文章先更一下流程，后面会再更一篇 安装 Linux 子系统(已更新)，安装 anconda 和一些分析软件(已更新)的流程，凑一个真正完整的入门生信的操作流程。

火山图、热图在 R 中可视化部分已更新

我的电脑配置，真的是战五渣。

电脑配置

但还是在一天内跑完了整个流程

运行环境python2.7

原始数据如下：

原始文件

包括四个文件：

bowtie2_hg19 index 文件(这里已经提前使用bowtie2建立好了index文件可以直接使用)

raw_data illumina 双端测序原始文件(对照组和实验组各两个，就是八条测序文件)

rRNA rRNA index 序列文件(用于去除 rRNA 的影响)

分析流程

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(也就是我们所需要的质控后的结果)

输出结果如下

tophat2输出结果

.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输出结果如下：

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(最小重复)

然后比对到参考转录本上

运行结果如下：

cuffdiff 结果

主要用到genes_exp.diff文件后续分析

以上文件含义查看 cuffdiff 输出文件的笔记内容(有时间我补充上来)

至此 rnaseq 分析结束一部分，可视化会另外做