叔本华:作为意志和表象的世界(一)

第一篇  世界作为表象初论

 服从充分根据律的表象

经验和科学的客体

跳出童年时代吧,朋友,觉醒呵!
——J·J·卢梭

“世界是我的表象”:这是一个真理,是对于任何一个生活着和认识着的生物都有效的真理;不过只有人能够将它纳入反省的,抽象的意识罢了。并且,要是人真的这样做了,那么,在他那儿就出现了哲学的思考。于是,他就会清楚而确切地明白,他不认识什么太阳,什么地球,而永远只是眼睛,是眼睛看见太阳;永远只是手,是手感触着地球就会明白围绕着他的这世界只是作为表象而存在着的,也就是说这世界的存在完全只是就它对一个其他事物的,一个进行“表象者”的关系来说的。这个进行“表象者”就是人自己。如果有一真理可以先验地说将出来,那就是这一真理了,因为这真理就是一切可能的、可想得到的经验所同具的那一形式的陈述。它比一切,比时间、空间、因果性等更为普遍,因为所有这些都要以这一真理为前提。我们既已把这些形式都认作根据律的一些特殊构成形态,如果其中每一形式只是对一特殊类型的表象有效,那么,与此相反,客体和主体的分立则是所有那些类型的共同形式。客体主体分立是这样一个形式:任何一个表象,不论是哪一种,抽象的或直观的,纯粹的或经验的,都只有在这一共同形式下,根本才有可能,才可想象。因此,再没有一个比这更确切,更不依赖其他真理,更不需要一个证明的真理了;即是说:对于“认识”而存在着的一切,也就是全世界,都只是同主体相关联着的客体,直观者的直观;一句话,都只是表象。当然,这里所说的对于现在,也对于任何过去,任何将来,对于最远的和近的都有效;因为这里所说的对于时间和空间本身就有效;而又只有在时间、空间中,所有这些[过去、现在、未来、远和近]才能区别出来。一切一切,凡已属于和能属于这世界的一切,都无可避免地带有以主体为条件[的性质],并且也仅仅只是为主体而存在。世界即是表象。

这个真理决不新颖。它已包含在笛卡儿所从出发的怀疑论观点中。不过贝克莱是断然把它说出来的第一人;尽管他那哲学的其余部分站不住脚,在这一点上,他却为哲学作出了不朽的贡献。康德首先一个缺点就是对这一命题的忽略,这在本书附录中将有详尽的交代。与此相反,吠檀多哲学被认为是毗耶舍的作品,这里所谈的基本原理在那里就已作为根本命题出现了,因此印度智者们很早就认识这一真理了。威廉·琼斯在他最近《论亚洲哲学》(《亚洲研究》,第四卷第164页)一文中为此作了证,他说:“吠檀多学派的基本教义不在于否认物质的存在,不在否认它的坚实性、不可入性、广延的形状(否认这些,将意味着疯狂),而是在于纠正世俗对于物质的观念,在于主张物质没有独立于心的知觉以外的本质,主张存在和可知觉性是可以互相换用的术语。”这些话已充分地表出了经验的实在性和先验的观念性两者的共存。

在这第一篇里,我们只从上述的这一方面,即仅仅是作为表象的一面来考察这世界。至于这一考察,虽无损于其为真理,究竟是片面的,从而也是由于某种任意的抽象作用引出来的,它宣告了每一个人内心的矛盾,他带着这一矛盾去假定这世界只是他的表象,另一方面他又再也不能摆脱这一假定。不过这一考察的片面性就会从下一篇得到补充,由另一真理得到补充。这一真理,可不如我们这里所从出发的那一个,是那么直接明确的,而是只有通过更深入的探讨,更艰难的抽象和“别异综同”的功夫才能达到的。它必然是很严肃的,对于每一个人纵不是可怕的,也必然是要加以郑重考虑的。这另一真理就是每人,他自己也能说并且必须说的:“世界是我的意志。”

在作这个补充之前,也就是在这第一篇里,我们必须坚定不移地考察世界的这一面,即我们所从出发的一面,“可知性”的一面:因此,也必须毫无抵触心情地将当前现成的客体,甚至自己的身体(我们就要进一步谈到这点)都仅仅作为表象看,并且也仅仅称之为表象。我们希望往后每一个人都会确切明白我们在这样做的时候,只仅仅是撇开了意志;而意志就是单独构成世界另外那一面的东西;因为这世界的一面自始至终是表象,正如另一面自始至终是意志。至于说有一种实在,并不是这两者中的任何一个方面,而是一个自在的客体(康德的“自在之物”可惜也不知不觉的蜕化为这样的客体),那是梦呓中的怪物;而承认这种怪物就会是哲学里引人误入迷途的鬼火。

转载:叔本华哲学智慧

GSEA命令行

前边已经说到MSigDB数据库的entrez ID和symble ID都是人源的,其他物种要想做GSEA的话,就必须要ID转换,做成属于自己研究物种的数据库,这一步是关键。想要根据简单的Nr注释或其他注释来转换是行不通的,比如我研究的物种为硬骨鱼类,Nr注释结果分值最高的前20个一般也都是鱼类,就算注释到了人源的基因,也可能存在版本不对的问题。因此我们要知道MSigDB数据库基因ID的版本,然后下载对应的全基因组蛋白序列,然后比对卡阈值,做成自己的数据库。

GSEA的运行可以界面操作和命令行操作,界面操作教程很多就不多做赘述,它的不足之处除了操作繁琐之外,里面几个较大的库需要的内存较高,像我的4G渣渣电脑,好几个都运行失败,不得不转到服务器上运行;而且界面也不能批量运行。

java -cp gsea-3.0.jar -Xmx10000m xtools.gsea.Gsea -res ma/normalized_counts_all.gct -cls ma/samples.cls#COS1_versus_COS0 -gmx ma/ma.h.all.v6.1.symbols.gmt -collapse false -mode Max_probe -norm meandiv -nperm 1000 -permute gene_set -rnd_type no_balance -scoring_scheme weighted -rpt_label COS1vsCOS0_h -metric Ratio_of_Classes -sort real -order descending -create_gcts false -create_svgs false -include_only_symbols true -make_sets true -median false -num 100 -plot_top_x 20 -rnd_seed timestamp -save_rnd_lists false -set_max 2000 -set_min 5 -zip_report false -out multiple -gui false

-Xmx10000m:设置最大内存10000兆,之前界面操作失败就因为这个参数
-res:表达量文件
-cls:样本信息
-gmx:数据库文件
-nperm:Number of permutations
-rpt_label:输出文件的前缀
-metric:排序的方法,如果有重复,可以考虑使用T-test;无重复,可以考虑使用Ratio of calsses(差异倍数)或Diff of classes(差异绝对值)
-plot_top_x:默认20,代表富集分析排名最高的20个通路
-set_max:通路基因的最大数量,默认500,但由于某些通路基因数大于500,建议提高阈值
-set_min:通路基因的最小数量,默认15
-out:输出的文件夹
其他参数均为默认

版权声明:本文为博主原创文章,未经博主允许不得转载。

GSEA分析流程

什么是GSEA?

基因集富集分析 (Gene Set Enrichment Analysis, GSEA)用来评估一个预先定义的基因集的基因在与表型相关度排序的基因列表中的分布趋势,从而判断其对表型的贡献。基因集合富集分析检测基因集合而不是单个基因的表达变化,因此可以包含这些细微的表达变化,预期得到更为理想的结果。

GSEA原理

给定一个排序的基因表L和一个预先定义的基因集S (比如编码某个代谢通路的产物的基因, 基因组上物理位置相近的基因,或同一GO注释下的基因),GSEA的目的是判断S里面的成员s在L里面是随机分布还是主要聚集在L的顶部或底部。这些基因排序的依据是其在不同表型状态下的表达差异,若研究的基因集S的成员显著聚集在L的顶部或底部,则说明此基因集成员对表型的差异有贡献,也是我们关注的基因集。

GSEA分析中几个关键的概念

计算富集得分 (ES, enrichment score),ES反应基因集成员s在排序列表L的两端富集的程度。计算方式是,从基因集L的第一个基因开始,计算一个累计统计值。当遇到一个落在s里面的基因,则增加统计值。遇到一个不在s里面的基因,则降低统计值。每一步统计值增加或减少的幅度与基因的表达变化程度(更严格的是与基因和表型的关联度)是相关的。富集得分ES最后定义为最大的峰值。正值ES表示基因集在列表的顶部富集,负值ES表示基因集在列表的底部富集。

评估富集得分(ES)的显著性,通过基于表型而不改变基因之间关系的排列检验 (permutation test)计算观察到的富集得分(ES)出现的可能性。若样品量少,也可基于基因集做排列检验 (permutation test),计算p-value。

多重假设检验矫正,首先对每个基因子集s计算得到的ES根据基因集的大小进行标准化得到Normalized Enrichment Score (NES)。随后针对NES计算假阳性率。(计算NES也有另外一种方法,是计算出的ES除以排列检验得到的所有ES的平均值)

Leading-edge subset,对富集得分贡献最大的基因成员。

MSigDB数据库

MSigDB(Molecular Signatures Database)数据库中定义了已知的基因集合:http://software.broadinstitute.org/gsea/msigdb

包括H和C1-C7八个系列(Collection),每个系列内容为:H: hallmark gene sets(效应)特征基因集合,共50组(比如细胞凋亡特征基因集);C1: positional gene sets 位置基因集合,根据染色体位置,共326个;C2: curated gene sets:(专家)共识基因集合,基于通路、文献等(这部分包括我们熟悉的KEGG信号通路);C3: motif gene sets:模式基因集合,主要包括microRNA和转录因子靶基因两部分;C4: computational gene sets:计算基因集合,通过挖掘癌症相关芯片数据定义的基因集合;C5: GO gene sets:Gene Ontology 基因本体论,包括BP(生物学过程biological process,细胞原件cellular component和分子功能molecular function三部分);C6: oncogenic signatures:癌症特征基因集合,大部分来源于NCBI GEO 未发表芯片数据;C7: immunologic signatures: 免疫相关基因集合。

软件下载:http://software.broadinstitute.org/gsea/downloads.jsp

数据格式:http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats

一:准备数据(GSEA 排序需要三个文件,如下)

1.表达量文件(*.gct)

2.说明文件(*.cls)

3. 基因集合(.gmt)

or

MSigDB提供的数据库是人源的,所以你的研究物种不是人的话,就需要自己制作gmt文件,这也是GSEA分析的难点和关键所在

DESeq分析流程

source("https://bioconductor.org/biocLite.R")
biocLite("DESeq")
library("DESeq")
database <- read.table(file = "macrophage_genes.count_table.matrix", sep = "\t", header = T, row.names = 1)
countData <- database[,1:3]
condition <- factor(c("A","B","C"))#无生物学重复
#type <- factor(c(rep("A",3), rep("B",3))) 有生物学重复
database <- round(as.matrix(countData))#取整数型
dds <- newCountDataSet(database,condition)
dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds, method="blind", sharingMode="fit-only" )#无生物学重复
dds <- estimateDispersions(dds) # 有生物学重复
resAvsB <- nbinomTest(dds,"A","B")
resAvsC <- nbinomTest(dds,"A","C")
table(resAvsB$pval <0.05)
table(resAvsC$pval <0.05)
resAvsB <- resAvsB[order(resAvsB$pval),]#排序
resAvsC <- resAvsC[order(resAvsC$pval),] 
plotMA(resAvsB, ylim = c(-5,5), col = ifelse(resAvsB$pval>=0.05, "gray32", "red3"),linecol = "#ff000080" )#画MA图
norcounts <- counts(dds, normalized=T)#提取标准化后的counts
write.table(as.data.frame(norcounts),
file="normalized_counts_ABC.txt",
sep="\t",
quote = F)
resAvsBup <- subset(resAvsB, pval < 0.05 & log2FoldChange >0 )
resAvsBdown <- subset(resAvsC, pval < 0.05 & log2FoldChange <0 )

Python replace()方法

描述

Python replace() 方法把字符串中的 old(旧字符串) 替换成 new(新字符串),如果指定第三个参数max,则替换不超过 max 次。

语法

replace()方法语法:

str.replace(old, new[, max])

参数

  • old — 将被替换的子字符串。
  • new — 新字符串,用于替换old子字符串。
  • max — 可选字符串, 替换不超过 max 次

返回值

返回字符串中的 old(旧字符串) 替换成 new(新字符串)后生成的新字符串,如果指定第三个参数max,则替换不超过 max 次。

实例

以下实例展示了replace()函数的使用方法:

#!/usr/bin/python

str = "this is string example....wow!!! this is really string";
print str.replace("is", "was");
print str.replace("is", "was", 3);

以上实例输出结果如下:

thwas was string example....wow!!! thwas was really string
thwas was string example....wow!!! thwas is really string

转载:http://www.runoob.com/python/att-string-replace.html

ggplot2:气泡图及散点图小结

气泡图也属于散点图的一种,在散点图的基础上改变点的形状,大小和颜色

1.如何改变点的形状

用自带的mtcars演示

p = ggplot(mtcars,aes(wt,mpg))
p + geom_point(shape=17,size=4)
p + geom_point(aes(shape = factor(cyl)),size=4)

2.改变点的大小:气泡图

加载包、数据的读入及数据格式

library(ggplot2)
pathway = read.table("R0-vs-R3.path.richFactor.head20.tsv",header=T,sep="\t")

3.初始化数据

pp = ggplot(pathway,aes(richFactor,Pathway))
pp + geom_point()

气泡图(三维数据),size 可以是一个值,也可以是数据中的一列

pp + geom_point(aes(size=R0vsR3))

添加颜色变化——四维数据

scale_colour_gradient:自定义一个连续型的配色,常用于热图;

pbubble = pp + geom_point(aes(size=R0vsR3,color=-1*log10(Qvalue)))
pbubble
pbubble + scale_colour_gradient(low="green",high="red")


expression函数改变样式,[]是用来添加下标,^是用来添加上标

pr = pbubble + scale_colour_gradient(low="green",high="red") + labs(color=expression(-log[10](Qvalue)),size="Gene number",x="Rich factor",y="Pathway name",title="Top20 of pathway enrichment")
pr
pr + theme_bw()

散点图小结

• 颜色: color/colour (可渐变)
• 形状:shape
• 大小:size (可渐变)
• 透明度:alpha (练习)
• 改变坐标轴的范围: xlim()和ylim()函数
• 改变标题,坐标轴标题和legend标题:labs()函数
• 保存图片:ggsave()函数
支持多种图片格式:png,pdf,svg等

ggsave("pathway_enrichment.png")
ggsave("pathway_enrichment_gray.png",pr)
ggsave("volcano4.pdf",volcano,width=4,height=4)
ggsave("volcano8.pdf",volcano,width=8,height=8)

绘制直线小结

• geom_abline(slope=xx,intercept=xx): 绘
制直线
• geom_hline(yintercept=xx):绘制水平线
• geom_vline(xintercept=xx):绘制垂直线
• size: 线条粗细
• color: 线条颜色
• linetype: 线条类型

详细讲解参考ggplot2说明文档基迪奥在线视频

ggplot2:火山图

火山图是散点图的一种,也是描点,用 geom_point()绘图

1.加载包,读入数据及数据格式

library(ggplot2)
data = read.table("R0-vsR3.isoforms.filter.tsv",header=T,row.names=1)

2.绘图

r03 = ggplot(data,aes(log2FC,-1*log10(FDR)))
r03 + geom_point()

3.改变点的颜色

r03 + geom_point(color="red")
r03 + geom_point(aes(color="red"))
r03 + geom_point(aes(color=significant)) # 按照“ significant”这一列定义点的颜色;

4.设置坐标轴范围和标题

函数xlim(),ylim(),规定X、 Y轴范围;labs(title=“..”,x=“..”,y=“..”),确定标题和坐标轴标签;expression函数改变样式,[]是用来添加下标,^是用来添加上标

r03xy = r03 + geom_point(aes(color=significant)) + xlim(-4,4) + ylim(0,30)
r03xy + labs(title="Volcano plot",x="log2(FC)")
r03xy + labs(title="Volcano plot",x=expression(log[2](FC)),y=expression(-log[10](FDR)))

5.自定义颜色

scale_color_manual(): 自定义颜色配色;

r03xyp = r03xy + labs(title="Volcano plot",x=expression(log[2](FC)),y=expression(-log[10](FDR)))
r03xyp + scale_color_manual(values = c("green","black", "red"))
volcano = r03xyp + scale_color_manual(values = c("#00ba38","#619cff", "#f8766d"))

6.添加阈值线

geom_hline()添加水平线; geom_vline()添加垂直线

volcano+geom_hline(yintercept=1.3)+geom_vline(xintercept=c(-1,1))

7.改变线条类型

linetype参数0 = blank, 1 = solid, 2 = dashed, 3 = dotted, 4 = dotdash, 5 = longdash, 6 = twodash

volcano+geom_hline(yintercept=1.3,linetype=4)+geom_vline(xintercept=c(-1,1),linetype=4)

详细讲解查看帮助文档基迪奥在线视频

ggplot2:散点图

1.加载ggplot2包,读取数据及数据格式

library(ggplot2)
fpkm = read.table("all.fpkm",header=T,row.names=1)
head(fpkm,10)

2.绘图

mp = ggplot(fpkm,aes(C1_FPKM,C2_FPKM))
mp + geom_point()

3.取log10后重新作图

mp = ggplot(fpkm,aes(log10(C1_FPKM),log10(C2_FPKM)))
mp + geom_point()

4.添加直线

slope设置斜率,intercept设置截距,color设置线条颜色,size设置线条粗细

mp + geom_abline(slope=1,intercept=0,color="red",size=2) + geom_point()

ggplot2的核心思想就是图层,可以像PS一样一层层的添加图层,以加号连接,更详细参数可参考ggplot2的帮助文档http://docs.ggplot2.org/current/

视频教程:http://www.omicshare.com/class/home/index/singlev?id=34

WGCNA分析流程

WGCNA是一个R包,对一个完全不会R的人来说,确实费了一番功夫,不过也将我对R的学习提前提上日程。

分析步骤:

1.WGCNA安装

source("http://bioconductor.org/biocLite.R")
biocLite(c("AnnotationDbi", "impute", "GO.db", "preprocessCore"))
install.packages("WGCNA")

2.输入数据准备

原始数据https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE48213,共56个样本;对数据的合并可参考http://www.biotrainee.com:8080/thread-603-1-1.html,输入文件格式如下:

 

setwd('WGCNA/')
fpkm=read.table('test.txt',sep = '\t',stringsAsFactors = F)
head(fpkm)
dim(fpkm) ##共56个细胞系,36953个基因的表达矩阵
names(fpkm)##细胞系名称

3.数据预处理

library(reshape2)
library('WGCNA')
RNAseq_voom <- fpkm
WGCNA_matrix = t(RNAseq_voom[order(apply(RNAseq_voom,1,mad), decreasing = T)[1:3000],]) ##top 3000 mad genes,t参数将矩阵变为符合WGCNA要求的形式:行名为gene,列名为样品
datExpr <- WGCNA_matrix
save(datExpr, file = "AS-green-FPKM-01-dataInput.RData")

4.确定最佳beta值

 

powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# Plot the results:
##sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")

最佳的beta值就是sft$powerEstimate,此例为6

5.一步法构建共表达矩阵

net = blockwiseModules(datExpr, power = 6, maxBlockSize = 6000,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "AS-green-FPKM-TOM",
                       verbose = 3)
table(net$colors)
#绘图结果展示
mergedColors = labels2colors(net$colors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
#结果保存
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
table(moduleColors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree,
     file = "AS-green-FPKM-02-networkConstruction-auto.RData")

#对样本做个系统聚类树
#明确基因数和样本数
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
datExpr_tree<-hclust(dist(datExpr), method = "average")
par(mar = c(0,5,2,0))
plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2, 
     cex.axis = 1, cex.main = 1,cex.lab=1)

6.导出网络到Cytoscape

# Recalculate topological overlap if needed
TOM = TOMsimilarityFromExpr(datExpr, power = 6);
# Read in the annotation file
# annot = read.csv(file = "GeneAnnotation.csv");
# Select modules需要修改,选择需要导出的模块颜色
modules = c("turquoise", "blue");
# Select module probes选择模块探测
probes = colnames(datExpr)
inModule = is.finite(match(moduleColors, modules));
modProbes = probes[inModule];
#modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)];
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
# Export the network into edge and node list files Cytoscape can read
cyt = exportNetworkToCytoscape(modTOM,
                               edgeFile = paste("AS-green-FPKM-One-step-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
                               nodeFile = paste("AS-green-FPKM-One-step-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
                               weighted = TRUE,
                               threshold = 0.02,
                               nodeNames = modProbes,
                               #altNodeNames = modGenes,
                               nodeAttr = moduleColors[inModule]);

7.分析网络可视化

用heatmap可视化权重网络,heatmap每一行或列对应一个基因,颜色越深表示有较高的邻近

options(stringsAsFactors = FALSE);
lnames = load(file = "AS-green-FPKM-01-dataInput.RData");
lnames
lnames = load(file = "AS-green-FPKM-02-networkConstruction-auto.RData");
lnames
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
#可视化全部基因网络
# Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# calculated during module detection, but let us do it again here.
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM = dissTOM^7;
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA;
# Call the plot function
#sizeGrWindow(9,9)
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
#随便选取1000个基因来可视化
nSelect = 1000
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There's no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
#sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

8.多步法构建网络

#2.选择合适的阀值,同上
#3. 网络构建:(1) Co-expression similarity and adjacency
softPower = 6;
adjacency = adjacency(datExpr, power = softPower);
#(2) 邻近矩阵到拓扑矩阵的转换,Turn adjacency into topological overlap
TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM
# (3) 聚类拓扑矩阵
#Call the hierarchical clustering function
geneTree = hclust(as.dist(dissTOM), method = "average");
# Plot the resulting clustering tree (dendrogram)
#sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
     labels = FALSE, hang = 0.04);
#(4) 聚类分支的休整dynamicTreeCut
# We like large modules, so we set the minimum module size relatively high:
minModuleSize = 30;
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                            deepSplit = 2, pamRespectsDendro = FALSE,
                            minClusterSize = minModuleSize);
table(dynamicMods)
#4. 绘画结果展示
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
# Plot the dendrogram and colors underneath
#sizeGrWindow(8,6)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05,
                    main = "Gene dendrogram and module colors")
#5. 聚类结果相似模块的融合,Merging of modules whose expression profiles are very similar
#在聚类树中每一leaf是一个短线,代表一个基因,
#不同分之间靠的越近表示有高的共表达基因,将共表达极其相似的modules进行融合
# Calculate eigengenes
MEList = moduleEigengenes(datExpr, colors = dynamicColors)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");
# Plot the result
#sizeGrWindow(7, 6)
plot(METree, main = "Clustering of module eigengenes",
     xlab = "", sub = "")
#选择有75%相关性的进行融合
MEDissThres = 0.25
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
# Call an automatic merging function
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors;
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEs;
#绘制融合前(Dynamic Tree Cut)和融合后(Merged dynamic)的聚类图
#sizeGrWindow(12, 9)
#pdf(file = "Plots/geneDendro-3.pdf", wi = 9, he = 6)
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
                    c("Dynamic Tree Cut", "Merged dynamic"),
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
#dev.off()
# 只是绘制融合后聚类图
plotDendroAndColors(geneTree,mergedColors,"Merged dynamic",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
#5.结果保存
# Rename to moduleColors
moduleColors = mergedColors
# Construct numerical labels corresponding to the colors
colorOrder = c("grey", standardColors(50));
moduleLabels = match(moduleColors, colorOrder)-1;
MEs = mergedMEs;
# Save module colors and labels for use in subsequent parts
save(MEs, moduleLabels, moduleColors, geneTree, file = "AS-green-FPKM-02-networkConstruction-stepByStep.RData")
#6. 导出网络到Cytoscape
# Recalculate topological overlap if needed
TOM = TOMsimilarityFromExpr(datExpr, power = 6);
# Read in the annotation file
# annot = read.csv(file = "GeneAnnotation.csv");
# Select modules需要修改
modules = c("brown", "red");
# Select module probes
probes = colnames(datExpr)
inModule = is.finite(match(moduleColors, modules));
modProbes = probes[inModule];
#modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)];
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)
# Export the network into edge and node list files Cytoscape can read
cyt = exportNetworkToCytoscape(modTOM,
                               edgeFile = paste("AS-green-FPKM-Step-by-step-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
                               nodeFile = paste("AS-green-FPKM-Step-by-step-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
                               weighted = TRUE,
                               threshold = 0.02,
                               nodeNames = modProbes,
                               #altNodeNames = modGenes,
                               nodeAttr = moduleColors[inModule]);
#=====================================================================================
#  分析网络可视化,用heatmap可视化权重网络,heatmap每一行或列对应一个基因,颜色越深表示有较高的邻近
#=====================================================================================
options(stringsAsFactors = FALSE);
lnames = load(file = "AS-green-FPKM-01-dataInput.RData");
lnames
lnames = load(file = "AS-green-FPKM-02-networkConstruction-stepByStep.RData");
lnames
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
#1. 可视化全部基因网络
# Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# calculated during module detection, but let us do it again here.
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM = dissTOM^7;
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA;
# Call the plot function
#sizeGrWindow(9,9)
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
#随便选取1000个基因来可视化
nSelect = 1000
# For reproducibility, we set the random seed
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# There's no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
# Open a graphical window
#sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")


根据基因间表达量进行模块聚类(前边得到)和所得到的各模块间的相关性

MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
MET = orderMEs(MEs)
sizeGrWindow(7, 6) 
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)


脚本来自:
http://tiramisutes.github.io/2016/09/14/WGCNA.html

https://bioconductor.org/packages/devel/bioc/vignettes/CVE/inst/doc/WGCNA_from_TCGA_RNAseq.html

http://www.bio-info-trainee.com/2535.html

WGCNA原理及应用

WGCNA介绍:

WGCNA(weighted gene co-expression network analysis,权重基因共表达网络分析)是一种分析多个样本基因表达模式的分析方法,可将表达模式相似的基因进行聚类,并分析模块与特定性状或表型之间的关联关系,因此在疾病以及其他性状与基因关联分析等方面的研究中被广泛应用。

WGCNA算法是构建基因共表达网络的常用算法(详解:http://www.jianshu.com/p/94b11358b3f3)。WGCNA算法首先假定基因网络服从无尺度分布,并定义基因共表达相关矩阵、基因网络形成的邻接函数,然后计算不同节点的相异系数,并据此构建分层聚类树(hierarchical clustering tree),该聚类树的不同分支代表不同的基因模块(module),模块内基因共表达程度高,而分属不同模块的基因共表达程度低。最后,探索模块与特定表型或疾病的关联关系,最终达到鉴定疾病治疗的靶点基因、基因网络的目的。在该方法中module被定义为一组具有类似表达谱的基因,如果某些基因在一个生理过程或不同组织中总是具有相类似的表达变化,那么我们有理由认为这些基因在功能上是相关的,可以把他们定义为一个模块(module)。这似乎有点类似于进行聚类分析所得到结果,但不同的是,WGCNA的聚类准则具有生物学意义,而非常规的聚类方法(如利用数据间的几何距离),因此该方法所得出的结果具有更高的可信度。当基因module被定义出来后,我们可以利用这些结果做很多进一步的工作,如关联性状,代谢通路建模,建立基因互作网络等。

 

WGCNA的用处:

 

  • 这类处于调控网络中心的基因称为核心基因(hub gene),这类基因通常是转录因子等关键的调控因子,是值得我们优先深入分析和挖掘的对象。
  • 在网络中,被调控线连接的基因,其表达模式是相似的。那么它们潜在有相似的功能。所以,在这个网络中,如果线条一端的基因功能是已知的,那么就可以预测线条另一端的功能未知的基因也有相似的功能。

 

下面的问答来自基迪奥,也能加深对WGCNA的理解

问1、调控网络和共表达网络有什么区别?
答:调控网络是个更广泛的概念,而共表达网络是调控网络的一种。
理论上我们可以利用各类信息构建调控网络(表达相关性,序列靶向关系、蛋白互作关系),另外调控网络构建的信息既可以来源真实的实验验证的关系,也可以来源生物信息的预测。而共表达网络特指利用基因间的表达相关性预测基因间调控关系的方法,而WGCNA又是共表达网络分析中最有效的方法之一。

问2、WGCNA分析适合的生物物种范围有规定么?
答:没有限制。对于任何物种中心法则都是存在的,调控关系对于任何物种都是存在的,所以WGCNA没有物种限定。

问3、同一物种,不同来源的转录组数据(比如不同文章/资料来源的),可以放在一起做WGCNA分析吗?
答:只要样本间有相似的生物学意义,是可以合并在一起做分析的。但要注意,不同批次之间的样本是有批次效应的,所以可能会带来一些误差,但是是可以放在一起分析的。

问4、相同材料不同处理之间,可以放在一起做WGCNA分析吗?比如重金属和盐碱处理。
答:可以的。这也正式WGCNA强大的地方,其可以将不同处理的样本,合并在一起做分析。其他方法则不一定有这么强大的能力,比如做基因表达趋势分析时,如果样本涉及到多个处理不同时期的时候,就不好合并分析(或合并后难以解读)。但WGCNA的方法关注的是调控关系,所以不管是多少个处理组,都可以很好的整合在一起做分析。

问5、不同批次的数据能放一起做WGCNA吗?
答:可以的。虽然有批次的干扰,但是干扰对WGCNA网络没有太大影响。因为WGCNA不是做差异分析,而是基因的共表达。因为批次效应理论上不影响相关性。

问6、不同类型的材料,比如亲本和F1,适合放一起进行WGCNA么?
答:如果是一个作图群体,当然亲本与F1是可以放在一起分析的,因为你只关心基因的表达模式,所以把亲本加进来是没有问题的。

问7、没有生物学重复,共3组,每组5个时间点能够做吗?
答:理论上有15个样本,是可以做WGCNA分析的。并且,分析出来的结果对你的研究应该是非常有用的。至少他会比趋势分析更有意义,更加准确。

问8、一般说WGCNA的样品不少于15个,15个样品考虑重复吗?不同倍性的材料呢?
答:15个样本这个是包含了生物学重复,比如5个时间点3个重复;在RNA-seq里面建议不要用不同倍性材料加进来。除非是有参考的多倍体,如果是无参的多倍体,不同倍性之间差异太大,会让调控网络不准确。所以用单一倍性的材料做调控网络会更加准确。

问9、可以将RNA-seq数据与蛋白组数据,甲基化数据放一起做WGCNA分析?
答:不能与蛋白数据一起分析。因为WGCNA是基于相关系数的算法。所以最好一起分析的数据变异度是类似的,RNAseq变异非常大,而蛋白的数据变异很小,两者的变化不在一个数量级上面。所以两种数据放在一起分析不合理。
但RNA数据可以尝试跟甲基化数据一起分析。当然我们也建议RNA数据与代谢组数据一起分析,因为代谢组的数据变异也非常大。

问10、表达量和表达的基因数目差异太大的样品可以一起分析吗?比如样品A有2k个gene表达 而样品B有2w个gene表达了 AB可以一起分析吗?
答:做WGCNA分析的时候,不能脱离生物学意义,既然要分析调控网络,那么应该分析有相似生物学意义的一组基因,比如说拿相似组织来一起做分析,比如不应该拿大脑的样本与脚趾的样本合并在一起做分析,因为很显然,这两个组织没有关联。如果两个样本之间是有相关联的生物学意义,哪怕表达的基因数不一样,或表达模式差异很大,那依然可以放在一起分析;但如果样本之间完全没有生物学意义,那么分析就没有意义。

问11、实验设计是case3个时间点(各点都有三个重复),control同样的3个时间点(每点三个重复),WGCNA怎么做?3个时间点和case-control两个因素能同时考虑进来分析吗?
答:可以的。做WGCNA是更加合理的,因为有两个梯度的样本,如果只是做差异分析的话,逻辑可能非常复杂,做WGCNA分析是对样本特性更好的解析,可以直观看到基因在六个处理组里面是怎样表达的。

问12、可以拿混合样本分析吗?比如一个病原细菌跟人类细胞的基因,能说明细菌跟人类细胞基因有调控关系吗?
答:可以。前提是病原菌有足够的数据并定量准确,并且这个分析是非常有意义的,最后可以说明这些病原菌可以调控哪些宿主基因。

问13、但是病原宿主混合分析的话,宿主蛋白不能分泌到宿主体内岂不是WGCNA生物学上也没有意义吗?
答:依然有意义。即使病原的基因没有分泌到宿主里面,但是病原的蛋白是会影响宿主基因的调控的,比如某个细菌感染某个植物,虽然细菌的蛋白不能直接分泌到植物体内,但会影响植物蛋白的分泌。混在一起分析依然是有意义,可以看到植物里面到底哪个基因对细菌蛋白产生应答作用。

问14、芯片数据两分类,每组20个样本,能否每组单独做WGCNA?
答:可以。WGCNA还有一种重要功能是做两个网络的比较,比如病人20个样本做一个调控网络,健康人做一个调控网络,然后两个网络做比较。

问15、WGCNA可以用来分析lncRNA对下游基因的调控分析吗?
答:可以。WGCNA网络有利于预测lncRNA的潜在功能。

问16、构建网络是用所有表达基因还是差异基因?
答:这个是具体问题具体分析。如果使用所有的基因分析,会导致运算量非常大。而也不是所有的基因在这个实验中都有生物学意义,所以我们会提前做一些过滤。
但用于分析的基因不一定是差异表达基因,有时可以用差异表达基因做一个并集,或通过计算变异系数将变异系数低的基因以及低表达的基因去除。但注意,如果你有关心的特定目标基因的话,应该尽量给予保留。

问17、关注某一个pathway上的基因以及调控因子之间的相关性,构建WGCNA网络的时候属于这个pathway的基因数量太少会不会影响结果呢?
答:这不是问题。在一个调控网络里面,样本的某个pathway上,并不是所有基因参与调控(或存在差异性),所以在做WGCNA分析的时候,会做一些过滤,将有变化的基因挑出来再做分析。即分析的是某个pathway上有变化的基因,不需要分析pathway上所有的基因,只需要分析那些变化的基因就够了。

问18、前期筛选的时候,要选出在所有样本中变异系数比较大的基因呢?还是直接用差异表达的基因取并集?用基因还是转录本,哪个好呢?
答:两则都可以,我推荐使用变异系数,选择那些变异较大的基因,来做下面的分析。然后建议用基因不要用转录本,因为转录本的定量是不准确的。

问19、变异系数一般取多大?

答:具体问题具体分析。例如,没有特定目标的时候,可以先计算变异系数,将变异系数的百分之前50来做分析,把变异系数偏低的后面一半过滤掉。

问20、输入数据用FPKM合适吗?
答:可以。

问21、RNA seq数据是RSEM值怎么办?
答:RSEM值原始输出结果为reads数,如果是RSEM值建议做一个RPKM校正再做分析。

问22、除了RPKM值以外,做WGANA是否还需要其他数据?TCGA数据可否来做WGCNA分析?
答:在做WGCNA分析必须要用表达量数据,但TCGA的数据某些层级没有表达量数据,没有表达量数据自然就无法做WGCNA分析。

问23、请问输入的基因样本的矩阵的时候,要不要对数据标准化?
答:做WGCNA分析的时候,不需要对数据进行标准化,输入RPKM值就足以做这个分析。虽然一些文章会做log2处理,但我认为取了LOG2后,会让一些表达关系没有那么丰富。

问24、每个样本有3个生物学重复,不需要对三个重复的表达量求平均值代表该样本吗?
答:注意,做WGCNA的时候每个样本是独立的,三个生物学重复样本是全部导入做分析,不是取均值再做分析,每个样本都是独立的。

问25、如果3个生物学重复,做WGCNA的时候是取三个值,还是用cuffdiff处理后取一个值?
答:如果是生物学重复样本进行调控网络分析,每个样本独立使用,而不是取均值。

问26、请问将样本信息同模块特征值进行相关性分析的时候,样本信息是怎么处理的呢?比如不同取样点、不同性别什么的,这不是数量性状信息的,这种情况应该怎么处理呢?
答:样本的任何信息都可以做模块相关性分析。比如相关时间点,可以按照先后量化为12134567。又如不同性别,男与女,可以定义为1,-1。任何性状量化为数字后,都可以进行相关性分析。

问27、怎么将模块与性状对应起来呢有些性状不好量化,如果直接将模块与分组对应,如何实现, 不需要量化指标么?
答:首先需要将性状量化,如果无法将性状量化,那么就无法分析。至于分组信息,也可以量化为类似00001111000(1代表一种组别,2代表另一组组别),实现分组信息的数字化。

问28、基因数量为3w左右时,modules数量为多少结果较为理想?怎么评价聚类效果的好坏?
答:modules数量没有标准,modules数量无法评估模块分的好坏,分组是否合理应该看树的树形图,比如树的分支很清晰就说明模块式清晰的。modules数量数由生物性状决定的。比如样本表达信息很丰富的时候,modules数量会很多;如果样本的基因表达相对单一,modules数量就会比较少。

问29、我运行例子的时候,得出来基因之间的direction全是undirected,这和前面的几种关系有什么区别?
答:WGCNA是一个undirected的方法,它的网络是无方向的,有相关关系但是无方向。

问30、如果做有向网络的构建,您推荐那些方法?
答:很多方法,例如贝叶斯的方法。

问31、非模式物种可以得出基因之间的相互关系类型么?得出的结果也是undirected么?
答:WGCNA是基于表达两处理的,所以即使是非模式生物,当然也可以他们之间关系,并且关系也是一个无向网络。

问32、选择几个表型数据进行结合分析比较好
答:越多越好,看实验设计。

问33、感染小鼠,5个时间点,3个重复,找不到合适的表型怎么办?
答:如果找不到合适表型,可以找某个时间点应答的基因,本身基因的表达趋势已经有某种生物学意义的。没有找到合适表型,也可以看变化趋势。不一定要做表型的相关分析,其他分析也是很有趣的。例如,可以对模块功能的富集分析,其实都是可以帮助你找到特定模块的。所以不用纠结于做某个表型的关联分析。

问34、weight就是tom值吗?
答:是的。

问35、剪模块是怎么做的?是根据TOM划分吗?需要自己设定,还是R自动的?
答:剪模块是R中自动完成的,不需要划分,但合并的时候你可以设定一个指标,比如差异度是0.25。

问36、看WGCNA说明是用相异矩阵D(D=1-TOM)去做聚类,然后动态剪切?
答:用TOM值来构建矩阵,TOM值就是两个样本的相似度,1-TOM值就是两个样本的差异度,相似度与差异度可以理解为一个东西,并不矛盾。

问37、模块特征值和样本性状相关分析的具体方法是?
答:R包用的是计算相关系数的方法。

问38、WGCNA里面一般会提到hubgene,如何确定hubgene?
答:在WGCNA分析里面,每个基因都会计算连通性,连通性高的就是hubgene。

问39、在R中安装“”WGCNA“”说不适合R3.3.1,那适合哪个版本?
答:WGCNA应该是所有版本都适合,如果版本没有可以考虑降低R软件的版本,这个对分析没有影响。因为不同R版本是一样的。

问40、用STEM分析的时候拟合多少个模型合适?
答:建议不要超过20个。模块太多不好分析。

参考网站:

http://tiramisutes.github.io/2016/09/14/WGCNA.html

http://www.jianshu.com/p/94b11358b3f3

http://www.omicshare.com/class/home/index/classdetail?id=20