rna-seq技术进行转录组学分析的原理,转录组RNA-seq

　　前段时间跟着单饼干大神视频学习，在自己的小破笔记本上完成了RNAseq差异表达(tophat2 cufflink cuffdiff)的整个分析过程。虽然这个流程比较老，但是现在也常用HTseq DESeq2等其他流程进行分析，但是作为入门知识还是很容易理解的。本文将先描述过程，再描述安装Linux子系统(更新)、anconda和一些分析软件(更新)的过程，为入门做一个真正完整的操作流程。

　　更新了R中火山地图和热图的可视化部分。

　　我的电脑配置真的没什么。

　　计算机配置

　　但是我一天就完成了整个过程。

　　运行环境python2.7

　　原始数据如下：

　　正本单据

　　包括四个文件：

　　Bowtie2_hg19索引文件(这里已经提前使用Bowtie2建立了索引文件，可以直接使用)

　　Raw_data illumina双端测序原始文件(对照组2个，实验组2个，即8个测序文件)

　　RRNA RRNA索引序列文件(用于去除rRNA的影响)

　　分析处理

　　RNA-seq原始数据质量的评价

　　原始数据的过滤和不可信数据的清除(干净读取)

　　阅读回复基因组和转录组(比对)

　　计数(计数)

　　基因差异分析

　　数据的下游分析(这次先不做这个，下次再单独写)

　　让我们开始正式分析。

　　1.fastqc质量控制

　　Dirfastqc _ result.raw #(建立输出文件夹)

　　fastqc-q-T3-o./fastqc _ result.raw/* . FQ . gz #(使用fast QC质量控制软件对所有raw_data进行质量控制测试)

　　2.multiqc集成质量控制文档(因为获得了许多质量控制测试结果，所以需要集成它们)

　　Multiqc。#(这一步是集成fastqc_reault文件夹下的所有文件)

　　整合后文档

　　3.根据质控结果，使用fastx_tolltik去除不良序列。

　　zcat HeLa _ ctrl _ re P1 _ R1 . FQ . gz fastx _ trimmer-f 11-l 140-z-o out _ re P1 _ R1 . FQ . gz

　　zcat HeLa _ ctrl _ re p2 _ R1 . FQ . gz fastx _ trimmer-f 11-l 140-z-o out _ re p2 _ R1 . FQ . gz

　　zcat HeLa _ ctrl _ re p2 _ R2 . FQ . gz fastx _ trimmer-f 11-l 140-z-o out _ re p2 _ R2 . FQ . gz

　　zcat HeLa _ ctrl _ re P1 _ R2 . FQ . gz fastx _ trimmer-f 11-l 140-z-o out _ re P1 _ R2 . FQ . gz

　　zcat HeLa _ treat _ re P1 _ R1 . FQ . gz fastx _ trimmer-f 11-l 140-z-o out _ t _ re P1 _ R1 . FQ . gz

　　zcat HeLa _ treat _ re P1 _ R2 . FQ . gz fastx _ trimmer-f 11-l 140-z-o out _ t _ re P1 _ R2 . FQ . gz

　　zcat HeLa _ treat _ re p2 _ R1 . FQ . gz fastx _ trimmer-f 11-l 140-z-o out _ t _ re p2 _ R1 . FQ . gz

　　zcat HeLa _ treat _ re p2 _ R2 . FQ . gz fastx _ trimmer-f 11-l 140-z-o out _ t _ re p2 _ R2 . FQ . gz

　　因为当时不会写shell脚本，也不懂正则表达式，所以一个一个写，但是shell就是生产力。这一步也可以放在后台。当时我一个个跑了，看结果。

　　Zcat解压，文件名，fastx_trimmer -f x -l y保留自x-y的sequence -z压缩命令-o输出结果后台运行

　　Trimmer可以使所有序列的长度相同。

　　4.cutadaptor移除read两端的适配器。

　　nohup cutadapt-times 1-e 0.1-o 3-quality-cut off 6-m 50-A AGATCGGAAGAGC-A AGATCGGAAGAGC-o cut _ out _ c _ re P1 _ R1 . FQ . gz-p cut _ out _ c _ re P1 _ R2 . FQ . gz out _ c _ re P1 _ R1 . FQ . gz out _ c _ re P1 _ R2 . FQ . gz out _ c _ re P1 _ R2 . FQ . gz

　　nohup cut adapt-times 1-e 0.1-o 3-quality-cut off 6-m 50-A AGATCGGAAGAGC-A AGATCGGAAGAGC-o cut _ out _ c _ re p2 _ R1 . FQ . gz-p cut _ out _ c _ re p2 _ R2 . FQ . gz out _ c _ re p2 _ R1 . FQ . gz out _ c _ re p2 _ R2 . FQ . gz out _ c _ re p2 _ R2 . FQ . gz

　　nohup cut adapt-times 1-e 0.1-o 3-quality-cut off 6-m 50-A AGATCGGAAGAGC-A AGATCGGAAGAGC-o cut _ out _ t _ re P1 _ R1 . FQ . gz-p cut _ out _ t _ re P1 _ R2 . FQ . gz out _ t _ re P1 _ R1 . FQ . gz out _ t _ re P1 _ R2 . FQ . gz out _ t _ re P1 _ R2 . FQ . gz

　　nohup cut adapt-times 1-e 0.1-o 3-quality-cut off 6-m 50-A AGATCGGAAGAGC-A AGATCGGAAGAGC-o cut _ out _ t _ rep 2 _ R1 . FQ . gz-p cut _ out _ t _ rep 2 _ R2 . FQ . gz out _ t _ rep 2 _ R1 . FQ . gz out _ t _ rep 2 _ R2 . FQ . gz out _ t _ rep 2 _ R2 . FQ . gz

　　-乘以1一个序列只去适配器；一次；

　　-e 0.1匹配时可以有10%的错误率；

　　-o 3衔接子序列必须与测序序列重叠3个碱基以上。

　　-质量-截止常用6；

　　-m 50处理后低于50就扔掉序列，短序列的测序质量可能不太好；

　　-a和-A是Illumina常用的常用引物。之所以输入两个引物，是因为它们是一个双端测序的结果，需要分别去掉两个文件的内容。-a对应于读取1，而-A对应于读取2。

　　-o前一输出快速q1 -p前一输出快速q2

　　最后，编写日志文件。

　　//其中nohup:不挂机运行命令。

　　在后台

　　//2

　　1是传递给shell脚本的第一个参数；

　　5.使用bowtie2将读数与rRNA进行比较，并消除rRNA的影响。

　　nohup bow tie 2-x $ rRNA _ index-1 cut _ out _ c _ re P1 _ R1 . FQ . gz-2 cut _ out _ c _ re P1 _ R2 . FQ . gz-S Sam _ out _ c _ re P1 _ bow tie-p 2-un-conc-gz fastq _ unmap _ c _ re P1 _ R1

　　nohup bow tie 2-x $ rRNA _ index-1 cut _ out _ c _ rep 2 _ R1 . FQ . gz-2 cut _ out _ c _ rep 2 _ R2 . FQ . gz-S Sam _ out _ c _ rep 2 _ bow tie-p 2-un-conc-gz fastq _ unmap _ c _ rep 2 _ R1

　　nohup bow tie 2-x $ rRNA _ index-1 cut _ out _ t _ rep 2 _ R1 . FQ . gz-2 cut _ out _ t _ rep 2 _ R2 . FQ . gz-S Sam _ out _ t _ rep 2 _ bow tie-p 2-un-conc-gz fastq _ unmap _ t _ rep 2 _ R1

　　nohup bow tie 2-x $ rRNA _ index-1 cut _ out _ t _ re P1 _ R1 . FQ . gz-2 cut _ out _ t _ re P1 _ R2 . FQ . gz-S Sam _ out _ t _ re P1 _ bow tie-p 2-un-conc-gz fastq _ unmap _ t _ re P1 _ R1

　　(这里，我们需要先为rRNA建立一个索引库。如果别人建了库，可以直接下载使用。)

　　rRNA _ index=/mnt/c/Users/wt/Desktop/test _ data/RNA seq _ test _ date/rRNA/bow tie 2 _ rRNA _ human

　　一般用map与rRNA的比值来衡量数据库构建的质量。一般rRNA的比例不要超过10%。

　　-x对应于rRNA的索引序列；

　　-1，-2是reads1和reads2只是输出；

　　-S是比较结果的输出文件；

　　-p 2用2核计算(P是处理器！);

　　-un-conc-gz输出不匹配的结果；(你配不上的才是你需要的)

　　输出结果如下

　　输出结果

　　实际上，还应该输出日志文件来检查运行的进程。如果报告了错误，很容易检查和处理它。但是第一次做的时候输出是空白的，不知道哪里出了问题。这里的sam_out文件有问题。虽然没什么用，但是本来应该输出的是比较比率，是一个日志文件。我们稍后会看到这里发生了什么。

　　质量控制到此结束！

　　6.使用tophat2将过滤的rRNA读数与ref(参考)基因组进行比较。

　　(如果mRNA直接和人类DNA比对，可能会有问题，可能会跨越一个内含子。tophat2考虑到了这个问题，它根据注释文件把读操作分成短序列，并重新对齐；)

　　需要先建立一个数据库(别人已经在这里建立了)

　　hg19 _ index=/mnt/c/Users/wt/Desktop/test _ data/RNA seq _ test _ date/bow tie 2 _ hg19/hg19 _ only _ chronous

　　nohup top hat 2-p2-o top _ out 1 $ hg19 _ index fastq _ unmap _ c _ re P1 _ R1 . FQ . 1.1 . gz fastq _ unmap _ c _ re P1 _ R1 . FQ . 1.2 . gz

　　nohup top hat 2-p2-o top _ out 2 $ hg19 _ index fastq _ unmap _ c _ rep 2 _ R1 . FQ . 2.1 . gz fastq _ unmap _ c _ rep 2 _ R1 . FQ . 2.2 . gz

　　nohup top hat 2-P1-o top _ out 3 $ hg19 _ index fastq _ unmap _ t _ re P1 _ R1 . FQ . 1.1 . gz fastq _ unmap _ t _ re P1 _ R1 . FQ . 1.2 . gz

　　nohup top hat 2-P1-o top _ out 4 $ hg19 _ index fastq _ unmap _ t _ rep 2 _ R1 . FQ . 2.1 . gz fastq _ unmap _ t _ rep 2 _ R1 . FQ . 2.2 . gz

　　-o top_out1，结果输出到这个文件夹。

　　Hg19是人类基因组的第19个版本。hg19和hg38是常用的，hg38将取代hg38。hg19包含丰富的信息。

　　使用的序列是上一步没有比对的序列unmap(也就是我们需要的质控结果)。

　　输出结果如下

　　tophat2的输出结果。bam是最终的比较结果；txt是对比的总结；

　　3 bed文件，软件检测到缺失的删除，插入的插入和连接区域。info有一些reads没有直接映射到连续的DNA基因组上，所以需要切割一些reads，处理reads的进程文件保存在info中；

　　Unmapped.bam未映射，并且层之间没有比较。可能是未注释的基因组，测序问题等。

　　7.用袖扣评估表达水平(正常情况下这一步没用)

　　袖扣-o袖扣_out1 -p 4 -G./ref seq _ genes _ hg19 . GTF top _ out 1/accepted _ hits . bam

　　袖扣-o袖扣_out2 -p 4 -G./ref seq _ genes _ hg19 . GTF top _ out 2/accepted _ hits . bam

　　袖扣-o袖扣_out3 -p 4 -G./ref seq _ genes _ hg19 . GTF top _ out 3/accepted _ hits . bam

　　袖扣-o袖扣_out4 -p 4 -G./ref seq _ genes _ hg19 . GTF top _ out 4/accepted _ hits . bam

　　-G转录组的参考文献

　　袖扣的输出如下：

　　袖扣的输出结果

　　Gene.fpkm基因的fpkm计算结果(基因表达水平)

　　异构体的fpkm计算结果(转录表达)fpkm异构体(可以理解为基因的每个外显子)

　　skipped.gtf的跳过基因的转录信息

　　转录组的gtf，包含袖扣的装配结果同种型。

　　Fpkm是衡量基因表达的数值。一个基因有不同的内含子和外显子，不同的外显子之间可以形成不同的转录本。每个转录物可以翻译成不同的蛋白质，这些蛋白质是同种型(亚型)。对于不同的转录本，一个基因有一个表达，就是基因的fpkm和同工型的fpkm。

　　8.cuffdiff计算基因表达的差异。

　　ctrl _ bam=top _ out 1/accepted _ hits . bam，top_out2/accepted_hits.bam

　　treat _ bam=top _ out 3/accepted _ hits . bam，top_out4/accepted_hits.bam

　　label=hela_ctrl，hela_treat

　　cuffdiff-o diff _ out 1-p 7-labels $ label-min-reps-for-js-test 2./ref seq _ genes _ hg19 . GTFctrl _ bamtreat _ bam

　　- lables是文件的输入顺序，如上图label=hela.ctrl，hela _ treat

　　- min每次treat有几次重复，最上面有两个ctrl_bam，应该和treat_bam的次数一致，且=2(最小重复)

　　然后将其与参考抄本进行比较。

　　运行结果如下：

　　Cuffdiff结果

　　主要使用genes_exp.diff文件进行后续分析。

　　以上文件含义，查看cuffdiff输出文件的注释内容(有时间我再补充)

　　此时，作为rnaseq分析的一部分，可视化将单独进行。