RNA-seq

1.什么是RNA-seq？

RNA-seq即RNA测序技术。又称转录组测序技术，是一种通过观察基因表达来分析整个基因组的技术，主要测序对象包括信使RNA（mRNA）、微RNA（smallRNA）和非编码RNA（noncoding RNA），用高通量测序技术进行测序分析，反映出它们的表达水平.

2.什么是高通量转录组测序：

高通量测序技术（High-throughput sequencing)是什么?

顾名思义，就是通量很高（对比sanger测序）的测序，又称“下一代”测序技术（Next-generation sequencing technology），一次对几十万到几百万条DNA分子进行序列测定，同时高通量测序使得对一个物种的转录组和基因组进行细致全貌的分析成为可能，所以又被称为深度测序（deep sequencing）。

转录组(transcriptome)是什么?

广义上指某一生理条件下,细胞内所有转录产物的集合,包括信使RNA,核糖体RNA,转运RNA以及非编码RNA;狭义上指所有mRNA的集合.(你的样品经过了某种处理,然后拿去提了总RNA,这个总RNA就是一个转录组.)

RNA-seq的应用

RNA-seq可以帮助我们了解：各种比较条件下，所有基因的表达情况的差异。它可以检测的差异有：正常组织和肿瘤组织之间的差异；也可以检测药物治疗前后，基因表达的差异；还可以检测在发育过程中，不同的发育阶段，不同的组织之间的基因表达差异等。在所有检测的差异类型中，最常见的就是检测所有mRNA的表达量的差异。同时，还可以检测RNA结构上的差异，例如：mRNA的剪接方式的差异，即“可变剪接”，还可以检测“融合基因”，同时还可以检测基因单点突变导致的SNP（Single Nucleotide Polymorphisom）.

简单理解RNA-seq（主要分3步）（参考b站statQuest）

背景1

假设有一群正常的神经细胞（蓝色）和一群变异的神经细胞（红色）

dcba97d2847d1e83e92d48e01daf660.png

那么，为什么它们会出现差别呢？是什么遗传机制导致了这个事情呢？ 因此，我们需要看一看它们的基因表达差异。

背景2

我们知道，每个细胞都是由一堆染色体组成，每个染色体由一堆基因组成，当然并不是所有的基因都是活跃的，只有一部分基因是可以表达，而表达的中间过程就要经历mRNA转录本，通过高通量测序，我们就能得知：哪些基因是活跃可以表达的，并且产生了多少转录本（也就是衡量基因表达量的指标） 74e365bfd928572cd141e9f2cd3c396.png

背景3

将正常的细胞测一遍，再将变异的细胞测一遍，得到它们的表达量，我们后来就是比较它们的表达量差异 22823d0b50a0b078b94f71303b382aa.png

可以看出，基因1在两组样本中差异不大或者没有差异；基因2在正常组织中表达，而在变异组中表达量很高，二者差别甚大；基因3有差别但比较小。

step1.准备工作

1.构建序列文库

注：我们建的文库用用的是cDNA而不是rna. 除了单链RNA不稳定外，还有一小部分原因是DNA的建库流程已经确定了，只要把RNA变成DNA后面的流程完全一样。

bd939af1e552072fb4ebfe9d1b6ddf8.png 分离RNA→将RNA打断成小片段→将小RNA片段反转录成DNA→加接头→PCR扩增（只有加上接头的测序片段才能被扩增）→质量检查QC（看下文库的浓度和片段长度）

对文库进行测序

c2cc4c8cd4e3fcfcdc370591fe965a1.png

一块测序板上能包含多于400,000,000个片段，垂直于测序板排列。 测序仪有四种颜色的荧光探针A、T、C、G，与测序片段上碱基互补，结合上就“放烟花”表示庆祝🎉（就是闪一下自己带的荧光，比如A带红光，G带蓝光，C绿光，T橙光）。当然，这一切都逃不过测序仪自带的高精度照相机的法眼【测序仪为什么贵？就是在于它的高精度照相机，想想要分辨这么微小的亮光，密密麻麻】许许多多的测序片段中同一排的碱基测完了，就把原来荧光的那个碱基冲掉了，再放下一个荧光碱基进来结合、放光，测序就是这样，结果就得到了raw data，就是fastq数据。

fastq测序文件

这里先附一个下载fastq文件的方法。例如若想下载SRR1552445，直接在服务器运行下面命令即可：

In [ ]:
~/software/sratoolkit.2.10.7-centos_linux64/bin/fastq-dump SRR1552445
In [ ]:
#或者挂机下载，这样就不用担心服务器中断再重新下载的问题了，命令为：
nohup ~/software/sratoolkit.2.10.7-centos_linux64/bin/fastq-dump SRR1552445 &
top -c #查看下载是否完成

在illumina的测序文件中，采用双端测序（paired-end），一个样本得到的是seq_1.fastq.gz和seq_2.fastq.gz两个文件，每个文件存放一段测序文件。首先，我们拿到的原始序列文件就是fastq，那么如何去理解fastq文件呢？ 84545381ca0af212279ff557f42bbb2.png

在illumina公司测得的序列文件经过处理以fastq文件协议储存为*.fastq格式文件。上面是一个fastq文件的格式第一行就是测序的坐标信息，告诉你这条reads的名字是什么。（reads（读长）是测序仪单次测序所得到的碱基序列，也就是一连串的ATCGGGTA之类的，不同的测序仪器，reads长度不一样。对整个基因组进行测序，就会产生上千万个reads）；

第二行就是我们测到的序列

第三行就是一个加号（好像没什么用）

第四行是质量信息，对应着上面的各个碱基

如果要仔细了解，请打开链接https://zhuanlan.zhihu.com/p/20714540

本次RNA-seq所用数据来自网页：

Mouse mammary data (fastq files): https://figshare.com/s/f5d63d8c265a05618137 You should download these files and place them in your /data directory.

原文章链接也附上 https://bioinformatics-core-shared-training.github.io/RNAseq-R/align-and-count.nb.html

如何对测序质量这些东西进行直观化？

用fastqc软件（命令为：fastqc 文件名），可以对fastq文件进行质检，质检的结果具体怎么看，请看下面链接https://zhuanlan.zhihu.com/p/20731723

若有多个fastq文件，推荐使用multiqc软件（命令为：multiqc .），multiqc可以对几个fastq报告文件进行总结并汇总到一个报告文件中，以便更直观的展示。

在明白了自己的测序数据质量之后，我们就要将数据中低质量的部分全部剔除掉，这里使用的软件是trim—galore

2.注释文件和参考基因组的获取（以小鼠为例）

下载参考基因组

这里已经下载了参考基因组，命名为chr1.fa.

下载注释文件

然而参考基因组是一部无字天书，要想解读书中的内容，需要额外的注释信息协助。因此我们要下载注释文件。本文的注释文件下载网页：

ftp://ftp.ensembl.org/pub/release-100/gtf/mus_musculus/

4d7aa499cbe0d26020ebaa155fbba30.png

准备工作完成之后，我们就得到了一份高质量的原始数据（raw data），从而正式进分析工作

step2.上游分析

比对到参考基因组

先将大的基因组序列打断成许多小片段，然后为了方便接下来寻找这些片段，需要对它们进行构建索引index（目的是标注每个小片段的位置）。再将测序的reads和基因组一样，也是打断成小片段，然后把它的小片段比对到基因组的小片段上.

回帖

首先，我们的fastq文件储存的数据是一个零散的状态，那么要怎么样把它恢复成打断前的状态？

这里我们就需要一个模板，按照模板，我们把序列排序，大概就长这个样子： 24c78a2167c81f5bd8b93af7151664b.png 这里的ref就是模板，即参考基因组，而我们的fastq文件本质就是一条一条的小序列，在模板的指引下，我们得到了他们原本在基因组上应该在位置，这一步就是回帖的含义。

回帖完之后，我们的回帖信息会被输入到一个文本文件：SAM文件（二进制位bam文件）总的来看，所谓的sam/bam文件就是记录回帖的序列是什么，回帖上多少，回帖的质量行不行，回帖到什么位置。而后就是用featureCounts结合注释文件gff/gtf，将转录本构建出来得到raw count.

下面演示如何在R里进行读取比对和计数

In [1]:
#读取数据文件夹
setwd("D:/3219673")
In [2]:
#首先将Rsubread包加载到R中
library(Rsubread)

Error in library(Rsubread): there is no package called 'Rsubread'
Traceback:

1. library(Rsubread)
In [3]:
#读取fastq.gz文件
#用list.files命令在数据目录中搜索所有的.fastq.gz文件，pattern参数采用正则表达式，$表示字符串的结尾，我们仅获取以"fastq.gz"结尾的文件.
fast.files<- list.files(path="./",pattern="fastq.gz$",full.names=TRUE)
fast.files
  1. './SRR1552444.fastq.gz'
  2. './SRR1552445.fastq.gz'
  3. './SRR1552446.fastq.gz'
  4. './SRR1552447.fastq.gz'
  5. './SRR1552448.fastq.gz'
  6. './SRR1552449.fastq.gz'
  7. './SRR1552450.fastq.gz'
  8. './SRR1552451.fastq.gz'
  9. './SRR1552452.fastq.gz'
  10. './SRR1552453.fastq.gz'
  11. './SRR1552454.fastq.gz'
  12. './SRR1552455.fastq.gz'
In [ ]:
#建立索引
buildindex(basename="chr1_mm10",reference="chr1.fa")

Mapping

In [ ]:
#比对，这里以一个bam文件为例
?align
A<-file.path("./","SRR1552444.fastq.gz")
align(index="chr1_mm10",readfile=A,output_file="SRR1552444.bam")
In [ ]:
#用promapped函数看mapped的情况
bam.files<-list.files(path="./",pattern=".bam$",full.names=TRUE)
bam.files#这里应该有12个bam文件
props<- promapped(files=bam.files)
props

Quality control

In [ ]:
#用Rsubread中的qualityScores函数查看每个碱基的质量得分
#首先提取文件SRR1552450.fastq.gz的100reads的质量得分
qs<- qualityScores(filename="./SRR1552450.fastq.gz",nreads=100)
dim(qs)
head(qs)
In [ ]:
#查看100reads质量得分的总体分布，画一个箱线图
boxplot(qs)

Counting

In [ ]:
#使用featureCounts函数
gtf<-list.files("./",pattern=".gtf$",full.names=T)
fc <-featureCounts(bam.files,annot.ext=gtf,isGTFAnnotationFile=T)
summary(fc)

5fa1c8a5218d7e600ba81e6e311584e.png

In [ ]:
#Take a look at the dimensions to see the number of genes
dim(fc$counts)
#Take a look at the first 6 lines
head(fc$counts)

7b644e69a845ed3f7dad35e4c16e929.png