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

  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分析的一部分,可视化将单独进行。

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