看了老大比较基因组学视频后扩展的小内容

作者: 小梦游仙境 | 来源:发表于2019-11-21 20:57 被阅读0次

老大视频比较基因组在IGV

看了http://www.bio-info-trainee.com/2218.html

wes：探针捕获

rna-seq：不会出现外显子两边下降

chip-seq：不会出现整齐，得到的文件是peaks bed files

RNA-seq这十年

https://mp.weixin.qq.com/s/a3y46NNNO-wardO3XWwh0w

经过不断的技术开发和改进，以Roche公司的454技术、illumina公司的Solexa/HiSeq技术和ABI公司的SOLID技术为标记的第二代测序技术诞生了。
在实验室中，其标准流程就分为三步：

第一步是构建测序文库，这一步骤包括提取RNA，==富集mRNA或清除核糖体RNA==，合成 cDNA，加上接头。

第二步，在高通量平台（通常是Illumina平台）上对文库进行测序，每个样本的测序深度为10-30M读长数（读长这里就是前面说的reads）。

第三步是数据分析，具体的工作是：对测序得到的读长进行比对(aligning)和/或组装到转录组上，对这些覆盖了转录组的读长进行过滤，归一化(Normalization)，根据统计模型找出那些在不同样本之间有差异的转录本。早期的RNA-seq从大量的实验样本中产生了DGE数据，这充分说明了RNA-seq在广泛的生物体以及系统中的使用，这些生物体包括玉米(Zea mays), 拟南芥(Arabiodopsis thaliana), 酿酒酵母（Saccharomyces cerevisae），小鼠(Mus musculus)以及人类。虽然RNA-seq这个术语经常被用于那些完全不同的方法学方法和/或生物学，但是==DGE分析仍然是RNA-seq==（补充材料中的表1）的主要应用，并被视为常规研究工具

illumina测序原理

1）构建DNA测序文库
- 把一堆乱糟糟的DNA分子用超声波打断成一定长度范围的小片段。目前除了一些特殊的需求之外，基本都是打断为300bp-800bp长的序列片段，并在这些小片段的两端添加上不同的接头【注】，构建出单链DNA文库，以备测序之用；
【注】接头在illumina中一般分为P5和P7接头，其中一个带有和flowcell上的探针反向互补的序列，以完成待测序列和探针结合的作用，另外一个接头带有barcord序列以区分不同的样本
测序流动槽（flowcell）
flowcell是用于吸附流动DNA片段的槽道，也是核心的测序反应容器——所有的测序过程就发生在这里。当文库建好后，这些文库中的DNA在通过flowcell的时候会随机附着在flowcell表面的槽道（称为lane）上。每个flowcell有8个lane，每个lane的表面都附有很多接头，这些接头能和建库过程中加在DNA片段两端的接头相互配对，这就是为什么flowcell能吸附建库后的DNA的原因，并能支持DNA在其表面进行桥式PCR的扩增，理论上这些lane之间是不会相互影响的。
3）桥式PCR扩增与变性
这是NGS技术的一个核心特点。桥式PCR以flowcell表面所固定的序列为模板，进行桥形扩增，经过不断的扩增和变性循环，最终每个DNA片段都将在各自的位置上集中成束，每一个束都含有单个DNA模板的很多分拷贝，这一过程的目的在于实现将单一碱基的信号强度进行放大，以达到测序所需的信号要求。
4）测序
- 测序方法采用边合成边测序的方法。向反应体系中同时添加DNA聚合酶、接头引物和带有碱基特异荧光标记的4中dNTP（如同Sanger测序法）。这些dNTP的3’-OH被化学方法所保护，因而每次只能添加一个dNTP，这就确保了在测序过程中，一次只会被添加一个碱基。同时在dNTP被添加到合成链上后，所有未使用的游离dNTP和DNA聚合酶会被洗脱掉。接着，再加入激发荧光所需的缓冲液，用激光激发荧光信号，并有光学设备完成荧光信号的记录，最后利用计算机分析将光学信号转化为测序碱基。这样荧光信号记录完成后，再加入化学试剂淬灭荧光信号并去除dNTP 3’-OH保护基团，以便能进行下一轮的测序反应。
  
  Illumina的这种每次只添加一个dNTP的技术特点能够很好的地解决同聚物长度的准确测量问题，它的主要测序错误来源是碱基的替换，目前它的测序错误率在1%-1.5%左右。测序周期以人类基因组重测序为例，30x-50x测序深度对于Hisq系列需要3-5天时间，而对于2017年初最新推出的NovaSeq系列则只需要40个小时！

陈巍学基因

https://www.bilibili.com/video/av23072634

illumina官网

https://support.illumina.com/documentation.html

image-20190923132407448

Differences Between NGS and Sanger Sequencing

In principle, the concepts behind Sanger vs. next-generation sequencing (NGS) technologies are similar. In both NGS and Sanger sequencing (also known as dideoxy or capillary electrophoresis sequencing), DNA polymerase adds fluorescent nucleotides one by one onto a growing DNA template strand. Each incorporated nucleotide is identified by its fluorescent tag.
The critical difference between Sanger sequencing and NGS is sequencing volume. While the Sanger method only sequences a single DNA fragment at a time, NGS is massively parallel, sequencing millions of fragments simultaneously per run. This high-throughput process translates into sequencing hundreds to thousands of genes at one time. NGS also offers greater discovery power to detect novel or rare variants with deep sequencing.

Differences Between NGS and qPCR

When comparing next-generation sequencing (NGS) vs. qPCR technologies, the key difference is discovery power. While both offer highly sensitive and reliable variant detection, qPCR can only detect known sequences. In contrast, NGS is a hypothesis-free approach that does not require prior knowledge of sequence information. NGS provides higher discovery power to detect novel genes and higher sensitivity to quantify rare variants and transcripts.
NGS vs. qPCR technologies also differ in scalability and throughput. While qPCR is effective for low target numbers, the workflow can be cumbersome for multiple targets. NGS is preferable for studies with many targets or samples. A single NGS experiment can identify variants across thousands of target regions with single-base resolution.

Benefits of RNA-Seq vs. Microarray Technology

Ability to detect novel transcripts: Unlike arrays, RNA-Seq technology does not require species- or transcript-specific probes. It can detect novel transcripts, gene fusions, single nucleotide variants, indels (small insertions and deletions), and other previously unknown changes that arrays cannot detect.1,2
Wider dynamic range: With array hybridization technology, gene expression measurement is limited by background at the low end and signal saturation at the high end. RNA-Seq technology produces discrete, digital sequencing read counts, and can quantify expression across a larger dynamic range (>105 for RNA-Seq vs. 103 for arrays).1,2,3
Higher specificity and sensitivity: Compared to microarrays, RNA-Seq technology can detect a higher percentage of differentially expressed genes, especially genes with low expression.4-6
Simple detection of rare and low-abundance transcripts:** Sequencing coverage depth can easily be increased to detect rare transcripts, single transcripts per cell, or weakly expressed genes.

实时荧光定量 PCR、传统 PCR 与数字 PCR 概述

	数字 PCR	实时荧光定量 PCR	传统 PCR
概述	测量阴性平行样的比例以确定绝对拷贝数。	在发生 PCR 扩增时进行测量。	在 PCR 循环结束时测量积聚的 PCR 产物的量。
定量？	是，阴性 RCR 反应的比例适合泊松统计算法。	是，因为在 PCR 的指数（对数）期内采集数据，在此期间 PCR 产物的数量与核酸模板的量成正比。	否，但是比较凝胶上的扩增条带的强度与已知浓度的标准品可为您提供“半定量”结果。
应用	病毒载荷绝对定量核酸标准品绝对定量下一代测序文库绝对定量稀有等位基因检测基因表达绝对定量混合物富集和分离	基因表达定量芯片验证质量控制和检测验证病原体检测SNP 基因分型拷贝数变异microRNA 分析病毒定量siRNA/RNAi 实验	DNA 扩增用于：测序基因分型分子克隆
总结	数字 PCR 的优点：无需依赖参考物或标准品通过增加 PCR 平行样的总数可实现所需的精度高度耐受抑制剂能够分析复杂混合物与传统 qPCR 不同，数字 PCR 对出现的拷贝数呈线性反应，可以检测到较小倍数变化的差异。	实时荧光定量 PCR 的优点：提高检测的动态范围无需 PCR 后处理检测能够覆盖低至 2 倍的变化在 PCR 的指数期内采集数据报告基团荧光信号的增加与生成的扩增子数量成正比裂解探针提供扩增子的永久裂解记录	传统 PCR 的缺点：精确度差灵敏度低动态范围小 < 2 logs分辨率低非自动化仅限基于规格的区别分析结果未表示为数字使用溴化乙锭染色不足以定量PCR 后处理

三代测序技术比较

	公司	平台名称	测序方法	检测方法	大约读长(碱基数)	优点	相对局限性
第一代	ABI/生命技术公司	3130xL-3730xL	桑格-毛细管电泳测序法	荧光/光学	600-1000	高读长，准确度一次性达标率高，能很好处理重复序列和多聚序列	通量低；样品制备成本高，使之难以做大量的平行测序
第二代	Roche/454	基因组测序仪FLX系统	焦磷酸测序法	光学	230-400	在第二代中最高读长；比第一代的测序通量大	样品制备较难；难于处理重复和同种碱基多聚区域；试剂冲洗带来错误累积；仪器昂贵
第二代	Illumina	HiSeq2000,HiSeq2500/MiSeq	可逆链终止物和合成测序法	荧光/光学	2x150	很高测序通量	仪器昂贵；用于数据删节和分析的费用很高
第二代	ABI/Solid	5500xlSolid系统	连接测序法	荧光/光学	25-35	很高测序通量；在广为接受的几种第二代平台中，所要拼接出人类基因组的试剂成本最低	测序运行时间长；读长短，造成成本高，数据分析困难和基因组拼接困难；仪器昂贵
第三代	太平洋生物科学公司	PacBio RS	实时单分子DNA测序	荧光/光学	~1000	高平均读长，比第一代的测序时间降低；不需要扩增；最长单个读长接近3000碱基	并不能高效地将DNA聚合酶加到测序阵列中；准确性一次性达标的机会低（81-83%）；DNA聚合酶在阵列中降解；总体上每个碱基测序成本高（仪器昂贵）；
第三代	全基因组学公司	GeXP遗传分析系统	复合探针锚杂交和连接技术	荧光/光学	10	在第三代中通量最高；在所有测序技术中，用于拼接一个人基因组的试剂成本最低；每个测序步骤独立，使错误的累积变得最低	低读长；模板制备妨碍长重复序列区域测序；样品制备费事；尚无商业化供应的仪器
第三代	Ion Torrent/生命技术公司	个人基因组测序仪（PGM）	合成测序法	以离子敏感场效应晶体管检测pH值变化	100-200	对核酸碱基的掺入可直接测定；在自然条件下进行DNA合成（不需要使用修饰过的碱基）	一步步的洗脱过程可导致错误累积；阅读高重复和同种多聚序列时有潜在困难

一篇文章说清楚什么是“插入片段”？

huangshujia.me/2018/08/11/2018-08-11-What-is-insert-sequence.html

Transitioning from Microarrays to mRNA-Seq

chrome-extension://cdonnmffkdaoajfknoeeecmchibpmkmg/static/pdf/web/viewer.html?file=https%3A%2F%2Fwww.illumina.com%2Fcontent%2Fdam%2Fillumina-marketing%2Fdocuments%2Ficommunity%2Farticle_2011_12_ea_rna-seq.pdf

While microarrays are effective for identifying the expression of known genes and transcripts, the snapshot of the transcriptome they provide is incomplete. They cannot detect previously unidentified genes or transcripts, and will thus miss changes in their expression that may be associated with a phenotype of interest.

虽然微阵列可以有效地识别已知基因和转录本的表达，但它们提供的转录组快照是不完整的。它们无法检测到以前未被识别的基因或转录本，因此会错过它们表达的变化，而这些变化可能与感兴趣的表型有关。

As the cost of sequencing has gone down, researchers have begun to turn to mRNA sequencing (mRNA-Seq) for gene expression analysis. Offering more comprehensive coverage through ==single-read or paired-end== sequencing, this approach enables rapid profiling and deep investi-gation of the transcriptome.

随着测序成本的降低，研究人员开始转向mRNA测序(mRNA- seq)进行基因表达分析。通过单读或双端测序提供更全面的覆盖，这种方法支持快速分析和对转录组的深入研究。

区分两种测序方式
单端测序(Single-read)首先将DNA样本进行片段化处理形成200-500bp的片段，引物序列连接到DNA片段的一端，然后末端加上接头，将片段固定在flow cell上生成DNA簇，上机测序单端读取序列。
双端测序(Paired-end)方法是指在构建待测DNA文库时在两端的接头上都加上测序引物结合位点，在第一轮测序完成后，去除第一轮测序的模板链，用对读测序模块(Paired-End Module)引导互补链在原位置再生和扩增，以达到第二轮测序所用的模板量，进行第二轮互补链的合成测序。

要用英文描述RNA-seq优于microarray的地方

用limma包的voom方法来做RNA-seq 差异分析

结合老大说的，==关于芯片数据在用limma时新增加的voom算法==找到老大原文出处了

老大原文如下

大家都知道，这十几年来最流行的差异分析软件就是R==的limma包了，但是它以前只支持microarray的表达数据。==

考虑到大家都熟悉了它，它又发了一个==voom的方法，让它从此支持RNA-seq的count数据==啦！

大家都知道芯片数据跟RNA-seq数据的本质就是value的分布不一样，所以各种==针对RNA-seq的差异分析包也就是提出来一个新的normalization方法而已==。

而我们==limma本身就提出了一个voom的方法来对RNA-seq数据进行normalization==

使用这个包也需要三个数据：

表达矩阵
分组矩阵
差异比较矩阵

用法与limma一模一样的，只是多了一个normalization而已。

读取自己的数据

一般我们会自己读取表达矩阵和分组信息，下面是一个例子：

options(warn=-1)
suppressMessages(library(DESeq2))
suppressMessages(library(limma))
suppressMessages(library(pasilla))
data(pasillaGenes)
exprSet=counts(pasillaGenes)
head(exprSet)  ##表达矩阵如下
##             treated1fb treated2fb treated3fb untreated1fb untreated2fb
## FBgn0000003          0          0          1            0            0
## FBgn0000008         78         46         43           47           89
## FBgn0000014          2          0          0            0            0
## FBgn0000015          1          0          1            0            1
## FBgn0000017       3187       1672       1859         2445         4615
## FBgn0000018        369        150        176          288          383
##             untreated3fb untreated4fb
## FBgn0000003            0            0
## FBgn0000008           53           27
## FBgn0000014            1            0
## FBgn0000015            1            2
## FBgn0000017         2063         1711
## FBgn0000018          135          174
(group_list=pasillaGenes$condition)
## [1] treated   treated   treated   untreated untreated untreated untreated
## Levels: treated untreated
##这是分组信息，7个样本，3个处理的，4个未处理的对照！

另外一个例子是从airway里面得到表达矩阵和分组信息

library(airway)
data(airway)
exprSet=assays(airway)$counts  ## 表达矩阵
group_list=colData(airway)$dex ## 分组信息

第一步：构建分组矩阵

suppressMessages(library(limma))
design <- model.matrix(~factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
##              treated untreated
## treated1fb         1         0
## treated2fb         1         0
## treated3fb         1         0
## untreated1fb       1         1
## untreated2fb       1         1
## untreated3fb       1         1
## untreated4fb       1         1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$`factor(group_list)`
## [1] "contr.treatment"

第二步：根据分组信息和表达矩阵进行normalization

这里构建了一个对象 An object of class “EList”

v <- voom(exprSet,design,normalize="quantile")
## 下面的代码如果你不感兴趣不需要运行，免得误导你
## 就是看看normalization前面的数据分布差异
#png("RAWvsNORM.png")
exprSet_new=v$E
par(cex = 0.7)
n.sample=ncol(exprSet)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
par(mfrow=c(2,2))
boxplot(exprSet, col = cols,main="expression value",las=2)
boxplot(exprSet_new, col = cols,main="expression value",las=2)
hist(exprSet)
hist(exprSet_new)

dev.off()


可以很明显看到normalization前后数据分布差异非常大！

####第三步：做差异分析，提取差异分析结果

> 这里不需要差异比较矩阵了，因为我的分组矩阵表明样本就分成两组，直接提取coef=2的数据即可

```r
fit <- lmFit(v,design)
fit2 <- eBayes(fit)
tempOutput = topTable(fit2, coef=2, n=Inf)
DEG_voom = na.omit(tempOutput)
head(DEG_voom)
##                  logFC  AveExpr         t      P.Value    adj.P.Val
## FBgn0029167  2.1608689 7.880589  41.19142 5.195806e-08 0.0007518331
## FBgn0003137 -0.9560776 8.421903 -29.97211 2.920473e-07 0.0021129620
## FBgn0003748 -1.0385933 6.395990 -23.78787 1.020682e-06 0.0049230892
## FBgn0035085  2.5190255 5.221183  21.01339 1.993995e-06 0.0058809216
## FBgn0050185 -0.4886279 5.487547 -19.95422 2.635106e-06 0.0058809216
## FBgn0015568 -1.1040979 3.922438 -19.96954 2.624231e-06 0.0058809216
##                    B
## FBgn0029167 9.065020
## FBgn0003137 7.800063
## FBgn0003748 6.652695
## FBgn0035085 5.870585
## FBgn0050185 5.734162
## FBgn0015568 5.557560

==而直接用于芯片数据的代码如下==

design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
library(limma)
####构建比较矩阵,设置用来对比的组别
contrast.matrix<-makeContrasts(contrasts=c('tumor-normal'),levels = design)
######limma三部曲,只需要归一化后的数据、实验矩阵、比较矩阵的输入
fit <- lmFit(exprSet,design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1)
####上面得到了p值等统计的结果，topTable对p值校验，对基因排序
DEG_mRNA_dat_filter <- topTable(fit2, adjust="fdr", number=nrow(fit2))

用DESeq2包来对RNA-seq数据进行差异分析

老大原文如下

差异分析的套路都是差不多的，大部分设计思想都是继承limma这个包，DESeq2也不例外。

DESeq2是DESeq包的更新版本，看样子应该不会有DESeq3了，哈哈，它的设计思想就是针对count类型的数据。

可以是任意features的count数据，比如对==各个基因的count，或者外显子，或者CHIP-seq的一些feature，都可以用来做差异分析。==

使用这个包也是需要三个数据：

表达矩阵
分组矩阵
差异比较矩阵

总结起来三个步骤，我下面会一一讲解

重点就是构造一个dds的对象，
然后直接用DESeq函数进行normalization处理即可，
处理之后用results函数来提取差异比较结果

读取自己的数据

一般我们会自己读取表达矩阵和分组信息，下面是一个例子：

options(warn=-1)
suppressMessages(library(DESeq2))
suppressMessages(library(limma))
suppressMessages(library(pasilla))
data(pasillaGenes)
exprSet=counts(pasillaGenes)
head(exprSet)  ##表达矩阵如下
##             treated1fb treated2fb treated3fb untreated1fb untreated2fb
## FBgn0000003          0          0          1            0            0
## FBgn0000008         78         46         43           47           89
## FBgn0000014          2          0          0            0            0
## FBgn0000015          1          0          1            0            1
## FBgn0000017       3187       1672       1859         2445         4615
## FBgn0000018        369        150        176          288          383
##             untreated3fb untreated4fb
## FBgn0000003            0            0
## FBgn0000008           53           27
## FBgn0000014            1            0
## FBgn0000015            1            2
## FBgn0000017         2063         1711
## FBgn0000018          135          174
(group_list=pasillaGenes$condition)
## [1] treated   treated   treated   untreated untreated untreated untreated
## Levels: treated untreated
##这是分组信息，7个样本，3个处理的，4个未处理的对照！

第一步：构建dds对象

那么根据这两个数据就可以构造dds的对象了

colData <- data.frame(row.names=colnames(exprSet), 
                       group_list=group_list
                       )
## 这是一个复杂的方法构造这个对象！
dds <- DESeqDataSetFromMatrix(countData = exprSet,
                              colData = colData,
                              design = ~ group_list)

## design 其实也是一个对象，还可以通过design(dds)来继续修改分组信息，但是一般没有必要。
dds
## class: DESeqDataSet 
## dim: 14470 7 
## exptData(0):
## assays(1): counts
## rownames(14470): FBgn0000003 FBgn0000008 ... FBgn0261574
##   FBgn0261575
## rowData metadata column names(0):
## colnames(7): treated1fb treated2fb ... untreated3fb untreated4fb
## colData names(1): group_list

可以看到我们构造的dds对象有7个样本，共14470features

从基因名可以看出，是果蝇的测序数据

我们也可以直接从expressionSet这个对象构建dds对象！

library(airway)
data(airway)
suppressMessages(library(DESeq2))
dds  <- DESeqDataSet(airway, design = ~ cell+ dex)
design(dds) <- ~ dex  
## 构造好的对象还可以更改分组信息

第二步：做normalization

suppressMessages(dds2 <- DESeq(dds)) 
##直接用DESeq函数即可
## 下面的代码如果你不感兴趣不需要运行，免得误导你
## 就是看看normalization前面的数据分布差异
rld <- rlogTransformation(dds2)  ## 得到经过DESeq2软件normlization的表达矩阵！
exprSet_new=assay(rld)
par(cex = 0.7)
n.sample=ncol(exprSet)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
par(mfrow=c(2,2))
boxplot(exprSet, col = cols,main="expression value",las=2)
boxplot(exprSet_new, col = cols,main="expression value",las=2)
hist(exprSet)
hist(exprSet_new)


#### 第三步：提取差异分析结果

```R
resultsNames(dds2)
## [1] "Intercept"           "group_listtreated"   "group_listuntreated"
## 其实如果只分成了两组，并没有必要指定这个比较矩阵！
res <-  results(dds2, contrast=c("group_list","treated","untreated")) 

## 提取你想要的差异分析结果，我们这里是trt组对untrt组进行比较
resOrdered <- res[order(res$padj),]
resOrdered=as.data.frame(resOrdered)
head(resOrdered)
##              baseMean log2FoldChange     lfcSE      stat        pvalue
## FBgn0039155  453.2753      -4.281830 0.1919977 -22.30146 3.576174e-110
## FBgn0029167 2165.0445      -2.182745 0.1080670 -20.19807  1.017931e-90
## FBgn0035085  366.8279      -2.436860 0.1505280 -16.18875  6.054219e-59
## FBgn0029896  257.9027      -2.511257 0.1823764 -13.76964  3.881667e-43
## FBgn0034736  118.4074      -3.166392 0.2375506 -13.32933  1.562878e-40
## FBgn0040091  610.6035      -1.526400 0.1278555 -11.93848  7.457520e-33
##                      padj
## FBgn0039155 2.764025e-106
## FBgn0029167  3.933795e-87
## FBgn0035085  1.559769e-55
## FBgn0029896  7.500351e-40
## FBgn0034736  2.415897e-37
## FBgn0040091  9.606528e-30

差异分析结果很容易看懂啦！

看了老大比较基因组学视频后扩展的小内容

老大视频比较基因组在IGV

RNA-seq这十年

illumina测序原理

陈巍学基因

illumina官网

Differences Between NGS and Sanger Sequencing

Differences Between NGS and qPCR

Benefits of RNA-Seq vs. Microarray Technology

实时荧光定量 PCR、传统 PCR 与数字 PCR 概述

一篇文章说清楚什么是“插入片段”？

Transitioning from Microarrays to mRNA-Seq

要用英文描述RNA-seq优于microarray的地方

用limma包的voom方法来做RNA-seq 差异分析

读取自己的数据

第一步：构建分组矩阵

第二步：根据分组信息和表达矩阵进行normalization

dev.off()

用DESeq2包来对RNA-seq数据进行差异分析

读取自己的数据

第一步：构建dds对象

第二步：做normalization

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

测序知识

基因组

RNA CHIP ATAC -seq