Month: December 2017

  • 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…

  • 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…

  • 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))…

  • 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的基因数量太少会不会影响结果呢?…

  • 多个相同行列文件的合并

    #! /usr/bin/python import glob out_file1 = open(r”E:\\Git-2.11.1-64-bit\\Git\\localhost\\lianshou\\4WGCNA\\merge.txt”, “w”) MyDiction = {} def Readfile(File): file_obj = open(File) try: while True: line = file_obj.readline().strip(“\n”) if not line: break array = line.split(“\t”) if array[0] in MyDiction: MyDiction[array[0]].append(array[1]) else: MyDiction[array[0]] = [array[1]] finally: file_obj.close() return MyDiction def main(): list_dirs = glob.glob(“./*.txt”) for i in list_dirs: Readfile(i) for gene_name in…

  • pheatmap包绘制热图

    目前,绘制热图的软件有很多,我也尝试了几款软件,觉得R的pheatmap包最好用。 pheatmap安装: source(“http://bioconductor.org/biocLite.R”) biocLite(“pheatmap”) pheatmap绘制热图只需三步,导入数据结构如下图: 软件调用: library(pheatmap) 数据导入: data<- read.table(“456.txt”,head = T,row.names=1,sep=”\t”) 绘图: pheatmap(data,cluster_rows=T,cluster_cols=T,scale=”none”,border_color=”white”,color=colorRampPalette(rev(c(“red”,”white”,”blue”)))(102),file=”test.pdf”,width=3,height=5) pheatmap可调的参数有很多,可根据需要自行添加,调整