R语言DESeq 包介绍
- 格式:docx
- 大小:306.12 KB
- 文档页数:8
【转录组⼊门】7:差异基因分析作业要求:使⽤R语⾔,载⼊表达矩阵,然后设置好分组信息,统⼀⽤DEseq2进⾏差异分析,当然也可以⾛⾛edgeR或者limma的voom流程。
基本任务是得到差异分析结果,进阶任务是⽐较多个差异分析结果的异同点。
【1】安装DESeq21 # 下⾯是在R语⾔中操作2 # 载⼊安装⼯具3 > source("/biocLite.R")4 # 安装包5 > biocLite("DESeq2")6 # 载⼊包7 > library("DESeq2")DESeq2对于输⼊数据的要求:1.DEseq2要求输⼊数据是由整数组成的矩阵。
2.DESeq2要求矩阵是没有标准化的。
【2】DESeq2进⾏差异表达分析DESeq2分析差异表达基因简单来说只有三步:构建dds矩阵,标准化,以及进⾏差异分析。
# dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) #~在R⾥⾯⽤于构建公式对象,~左边为因变量,右边为⾃变量。
# dds <- DESeq(dds) #标准化# res <- results(dds, contrast=c("condition","treated","control")) #差异分析结果【3】构建dds矩阵1 > library(DESeq2) # 加载包2 > countData <- raw_count_filter[2:7] # 中括号中的数量要与condition中数量⼀致3 > condition <- factor(c("control","control","control","hypoxia","hypoxia","hypoxia"))4 > colData <- data.frame(s=colnames(countData),condition)5 # raw_count_filter:是所有样品的count按照gene id融合后⽣成的矩阵。
Package‘dupree’October13,2022Type PackageTitle Identify Duplicated R Code in a ProjectVersion0.3.0Author Russ Hyde,University of GlasgowMaintainer Russ Hyde<************************>Description Identifies code blocks that have a high level of similaritywithin a set of Rfiles.URL https:///russHyde/dupreeBugReports https:///russHyde/dupree/issuesLicense MIT+file LICENSEEncoding UTF-8Language en-GBLazyData trueSuggests testthat(>=2.1.0),knitr,rmarkdownImports dplyr,purrr,tibble,magrittr,methods,stringdist(>=0.9.5.5),lintr(>=2.0.0),rlangRoxygenNote7.1.0Collate'utils.R''dupree.R''dupree_classes.R''dupree_data_validity.R''dupree_code_enumeration.R''dups-class.R'NeedsCompilation noRepository CRANDate/Publication2020-04-2110:20:02UTCR topics documented:as.data.frame.dups (2)as_tibble.dups (2)dupree (3)12as_tibble.dups dupree_dir (4)dupree_package (5)EnumeratedCodeTable-class (5)print.dups (6)Index7 as.data.frame.dups as.data.frame method for‘dups‘classDescriptionas.data.frame method for‘dups‘classUsage##S3method for class dupsas.data.frame(x,...)Argumentsx any R object....additional arguments to be passed to or from methods.as_tibble.dups convert a‘dups‘object to a‘tibble‘Descriptionconvert a‘dups‘object to a‘tibble‘Usage##S3method for class dupsas_tibble(x,...)Argumentsx A data frame,list,matrix,or other object that could reasonably be coerced to a tibble....Other arguments passed on to individual methods.dupree3 dupree Detect code duplication between the code-blocks in a set offilesDescriptionThis function identifies all code-blocks in a set offiles and then computes a similarity score between those code-blocks to help identify functions/classes that have a high level of duplication,and could possibly be refactored.Usagedupree(files,min_block_size=40,...)Argumentsfiles A set offiles over which code-duplication should be measured.min_block_size dupree uses a notion of non-trivial symbols.These are the symbols/code-words that remain afterfiltering out really common symbols like<-,,,etc.Afterfiltering out these symbols from each code-block,only those blocks containingat least min_block_size symbols are used in the inter-block code-duplicationmeasurement....Unused at present.DetailsCode-blocks under a size threshold are disregarded before analysis(the size threshold is controlled by min_block_size);and only top-level code blocks are considered.Every sufficiently large code-block in the inputfiles will be present in the results at least once.If code-block X and code-block Y are present in a row of the resulting data-frame,then either X is the closest match to Y,or Y is the closest match to X(or possibly both)according to the similarity score;as such,some code-blocks may be present multiple times in the results.Similarity between code-blocks is calculated using the longest-common-subsequence(lcs)mea-sure from the package stringdist.This measure is applied to a tokenised version of the code-blocks.That is,each function name/operator/variable in the code blocks is converted to a unique integer so that a code-block can be represented as a vector of integers and the lcs measure is applied to each pair of these vectors.ValueA tibble.Each row in the table summarises the comparison between two code-blocks(block’a’and block’b’)in the inputfiles.Each code-block in the pair is indicated by:i)thefile(file_a/ file_b)that contains it;ii)its position within thatfile(block_a/block_b;1being thefirst code-block in a givenfile);and iii)the line where that code-block starts in thatfile(line_a/line_b).The pairs of code-blocks are ordered by decreasing similarity.Any match that is returned is either the top hit for block’a’or for block’b’(or both).4dupree_dirExamples#To quantify duplication between the top-level code-blocks in a fileexample_file<-system.file("extdata","duplicated.R",package="dupree")dup<-dupree(example_file,min_block_size=10)dup#For the block-pair with the highest duplication,we print the first four#lines:readLines(example_file)[dup$line_a[1]+c(0:3)]readLines(example_file)[dup$line_b[1]+c(0:3)]#The code-blocks in the example file are rather small,so if# min_block_size is too large,none of the code-blocks will be analysed#and the results will be empty:dupree(example_file,min_block_size=40)dupree_dir Run duplicate-code detection over all R-files in a directoryDescriptionRun duplicate-code detection over all R-files in a directoryUsagedupree_dir(path=".",min_block_size=40,filter=NULL,...,recursive=TRUE)Argumentspath A directory(By default the current working directory).Allfiles in this directory that have a".R",".r"or".Rmd"extension will be checked for code duplication.min_block_size dupree uses a notion of non-trivial symbols.These are the symbols/code-words that remain afterfiltering out really common symbols like<-,,,etc.Afterfiltering out these symbols from each code-block,only those blocks containingat least min_block_size symbols are used in the inter-block code-duplicationmeasurement.filter A pattern for use in grep-this is used to keep only particularfiles:eg,filter= "classes"would comparefiles with‘classes‘in thefilename ...Further arguments for grep.For example,‘filter="test",invert=TRUE‘would disregard allfiles with‘test‘in thefile-path.recursive Should we considerfiles in subdirectories as well?dupree_package5See Alsodupreedupree_package Run duplicate-code detection over allfiles in the‘R‘directory of apackageDescriptionThe function fails if the path does not look like a typical R package(it should have both an R/ subdirectory and a DESCRIPTIONfile present).Usagedupree_package(package=".",min_block_size=40)Argumentspackage The name or path to the package that is to be checked(By default the current working directory).min_block_size dupree uses a notion of non-trivial symbols.These are the symbols/code-words that remain afterfiltering out really common symbols like<-,,,etc.Afterfiltering out these symbols from each code-block,only those blocks containingat least min_block_size symbols are used in the inter-block code-duplicationmeasurement.See AlsodupreeEnumeratedCodeTable-classAn S4class to represent the code blocks as strings of integersDescriptionAn S4class to represent the code blocks as strings of integersSlotsblocks A tbl_df with columns‘file‘,‘block‘,‘start_line‘and‘enumerated_code‘6print.dups print.dups print method for‘dups‘classDescriptionprint method for‘dups‘classUsage##S3method for class dupsprint(x,...)Argumentsx an object used to select a method....further arguments passed to or from other methods.Indexas.data.frame.dups,2as_tibble.dups,2dupree,3dupree_dir,4dupree_package,5 EnumeratedCodeTable-class,5print.dups,67。
微生物群落组装过程 r语言微生物群落组装是指不同微生物个体在特定环境中相互作用并形成群落的过程。
R语言是一种功能强大的统计分析和可视化工具,可以用于微生物群落组装的研究。
以下是使用R语言进行微生物群落组装分析的一般步骤:1. 数据准备:将微生物群落数据整理成适合R语言处理的格式。
常见的数据格式包括OTU表(含有不同样本中各个OTU的丰度信息)和环境因子表(包含与微生物群落组成相关的环境因子数据)。
2. 数据预处理:根据需要对数据进行预处理,如去除低丰度的OTU、归一化样本间的差异等。
这可以使用R中的统计分析包进行操作,例如"vegan"、"phyloseq"等。
3. 多样性分析:使用R中的多样性分析包,如"vegan"和"phyloseq",计算微生物群落的Alpha多样性(单样本内的物种多样性)和Beta多样性(样本间的物种多样性)。
4. 群落结构分析:使用R中的群落结构分析包,如"vegan"和"ade4",可以进行聚类分析、主坐标分析(PCoA)等,以揭示样本间的群落结构差异。
5. 群落组装模型:使用R中的群落组装模型包,如"DESeq2"、"LEfSe"等,可以进行微生物群落组装的相关统计分析,并鉴定与环境因子显著相关的OTU或菌群。
6. 可视化:使用R中的可视化包,如"ggplot2"和"phyloseq",绘制各种统计图表,如柱状图、热图、PCoA图等,以展示微生物群落的组成和结构。
需要注意的是,微生物群落组装是一个复杂的过程,涉及到多个环境因子和微生物个体之间的相互作用。
正确的数据处理和分析方法对于结果的解释和解读至关重要。
因此,在使用R语言进行微生物群落组装研究时,建议参考相关文献和专业指导,确保分析方法的准确性和可靠性。
【R】R语⾔常⽤包汇总⼀、⼀些函数包⼤汇总时间上有点过期,下⾯的资料供⼤家参考基本的R包已经实现了传统多元统计的很多功能,然⽽CRNA的许多其它包提供了更深⼊的多元统计⽅法,下⾯要综述的包主要分为以下⼏个部分:1)多元数据可视化(Visualising multivariate data)绘图⽅法 基本画图函数(如:pairs()、coplot())和 lattice包⾥的画图函数(xyplot()、splom())可以画成对列表的⼆维散点图,3维密度图。
car 包⾥的scatterplot.matrix()函数提供更强⼤的⼆维散点图的画法。
cwhmisc包集合⾥的cwhplot包的pltSplomT()函数类似pair()画散点图矩阵,⽽且可以在对⾓位置画柱状图或密度估计图。
除此之外,scatterplot3d包可画3维的散点图,aplpack包⾥bagplot()可画⼆变量的boxplot,spin3R()可画可旋转的三维点图。
misc3d包有可视化密度的函数。
YaleToolkit包提供许多多元数据可视化技术,agsemisc也是这样。
更特殊的多元图包括:aplpack包⾥的faces()可画Chernoff’s face;MASS包⾥的parcoord()可画平⾏坐标图(矩阵的每⼀⾏画⼀条线,横轴表⽰矩阵的每列); graphics包⾥的stars()可画多元数据的星状图(矩阵的每⼀⾏⽤⼀个星状图表⽰)。
ade4包⾥的mstree()和vegan包⾥的spantree()可画最⼩⽣成树。
calibrate包⽀持双变量图和散点图,chplot包可画convex hull图。
geometry包提供了和qhull库的接⼝,由convexhulln()可给出相应点的索引。
ellipse包可画椭圆,也可以⽤plotcorr()可视化相关矩阵。
denpro包为多元可视化提供⽔平集树形结构(level set trees)。
基因组学数据分析中差异表达分析的使用方法差异表达分析是基因组学研究中常用的一种分析方法,用于比较不同条件或样本之间基因表达的差异。
这种分析方法可以帮助研究人员理解基因表达调控的机制、寻找与特定疾病相关的基因、发现新的分子标记物等。
本文将介绍差异表达分析的基本流程和相关的工具。
差异表达分析的基本流程通常包括:数据预处理、差异基因筛选和分析结果的验证。
首先,进行数据预处理,主要包括原始数据的读取、质量控制和归一化等。
其次,进行差异基因筛选,通过统计学方法或基因表达差异的幅度和显著性来确定差异基因。
最后,对选定的差异基因进行验证,可以使用实验室实验技术如RT-PCR、Western blot 等,或者进行功能富集分析、通路分析等。
在进行差异表达分析时,研究人员通常会使用一些专门的工具和软件。
下面以DESeq2为例,介绍差异表达分析的具体步骤。
首先,将原始数据导入到R或Python等编程环境中,利用相应的包或库进行数据预处理。
DESeq2是R语言中最常用的差异表达分析工具之一。
它可以处理RNA-seq数据,对基因表达的差异进行统计显著性分析。
使用DESeq2进行数据预处理时,首先需要进行数据的归一化处理,通常使用RPKM(reads per kilobase of exon per million mapped reads)或TPM(Transcripts Per Kilobase Million)方法将数据进行归一化处理。
接下来,进行差异基因筛选。
DESeq2使用负二项分布模型对差异表达进行建模,通过p值和调整的p值来判断差异的显著性。
DESeq2还结合了前沿理论来估计基因表达的离散性,提高了差异表达分析的准确性。
通过设定阈值,筛选出显著差异表达的基因。
此外,DESeq2还可以计算基因的表达倍数差异,帮助研究人员判断差异的生物学意义。
最后,对差异基因进行验证。
一种常用的验证方法是使用实验室技术如RT-PCR,验证差异基因在不同条件或样本中的表达水平。
转录组分析_R范文转录组分析是利用高通量测序技术对细胞或组织的所有转录本进行检测和定量的一种方法。
由于转录组包含了所有的转录本信息,因此转录组分析可以帮助我们了解基因的表达模式、功能注释以及基因调控网络等重要信息。
在转录组研究中,R语言是一个广泛应用的工具,用于数据分析和可视化。
下面是一些常见的R语言在转录组分析中的应用:1. 数据预处理:在进行转录组分析之前,我们需要对原始测序数据进行预处理。
这包括质量控制、去除低质量读段和适配序列以及对数据进行过滤和修剪。
R语言提供了许多用于转录组数据预处理的包,例如“ShortRead”和“fastqcr”等。
2. 差异表达基因分析:差异表达基因分析是转录组研究中最常见的任务之一、通过比较不同条件下的基因表达水平,我们可以识别特异表达的基因,并了解它们在生物过程中的功能。
在R语言中,可以使用各种包进行差异表达基因分析,例如“DESeq2”和“limma”等。
3. 功能注释和富集分析:功能注释和富集分析帮助我们理解差异表达基因的生物学功能。
R语言提供了许多包用于进行功能注释和富集分析,例如“clusterProfiler”和“enrichR”等。
这些包可以根据基因的注释信息,并通过各种数据库进行富集分析,例如基因本体论(Gene Ontology,GO)和KEGG等。
4. 转录因子识别和调控网络分析:在转录组分析中,我们还可以通过分析差异表达基因的共享调控因子来研究基因调控网络。
R语言提供了一些用于识别转录因子结合位点的包,例如“ChIPseeker”和“findPeaks”等。
此外,通过整合转录组数据和转录因子结合位点信息,可以构建基因调控网络并进行网络分析。
5. 可视化:数据可视化对于转录组分析来说非常重要,可以帮助我们更直观地理解数据和结果。
R语言提供了许多绘图包用于转录组数据的可视化,例如“ggplot2”和“pheatmap”等。
这些包可以生成各种图表,包括热图、散点图和折线图等。
r语言数据标准化在数据分析中,数据标准化是一个非常重要的步骤。
它可以帮助我们消除不同变量之间的量纲影响,使得数据更加可比较和可解释。
在R语言中,数据标准化可以通过多种方法实现,本文将介绍一些常用的数据标准化方法,并给出在R语言中的具体实现步骤。
一、Z-score标准化。
Z-score标准化是一种常用的数据标准化方法,它通过将原始数据减去均值,再除以标准差的方式来实现。
在R语言中,可以使用如下代码来对数据进行Z-score标准化:```R。
# 假设data为需要标准化的数据框。
data_standardized <scale(data)。
```。
这段代码中,scale()函数可以对数据进行Z-score标准化,得到的data_standardized即为标准化后的数据。
二、Min-max标准化。
Min-max标准化是另一种常用的数据标准化方法,它通过将原始数据减去最小值,再除以最大值和最小值之差的方式来实现。
在R语言中,可以使用如下代码来对数据进行Min-max标准化:```R。
# 假设data为需要标准化的数据框。
data_standardized <apply(data, 2, function(x) (x min(x)) / (max(x) min(x)))。
```。
这段代码中,apply()函数可以对数据进行Min-max标准化,得到的data_standardized即为标准化后的数据。
三、小数定标标准化。
小数定标标准化是另一种常用的数据标准化方法,它通过将原始数据除以一个适当的基数(10的幂次方)来实现。
在R语言中,可以使用如下代码来对数据进行小数定标标准化:```R。
# 假设data为需要标准化的数据框。
data_standardized <apply(data, 2, function(x) x /10^ceiling(log10(max(abs(x)))))。
```。
dep包蛋白组学一、引言随着生物信息学的发展,dep包作为一种高效、精确的蛋白质组学数据分析工具,逐渐在科研领域崭露头角。
本文将从dep包的作用、基本原理、应用案例、优缺点分析、使用方法与技巧等方面进行介绍,以期帮助读者更好地应用dep包开展蛋白质组学研究。
dep包,全称为“DESeq2”,是由Rossettini等人于2012年开发的一种基于RNA测序数据的差异表达分析工具。
近年来,随着质谱技术在蛋白质组学领域的广泛应用,dep包也被成功地应用于蛋白质组学数据分析。
通过对蛋白质表达量的精确测量和差异分析,dep包为研究者揭示了生物过程中的分子机制,为疾病诊断、治疗和预防提供了新的思路。
二、dep包的基本原理dep包基于RNA测序数据进行蛋白质组学分析,其主要原理包括以下两个方面:1.算法原理:dep包采用了一种基于负log似然估计的方法(LDA),将原始测序数据转化为蛋白质表达量估计值。
通过对比不同样本间的表达量差异,dep包可以识别出差异表达的蛋白质。
2.数据处理流程:首先,将原始质谱数据进行预处理,包括峰识别、峰匹配、定量等步骤;然后,利用dep包进行差异表达分析;最后,对分析结果进行统计检验和可视化展示。
三、dep包在蛋白组学分析中的应用案例1.蛋白质定量:dep包可以对蛋白质进行定量分析,从而为研究者提供有关生物过程中蛋白质表达变化的信息。
2.差异表达分析:通过对比不同样本或处理组的蛋白质表达量,dep包可以识别出差异表达的蛋白质,为进一步的生物实验和功能研究提供依据。
3.蛋白质相互作用网络构建:dep包可以结合其他生物信息学工具,如Cytoscape等,构建蛋白质相互作用网络,揭示生物体内的分子调控机制。
四、dep包的优缺点分析1.优点:(1)高效:dep包可以快速处理大量蛋白质组数据;(2)精确:dep包对蛋白质表达量的估计具有较高的准确性;(3)可靠性:dep包的结果经过统计检验,具有较高的可靠性。
Bioconductor包的使用方法1. 什么是Bioconductor包?Bioconductor是一个用于生物信息学和计算生物学的开源软件项目,旨在提供生物学数据的分析和可视化工具。
Bioconductor项目提供了大量的R语言包,用于处理、分析和可视化生物学数据。
这些包涵盖了多个领域,包括基因表达、基因组学、蛋白质组学、代谢组学等。
Bioconductor包的使用方法可以帮助生物学家、生物信息学家和计算生物学家更好地利用Bioconductor项目提供的工具和资源,进行生物学数据的分析和解释。
2. 安装Bioconductor包要使用Bioconductor包,首先需要安装Bioconductor。
安装Bioconductor的方法如下:if (!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")BiocManager::install()上述代码会安装BiocManager包,然后使用BiocManager::install()函数安装Bioconductor。
3. 查找Bioconductor包Bioconductor提供了一个网站,可以用于查找和浏览Bioconductor包。
该网站的网址是,可以在搜索框中输入关键词来查找相关的包。
在R中,可以使用BiocManager::available()函数来列出所有可用的Bioconductor 包。
可以使用BiocManager::search()函数来搜索包含指定关键词的包。
# 列出所有可用的Bioconductor包BiocManager::available()# 搜索包含指定关键词的Bioconductor包BiocManager::search("gene expression")4. 安装Bioconductor包要安装Bioconductor包,可以使用BiocManager::install()函数。
导出差异基因的表达矩阵r语言在R语言中,差异基因表达分析是生物信息学领域的常见任务。
差异基因表达矩阵通常包含基因在不同条件下的表达水平,并且已经过统计测试来识别表达显著差异的基因。
以下是使用R 语言导出差异基因表达矩阵的基本步骤:1.数据导入与预处理:首先,你需要将基因表达数据导入R中。
这些数据通常以表格形式存在,如CSV或TSV文件。
你可以使用read.csv()或read.delim()函数来导入数据。
2.数据格式化:确保数据以适当的格式进行组织,通常是一个矩阵或数据框(data.frame),其中行代表基因,列代表样本或条件。
3.差异表达分析:使用适当的统计方法(如t检验、DESeq2、edgeR等)来分析基因在不同条件下的表达差异。
这些方法通常需要你将数据分组为比较的条件。
4.提取差异基因:根据统计测试的结果,提取显著差异表达的基因。
这通常涉及设置一个显著性阈值(如p值<0.05)和/或一个表达变化阈值(如log2倍数变化>1或<-1)。
5.导出差异基因表达矩阵:创建一个新的数据框或矩阵,其中包含显著差异基因的表达数据,并使用write.csv()或类似的函数将其导出为CSV文件。
以下是一个简化的示例代码,展示了这些步骤的基本框架:R# 1. 数据导入expression_data <- read.csv("gene_expression_data.csv", s = 1)# 2. 数据预处理(如果需要)# ...# 3. 差异表达分析(这里使用假设的函数进行演示)# 假设你已经有了条件分组信息,如条件A和条件Bcondition <- factor(c(rep("A", 3), rep("B", 3))) # 示例条件向量diff_expression_results <-perform_differential_expression_analysis(expression_data, condition)# 4. 提取差异基因# 假设diff_expression_results包含p值和log2倍数变化等信息significant_genes <-diff_expression_results[diff_expression_results$p_value < 0.05 &abs(diff_expression_results$log2FC) > 1,]# 5. 导出差异基因表达矩阵# 这里我们只导出显著差异基因的原始表达数据significant_gene_expression <-expression_data[rownames(significant_genes),]write.csv(significant_gene_expression,"significant_gene_expression_matrix.csv")请注意,上面的代码是一个框架示例,并不包含实际的差异表达分析函数(如perform_differential_expression_analysis),因为这需要根据你的具体数据和所使用的包来定制。
R语言DESeq包介绍 分析RNA序列数据的一个主要任务是探测基因的差异表达,DESeq包提供了测试差异表达的方法,应用负二项分布和收缩的分布方程估计。
1. 包的安装 输入如下命令,DESeq和相关的包就可以自动下载和安装。 > source("http://bioconductor.org/biocLite.R") > biocLite("DESeq") 相关的包会自动下载安装,安装的包如下:
中间会有个选择需要更新相关的包,选择更新全部,更新的包有: 另外还要安装一个数据包,供下面介绍包中的方法时使用,包名为pasilla. 2. 输入数据和准备 2.1 计数表 数据表的第i行第j列元素表示第j个样本的第i个基因有多少个reads。本文使用的数据来自于pasilla数据包,函数system.file告诉我们数据文件保存的路径。 > datafile = system.file( "extdata/pasilla_gene_counts.tsv", package="pasilla" ) > datafile [1] "D:/Program Files/R/R-2.15.3/library/pasilla/extdata/pasilla_gene_counts.tsv" 在R中读取这个文件,使用read.table函数。 > pasillaCountTable = read.table( datafile, header=TRUE, row.names=1 ) > head( pasillaCountTable )
2.2 元数据 没有元数据的数据是没有用的,元数据可以分为三组,分别是样本(行),特征(列)和整个实验的信息。 首先需要样本的描述信息,data.frame的列表示各种信息,行表示7个样本。 > pasillaDesign = data.frame( + row.names = colnames( pasillaCountTable ), + condition = c( "untreated", "untreated", "untreated", + "untreated", "treated", "treated", "treated" ), + libType = c( "single-end", "single-end", "paired-end", + "paired-end", "single-end", "paired-end", "paired-end" ) ) > pasillaDesign
这边简单地使用R代码进行设定,通常情况是从扩展表中读取这些数据。为了分析这些样本,我们需要解释single-end和paired-end的方法,这边首先简单的分析paired-end样本。 > pairedSamples = pasillaDesign$libType == "paired-end" > countTable = pasillaCountTable[ , pairedSamples ] > condition = pasillaDesign$condition[ pairedSamples ] 现在,我们有下面的数据输入 > head(countTable) > condition 对于自己的数据,可以简单的创建因子。 > #not run > condition = factor( c( "untreated", "untreated", "treated", "treated" ) ) 我们现在举例CountDataSet,DESeq包的核心数据结构 > library( "DESeq" ) > cds = newCountDataSet( countTable, condition )
2.3 规范化 函数estimateSizeFactors估计统计数据的大小因子 > cds = estimateSizeFactors( cds ) > sizeFactors( cds )
如果统计数据的每列除以这列的大小因子,这样统计值就变成同一规模,使它们具有可比性。函数counts可以做这个计算。 > head( counts( cds, normalized=TRUE ) )
3. 方差估计 DESeq推断依靠估计典型的数据间的方差和平均关系,或者等效的离差和均值。离差可以理解为生物变异系数的平方。估计离差可以使用以下命令。 > cds = estimateDispersions( cds ) 函数estimateDispersions做了三步,首先估计每条基因的离差,然后,通过估计匹配一条曲线,最后,每个基因分配一个离差,从每条基因估计值和匹配值选一个。为了让用户知道中间过程,fitInfo对象被储存下来。 > str( fitInfo(cds) ) 函数 plotDispEsts可以画出每条基因的估计值和平均正常统计值的关系。 > plotDispEsts( cds )
图1 经验的(黑点)和匹配的(红线)离差值与平均正常统计值的关系图 在任何情况下,可以被子序列测试使用的离差值被存储在cds的特征数据集里。 > head( fData(cds) )
4. 推断:称之为差异表达 4.1 两个实验条件下的标准比较 为了看在条件”untreated"和”treated"下是否有差异表达,我们简单的使用nbinomTest函数, > res = nbinomTest( cds, "untreated", "treated" ) > head(res) 各列表示的是: id feature identier baseMean mean normalised counts, averaged over all samples from both conditions baseMeanA mean normalised counts from condition A baseMeanB mean normalised counts from condition B foldChange fold change from condition A to B log2FoldChange the logarithm (to basis 2) of the fold change pval pvalue for the statistical signicance of this change padj pvalue adjusted for multiple testing with the Benjamini-Hochberg procedure (see the R functionp.adjust), which controls false discovery rate (FDR) 我们首先画log2折叠变换和平均正常统计量的关系,红色的点表示在10%FDR的基因。 > plotMA(res)
图2 log2折叠变换和平均正常统计量的关系图 图3 从nbinomTest的p值统计直方图 > hist(res$pval, breaks=100, col="skyblue", border="slateblue", main="") 我们可以通过FDR过滤有效的基因。 > resSig = res[ res$padj < 0.1, ] 列举最有效的差异表达基因 > head( resSig[ order(resSig$pval), ] ) 我们也想看到有效基因的最强下调。 > head( resSig[ order( resSig$foldChange, -resSig$baseMean ), ] )
或者是最强上调 > head( resSig[ order( -resSig$foldChange, -resSig$baseMean ), ] )
可以保存这些输出到指定CSV文件,命令如下 > write.csv( res, file="My Pasilla Analysis Result Table.csv" ) 对于弱表达基因,需要加更强的改变,这边以两个untreated基因为例。 > ncu = counts( cds, normalized=TRUE )[ , conditions(cds)=="untreated" ] Ncu这边是两列的矩阵 >plotMA(data.frame(baseMean = rowMeans(ncu), + log2FoldChange = log2( ncu[,2] / ncu[,1] )), + col = "black") 画出如图4所示的图像。
图4 两个untreated复制平均表达长度的log2折叠变换图 4.2 部分不复制的情况 如果只有一种条件的基因复制,而另一种条件下的不复制。那么只有复制的条件的可以用来估计离差。为了证明,我们截取仅有三个样本的数据对象。 > cdsUUT = cds[ , 1:3] > pData( cdsUUT ) 图5 “treated" vs. “untreated"对比的MvA图,用两个treated和一个untreated样本 现在我们可以像刚才那样分析了。 >cdsUUT = estimateSizeFactors( cdsUUT ) >cdsUUT = estimateDispersions( cdsUUT ) >resUUT = nbinomTest( cdsUUT, "untreated", "treated" ) 生成类似的图,用 >plotMA(resUUT)
4.3 没有任何复制的情况 DESeq包可以在没有复制的情况分析方法如下。 > cds2 = cds[ ,c( "untreated3", "treated3" ) ] > cds2 = estimateDispersions( cds2, method="blind", sharingMode="fit-only" ) > res2 = nbinomTest( cds2, "untreated", "treated" ) > plotMA(res2)