完整转录组RNAseq分析流程（tophat2+cufflink+cuffdiff）

2024-01-21 09:35| 来源: 网络整理| 查看: 265

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

我的电脑配置，真的是战五渣。但还是在一天内跑完了整个流程

运行环境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：$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 分析结束一部分，可视化会另外做

【本文地址】

公司简介

联系我们