基因表达量计算流程.md

基因表达量计算流程.md

Administrator 87 2020-07-08

基因表达量计算流程

表达量的定义

一个细胞中或者一定含量的RNA中，含有转录本的数目就是表达量

检测方法原理实际应用
qPCR（实时荧光定量）PCR扩增过程中增加荧光信号的实时监测，看看每个循环产物数量的变化常规基因表达验证（低通量）
Northern blotting放射性标记探针杂交半定量分析（用于少量基因的定性）（低通量）
FISH（荧光原位杂交）在组织水平上利用目的蛋白的抗体检测表达某个基因在特定组织中表达（低通量）
SAGE（基因表达系列分析）分析大量的EST（表达序列标签）寻找不同丰度的标签序列分析样本量较多（中通量）
Microarray（表达芯片或者微阵列）将寡核苷酸探针固定在芯片上，然后将待测样本mRNA加上荧光标记与芯片杂交，然后通过分析荧光信号来监测表达量高通量范围
RNA-seq直接就是高通量测序，将个体的RNA序列测出然后用比对的方法鉴定表达量高通量范围

RNA-seq表达量来源 大体是这样的：

  1. 提取mRNA

  2. 建库测序

  3. 拿到一定的原始数据量raw data【双端x测序深度x建库大小，比如测序深度20X，双端150bpillumina测序，得到的数据就是2x150bpx20M=6G的数据量，但这个6G是碱基量，和硬盘占用空间的gigabytes(g)不一样哦

  4. 质控过滤

  5. clean data

  • 研究物种是没有参考基因组的：从头组装转录本（de novo assembly）【工具包括：Trinity（可以预测更多的基因、转录本）；SOANdenovo-Trans（更好检测高表达转录本）；Oases（有效检测低表达基因以及长的亚型，N10-N50值高）】，然后再将reads比对到重装的转录本上，利用RSEM/Deseq2可以计算表达量

  • 研究物种是有参的，参考序列可以选基因组或者转录本，然后用featureCounts（推荐）/RSEM/htseq定量

  • 如果选基因组为参考，那么就要考虑到基因组含有的内含子和外显子。我们知道，转录的过程其实是把内含子去除了的，因此我们使用比对软件就要支持跨外显子拼接（exon-exon juction） 。最新的当属"Hisat2+Stringtie"，stringtie的merge算法会按照去掉内含子的基因组把所有样本的组装结果合并。

  • 如果选转录本为参考（前提：一般有基因组但不完整），那么转录本已经是去除内含子的序列了，因此计算压力就没那么大，比如常用的BWA、Bowtie2、STAR都可以；或者可以用非比对定量软件Sailfish【它原理是利用参考序列和k-mer长度建立索引，然后将reads直接比对到转录本，然后采用EM(Expectation maximazation，最大期望)算法估计转录本与基因极大似然估计】，这样只能对已知的转录本、基因定量，不能预测新的。

    差异分析三R包法宝：Limma、DESeq2、edgeR

    • 源自芯片的金标准Limma：芯片数据普遍认为符合正态分布 ，而正态分布的群体一般就是用t检验（两个样本）或者方差分析（多个样本）。芯片数据量是很大的，limma据此而生，并且支持多重检验校正来控制假阳性。limma采用贝叶斯模型（Empirical Bayesian model），更新的limma-voom适配了转录组数据

    关于二项分布、泊松分布、正态分布的关系：

    泊松分布，二项分布都是离散分布；正态分布是连续分布，方差分析其实就是看组内方差显著性与组间方差的对比，据此判断不同组（也就是不同样本）之间差异是否存在。

    RNA数据普遍认为符合泊松分布

    DESeq2：采用负二项分布算法（negative binomial distribution）

    RNA-seq中，技术误差是满足泊松分布的，因为期望和方差差不多。但是生物学重复之间的误差不能用泊松分布来描述，因为他的方差可能很大，所以要用负二项分布

    functionpackageframeworkoutputDESeq2 input function
    summarizeOverlapsGenomicAlignmentsR/BioconductorSummarizedExperimentDESeqDataSet
    featureCountsRsubreadR/BioconductormatrixDESeqDataSetFromMatrix
    tximporttximportR/Bioconductorlist of matricesDESeqDataSetFromTximport
    htseq-countHTSeqPythonfilesDESeqDataSetFromHTSeq
     
    #如果是自己构建矩阵的话
library(DESeq2)
exprSet <- as.matrix(原始表达量矩阵)# 每列表示一个样本，每行表示一个基因，必要时可以每行加上基因名
condition <- factor(c("control","case1","case2")) #表型
dds<- DESeqDataSetFromMatrix(countData = exprSet,colData = data.frame(condition),design = ~condition)
dds <- dds[ rowSums(counts(dds)) > CUT, ] #【选做，表示过滤一部分count值低于CUT数字的数据，count自定义】
dds2<-DESeq(dds) #开始运行
result <- results() #查看结果
mcols(result, use.names = TRUE) #可以看每一个结果的意义，利用其中的指标可以画图
summary(result) #看总体结果，多少上调，多少下调【一般认为foldChange绝对值 > 1以及p value < 0.05就是差异基因】

DEseq2经过三部计算：estimateSizeFactors、estimateDispersons、nbinomWaldTest

  • DESeq2用自己的算法进行归一化：
    第一步：estimateSizeFactors计算每个样本的sizeFactor
    相关的过程是这样的：将原始表达量矩阵先进行log转换，然后计算每个基因在所有样本的均值rowMeans(log(raw_counts)) ，再将要计算的单个样本中每个基因的表达量减去前面的均值，再取中位数，就得到了sizeFactor；
    第二步：有了每个样本的sizeFacor，再将样本原始表达量除以sizeFactor，就是归一化的表达量

  • edgeR:
    使用DEGList读取表达矩阵
    利用(count-per-million，CPM）严格过滤count值低的数据
    calcNormFactors函数使用TMM算法对矩阵标准化
    实验设计矩阵Design matrix，model.matrix(~0+group)
    估算离散值dispersionestimateDisp
    构建比较矩阵makeContrasts 、glmQLFTest
    提取差异基因decideTestsDGE 、glmTreat

生物信息 转录组 Linux 统计