文章目录 生信入门（四）——使用DESeq2进行RNA-seq数据分析一、学习目标二、实验数据1、数据来源2、建模计数数据3、转录本丰度4、salmon定量 三、用tximport导入R1、指定文件位置2、将转录本映射到基因3、tximport命令
今日学习内容DESeq2分析RNA-seq数据
一、学习目标 直观地评估 RNA-seq 数据的质量对 RNA-seq 数据进行基本的差异分析将 RNA-seq 结果与其他实验数据进行比较从 FASTQ 文件中量化转录表达将定量导入 R/Bioconductor执行质量控制和探索性数据分析执行差异表达与其他实验数据重叠构建动态报告整合 DESeq2 和 zinbwave 以获取单细胞 RNA-seq 数据 二、实验数据 1、数据来源
GEO 条目 GSE52778
本文将展示如何使用替代数据集（tximport小插图中使用的tximportData包）导入 RNA-seq 量化数据。后来我们将加载了数道数据集，其使用计数summarizeOverlaps从GenomicAlignments包。如下所述，我们推荐使用tximport管道来生成计数矩阵，但我们还没有包含气道数据集必要量化文件的 Bioconductor 包。 2、建模计数数据作为输入，基于计数的统计方法，例如DESeq2 , edgeR , limma with the voom method , DSS、EBSeq 和baySeq期望输入数据是从 RNA-seq 或其他高通量测序实验中获得的计数矩阵。第i行和j 中的值矩阵的第 -th 列表示样本j 中的基因i分配了多少读数（或片段，对于双端 RNA-seq）。类似地，对于其他类型的检测，矩阵的行可能对应于例如结合区域（使用 ChIP-Seq）、细菌种类（使用宏基因组数据集）或肽序列（使用定量质谱）。
3、转录本丰度在此工作流程中，我们将展示如何使用由Salmon (Patro et al. 2017)软件包量化的转录本丰度。Salmon和其他方法，例如Sailfish 、kallisto 或RSEM ，估计所有（已知的、带注释的）转录本的相对丰度，而没有对齐读取。因为估计转录本的丰度涉及推理步骤，所以估计计数. 大多数方法要么使用称为估计最大化或贝叶斯技术的统计框架来估计丰度和计数。在量化之后，我们将使用tximport （Soneson、Love 和 Robinson 2015）包来组装估计计数和偏移矩阵，以便与 Bioconductor 差异基因表达包一起使用。
将转录本丰度量词与tximport结合使用来生成基因级计数矩阵和归一化偏移的优点是：
这种方法纠正了样本中基因长度的任何潜在变化（例如，来自不同异构体的使用）与基于对齐的方法相比，其中一些方法要快得多，需要更少的内存和磁盘使用量可以避免丢弃那些可以与具有同源序列的多个基因对齐的片段（Robert and Watson 2015）。请注意，转录本丰度量词会跳过存储读取比对（SAM 或 BAM 文件）的大文件的生成，而是生成存储每个转录本的估计丰度、计数和有效长度的较小文件。
salmon量化转录本丰度教程
略~因为我看不懂
量化之后，我们可以使用tximport将数据导入 R 并使用 Bioconductor 包进行统计分析。通常，我们只需将tximport指向quant.sf我们机器上的文件。但是，因为我们将这些文件作为 R 包的一部分进行分发，所以我们必须执行一些额外的步骤，以确定 R 包以及这些文件在您的计算机上的位置。
R 函数system.file可用于找出包中的文件安装在计算机上的哪个位置。这里我们要求提供R 包存储外部数据的目录的完整路径，这是tximportData包的一部分。
我们需要创建一个指向量化文件的命名向量。我们将首先通过读取包含样本 ID 的表来创建文件名向量，然后将其与和结合。（我们对量化文件进行了 gzip 压缩，使数据包更小，这对于我们用于导入文件的 R 函数来说不是问题。）
samples<-read.table(file.path(dir,"samples.txt"),header = TRUE)files<-file.path(dir,"salmon",samples$run,"quant.sf.gz")names(files)<-paste0("samples",1:6)all(file.exists(files)) 2、将转录本映射到基因转录本需要与基因 ID 相关联，以便进行基因级汇总。因此，我们将构建一个名为data.frame 的tx2gene两列：1) 转录本 ID 和 2) 基因 ID。列名无关紧要，但必须使用此列顺序。转录本 ID 必须与丰度文件中使用的相同。这可以通过在下载转录组 FASTA 的同时下载tx2geneGTF 文件并使用 Bioconductor 的TxDb基础设施从 GTF 文件生成来最容易地完成。
创建tx2gene data.frame可以通过在TxDb对象上调用AnnotationDbi包中的select函数来完成。以下代码可用于构造这样的表：
#将转录本映射到基因library("TxDb.Hsapiens.UCSC.hg38.knownGene")txdb<-TxDb.Hsapiens.UCSC.hg38.knownGenetxdbk<-keys(txdb,keytype = "TXNAME")#k 创建tx2gene tx2gene<-select(txdb,k,"GENEID","TXNAME")#使用Gencode V27 chr 转录本构建salmon索引library("readr")tx2gene<-read_csv(file.path(dir,"tx2gene.gencode.v27.csv"))head(tx2gene) 3、tximport命令 #tximport 命令——将salmon转录本量化导入Rlibrary("tximport")library("jsonlite")txi<-tximport(files,type="salmon",tx2gene = tx2gene)#txi 对象是一个矩阵列表names(txi)txi$counts[1:3,1:3]#提取前三行前三列的txi矩阵中的countstxi$length[1:3,1:3]#提取前三行前三列的txi矩阵中的lengthtxi$abundance[1:3,1:3]#提取前三行前三列的txi矩阵中的abundancetxi$countsFromAbundance
如果我们继续使用 GEUVADIS 样本，我们将使用以下代码行创建DESeqDataSet。因为样本之间没有差异（相同的群体和相同的测序批次），我们指定了一个~1的设计公式，这意味着我们只能拟合一个截距项——所以我们不能对这些样本进行差异表达分析。
告一段落，明天——>探索性数据分析