|
||||||
|
> gtsdata_test=read.table(“gtsdata.txt”, header=T) > gtsenv=read.table(“gtsenv.txt”, header=T) > gtsdata_data_t<-t(gtsdata_data) > decorana(gtsdata_data_t) Call: decorana(veg = gtsdata_data_t) Detrended correspondence analysis with 26 segments. Rescaling of axes with 4 iterations. DCA1 DCA2 DCA3 DCA4 Eigenvalues 0.8634 0.4834 0.23788 0 Decorana values 0.8721 0.3793 0.07223 0 Axis lengths 5.3292 2.1115 1.80907 0 > gts.ca=cca(gtsdata_data_t) > gts.ca Call: cca(X = […] a.str <- matrix(c(‘1′,’2′,’3′,’5′,NA,’6′) + ,c(2,3),dimnames = list(c(‘g1′,’g2′),c(‘t1′,’t2′,’t3′))) a.str # t1 t2 t3 # g1 “1” “3” NA # g2 “2” “5” “6” a.num <- apply(a.str, c(1,2), as.numeric) a.num # t1 t2 t3 # g1 1 3 NA # g2 2 5 6 Note: 第一行，第一列位置要为空！！ # –Dot/Line graphs– # We might be interested in plotting the first two dimensions of a PCoA or # NMDS plot. Let’s do this with data generated in the Costello stool analysis # tutorial. The necessary file is in your folder. nmds<-read.table(file=”stool.final.an.thetayc.0.03.lt.nmds.axes”, header=T) plot(nmds$axis1, nmds$axis2) # or plot(nmds$axis2~nmds$axis1) # Looking at the group names in […] “基于Vegan 软件包的生态学数据排序分析 赖江山 米湘成 (中国科学院植物研究所植被与环境变化国家重点实验室，北京 100093) 摘要：群落学数据一般是多维数据，例如物种属性或环境因子的属性。多元统计分析是群落生态学常用的分析方法，排序（ordination）是多元统计最常用的方法之一。CANOCO是广泛使用的排序软件，但缺点是商业软件价格不菲，版本更新速度也很慢。近年来，R语言以其灵活、开放、易于掌握、免费等诸多优点，在生态学和生物多样性研究领域迅速赢得广大研究人员的青睐。R语言中的外在软件包“Vegan”是专门用于群落生态学分析的工具。Vegan能够提供所有基本的排序方法，同时具有生成精美排序图的功能，版本更新很快。我们认为Vegan包完全可以取代CANOCO，成为今后排序分析的首选统计工具。本文首先简述排序的原理和类型，然后介绍Vegan的基本信息和下载安装过程，最后以古田山24公顷样地内随机抽取40个20m×20m的样方为例，展示Vegan包内各种常用排序方法（PCA,RDA,CA和CCA）和排序图生成过程，希望能为R的初学者尽快熟悉并利用Vegan包进行排序分析提供参考。 gtsdata gtsenv.txt 赖江山.pdf > setwd(“/winxp_disk2/shenzy/R/Vegan”) > gtsdata=read.table(“gtsdata.txt”, header=T) > gtsenv=read.table(“gtsenv.txt”, header=T) > install.packages(“vegan”) Installing package(s) into ‘/home/shenzy/R/x86_64-pc-linux-gnu-library/2.15’ (as ‘lib’ is unspecified) 试开URL’http://cran.csiro.au/src/contrib/vegan_2.0-4.tar.gz’ Content type ‘application/x-gzip’ length 1576584 bytes (1.5 Mb) 打开了URL ================================================== downloaded 1.5 Mb * installing *source* package ‘vegan’ … ** 成功将‘vegan’程序包解包并MD5和检查 ** libs gfortran -fpic -O3 […] source(“http://www.bioconductor.org/biocLite.R”); biocLite(“affy”); biocLite(“Biobase”); library(affy); library(Biobase); >bac_4sampledata=read.csv(“/home/R_heatmap/4sample_R_cluster.csv”, sep=”\t”) > row.names(bac_4sampledata)<-bac_4sampledata$Group > bac_4sample_Datamatrix<-data.matrix(bac_4sampledata[,2:5]) > heatmap.2(bac_4sample_Datamatrix, distfun=dist,col=greenred(256), scale=”row”, key=TRUE, symkey=FALSE, density.info=”none”, trace=”none”, cexRow=0.5, cexCol=0.7,margin=c(7,30), keysize=1.5); 4sample_R_cluster_stdtop100 > heatmap.2(bac_4sample_Datamatrix, distfun = function(x) dist(x,method = ‘euclidean’),hclustfun = function(x) hclust(x,method = ‘centroid’),col=greenred(256), scale=”row”, key=TRUE, symkey=FALSE, density.info=”none”, trace=”none”, cexRow=0.5, cexCol=0.7,margin=c(7,30), keysize=1.5); […] 如果有人告诉你用显微镜实时观测单分子DNA聚合酶复制DNA，并用它来测序，你一定会 认为他异想天开，没有一点生物的sense。 我最初就是这样认为的，然而它不仅可以实现，而且已经实现了！这个就是被称为第三 代的测序技术，Pacific Biosciences公司推出的“Single Molecule Real Time (SMRT ™) DNA Sequencing”（单分子实时DNA测序）。 我有幸在NIH听到了这个技术发明人Stephen Turner博士的讲座，根据自己粗浅的理解 记录整理一下。 要实现单分子实时测序，有三个关键的技术。 第一个是荧光标记的脱氧核苷酸。显微镜现在再厉害，也不可能真的实时看到“单分子 ”。但是它可以实时记录荧光的强度变化。当荧光标记的脱氧核苷酸被掺入DNA链的时 候，它的荧光就同时能在DNA链上探测到。当它与DNA链形成化学键的时候，它的荧光基 团就被DNA聚合酶切除，荧光消失。这种荧光标记的脱氧核苷酸不会影响DNA聚合酶的活 性，并且在荧光被切除之后，合成的DNA链和天然的DNA链完全一样。 第二个是纳米微孔。因为在显微镜实时记录DNA链上的荧光的时候，DNA链周围的众多的 荧光标记的脱氧核苷酸形成了非常强大的荧光背景。这种强大的荧光背景使单分子的荧 光探测成为不可能。Pacific Biosciences公司发明了一种直径只有几十纳米的纳米孔[ zero-mode waveguides (ZMWs)]，单分子的DNA聚合酶被固定在这个孔内。在这么小的 孔内，DNA链周围的荧光标记的脱氧核苷酸有限，而且由于A，T，C，G这四种荧光标记 的脱氧核苷酸非常快速地从外面进入到孔内又出去，它们形成了非常稳定的背景荧光信 号。而当某一种荧光标记的脱氧核苷酸被掺入到DNA链时，这种特定颜色的荧光会持续 一小段时间，直到新的化学键形成，荧光基团被DNA聚合酶切除为止（见图）。 第三个是共聚焦显微镜实时地快速地对集成在板上的无数的纳米小孔同时进行记录。由 于我对显微原理的物理知识匮乏，而Pacific Biosciences公司又没有非常强调在这方 面的发明，不做进一步探讨。 他们还对这一技术进行进一步的优化。 第一个是把双链DNA环化反复测序。人们可以在双链DNA的两头连上发夹结构的DNA adaptor，从而使DNA环化。而DNA聚合酶就能够以环化的DNA作为模板滚环复制，反复测 一段DNA序列。这种反复测序，纠正了偶尔出现的复制错误，从而使测序精度非常高。 第二个是激发光中断测序法。DNA聚合酶虽然很稳定，但是在强大的激发光作用下酶也 是有一定寿命的。如果把激发光中断一段时间，在这段时间内DNA聚合酶继续复制DNA， 当激发光重新开启以后，人们就可以测到长DNA链后面的序列。 第三代测序技术非常可怕。1、它实现了DNA聚合酶内在自身的反应速度，一秒可以测10 个碱基，测序速度是化学法测序的2万倍。2、它实现了DNA聚合酶内在自身的 processivity（延续性，也就是DNA聚合酶一次可以合成很长的片段），一个反应就可 以测非常长的序列。 二代测序现在可以测到上百个碱基，但是三代测序现在就可以测 几千个碱基。这为基因组的重复序列的拼接提供了非常好的条件。3、它的精度非常高 ，达到99.9999%。 此外，它还有两个应用是二代测序所不具备的。 第一个是直接测RNA的序列。既然DNA聚合酶能够实时观测，那么以RNA为模板复制DNA的 逆转录酶也同样可以。RNA的直接测序，将大大降低体外逆转录产生的系统误差。 第二个是直接测甲基化的DNA序列。实际上DNA聚合酶复制A、T、C、G的速度是不一样的 。正常的C或者甲基化的C为模板，DNA聚合酶停顿的时间不同。根据这个不同的时间， […] Glancing at the code for heatmap.2 I’m fairly sure that the default is to use dist, and it’s default is in turn to use euclidean distances. The reason your attempt at passing distfun = dist(method = ‘euclidean’) didn’t work is thatdistfun (and hclustfun) are supposed to simply be name of functions. So if you want […] ## load the libraries required library(made4) require(graphics) library(cluster) library(stats) #### read data in from a comma-delimited file with header array_data<-read.csv(“Overall_Outsidetreat_sort_top100.csv”,header=T, row.names=1, check.names=T) ## array_data is a data frame type ## rows are taxa_strings and columns are samples # now transpose rows and columns # array_data<-t(array_data) ## now rows are chips and columns are taxa dim(array_data) […] 我们在分析了差异表达数据之后，经常要生成一种直观图－－热图(heatmap)。这一节就以基因芯片数据为例，示例生成高品质的热图。 比如 钢蓝渐白配色的热图 首先还是从最简单的heatmap开始。 > library(ggplot2) > library(ALL) #可以使用biocLite(“ALL”)安装该数据包 > data(“ALL”) > library(limma) > eset<-ALL[,ALL$mol.biol %in% c(“BCR/ABL”,”ALL1/AF4″)] > f<-factor(as.character(eset$mol.biol)) > design<-model.matrix(~f) > fit<-eBayes(lmFit(eset,design)) #对基因芯片数据进行分析，得到差异表达的数据 > selected <- p.adjust(fit$p.value[, 2]) <0.001 > esetSel <- eset[selected,] #选择其中一部分绘制热图 > dim(esetSel) #从这尺度上看，数目并不多，但也不少。如果基因数过多，可以分两次做图。 Features Samples 84 47 > library(hgu95av2.db) > data<-exprs(esetSel) > probes<-rownames(data) > symbol<-mget(probes,hgu95av2SYMBOL,ifnotfound=NA) > symbol<-do.call(rbind,symbol) > symbol[is.na(symbol[,1]),1]<-rownames(symbol)[is.na(symbol[,1])] > […] R当中的坐标中断一般都使用plotrix库中的axis.break(), gap.plot(), gap.barplot(), gap.boxplot()等几个函数来实现，例： gap plot > library(plotrix) > opar<-par(mfrow=c(3,2)) > plot(sample(5:7,20,replace=T),main=”Axis break test”,ylim=c(2,8)) > axis.break(axis=2,breakpos=2.5,style=”gap”) > axis.break(axis=2,breakpos=3.5,style=”slash”) > axis.break(axis=2,breakpos=4.5,style=”zigzag”) > twogrp<-c(rnorm(5)+4,rnorm(5)+20,rnorm(5)+5,rnorm(5)+22) > gap.plot(twogrp,gap=c(8,16,25,35), + xlab=”X values”,ylab=”Y values”,xlim=c(1,30),ylim=c(0,25), + main=”Test two gap plot with the lot”,xtics=seq(0,30,by=5), + ytics=c(4,6,18,20,22,38,40,42), + lty=c(rep(1,10),rep(2,10)), + pch=c(rep(2,10),rep(3,10)), + col=c(rep(2,10),rep(3,10)), + type=”b”) > gap.plot(21:30,rnorm(10)+40,gap=c(8,16,25,35),add=TRUE, + lty=rep(3,10),col=rep(4,10),type=”l”) > gap.barplot(twogrp,gap=c(8,16),xlab=”Index”,ytics=c(3,6,17,20), + ylab=”Group values”,main=”Barplot […] |
|
||||
|
Copyright © 2025 小生这厢有礼了(BioFaceBook Personal Blog) - All Rights Reserved Powered by WordPress & Atahualpa |
||||||
Recent Comments