学术咖

Mindblown: a blog about philosophy.

  • 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可调的参数有很多，可根据需要自行添加，调整

  • GO功能富集

    前几天实验室一个师兄给我一个质谱结果，让帮忙做下go的功能富集，数据格式大概是这样的： 由于之前做go和kegg时都是跑流程，像这种针对性的go富集还没做过，说到底，还是由于自己手上缺少数据，没有属于自己的项目，很多细节性的问题都没有经历过。但这不妨碍咱一颗求知的心，我们都是在学习中成长。由于没事的时候逛论坛逛的比较频繁，知道数据的第二列是UniPro数据库的accession，然后该怎么办呢？作为生信人，Google是少不了的，看到Google结果，瞬间明了。根据Google的指引我从网上下载了UniProt数据库里的idmapping.tb.gz文件（wget -c -t 10000 ftp://ftp.pir.georgetown.edu/databases/idmapping/idmapping.tb.gz），大概18G左右，数据结构如下： 一共有22列，依次分别是：UniProtKB accession，UniProtKB ID，EntrezGene，RefSeq，NCBI GI number，PDB，Pfam，GO，PIRSF，IPI，UniRef100，UniRef90，UniRef50，UniParc，PIR-PSD accession，NCBI taxonomy，MIM，UniGene，Ensembl，PubMed ID，EMBL/GenBank/DDBJ，EMBL protein_id；这就有意思了，数据的第八列就是我们想要的go信息。更有意思的是，有了这个数据库信息，我们就可以根据不同数据库的注释信息做go富集啦！ 下面要做的是写一个脚本，根据师兄给的结果调出对应的go号，对于会编程的人来说，这点自然不在话下，代码如下： import sys USAGE = “\nusage: python %s idmapping.tb.gz blastout outputfile outputfile2\n” % sys.argv[0] if len(sys.argv) != 5: print USAGE sys.exit() def parseIDmapping(filename): UniProt_GO = {} with open(filename, ‘r’) as f: for line in f: lsplit = line.rstrip().split(“\t”) if lsplit[7]:…

  • sed，awk

    sed工具 sed本身也是一个管道命令，可以分析standard input的，而且sed还可以将数据进行替换、删除、新增、选取特定行等的功能。 用法：sed [-nefr] 动作 参数： -n：使用安静模式。在一般sed的用法中，所有来自stdin的数据一般都会被列出到屏幕上。但如果加上-n参数后，则只有经过sed特殊处理的那一行（或者操作）才会被列出。 -e：直接在命令行模式上进行sed的动作编辑。 -f：直接将sed的动作写在一个文件内，-f filename则可以执行filename内的sed动作。 -r：sed的动作支持的是扩展型正则表达式的语法（默认是基础正则表达式语法）。 -i：直接修改读取的文件内容，而不是屏幕输出。 动作说明：[n1[,n2]]function n1,n2:不见得会存在，一般代表选择进行动作的行数，举例来说，如果我的动作是需要在10到20行之间进行的，则‘10,20[动作行为]’ function有下面这些参数： a：新增，a的后面可以接字符串，而这些字符串会在新的一行出现（目前的下一行）； c：替换，c的后面可以接字符串，这些字符串可以替换n1，n2之间的行； d：删除，因为是删除，所以d后面通常不接任何参数； i：插入，i的后面可以接字符串，而这些字符串会在新的一行出现（目前的上一行）； p：打印，将某个选择的数据打印出来。通常p会与参数sed -n一起运行； s：替换，可以直接进行替换工作。通常这个s的动作可以搭配正则表达式。 例一：删除第2到5行 注意sed后面接的动作，务必以”两个单引号括住！ 如果体型变换一下，删除第三到最后一行，则是‘nl /etc/passwd | sed ‘3,$d’’，这个$代表最后一行。 例二：将第2~5行的内容替换成为‘No 2-5 number’ sed 另一个强大的用法，部分数据的查找并替换：sed ‘s/要被替换的字符串/新的字符串/g’ awk:好用的数据处理工具 相比于sed常常作用于一整行的处理，awk则比较倾向于将一行分成数个‘字段’来处理。因此。awk适合小型的数据处理 用法：awk ‘条件类型 1{动作 1} 条件类型 2{动作 2} …’ filename awk的处理流程： 读入第一行，并将第一行的数据填入$0,$1,$2等变量当中； 依据条件类型的限制，判断是否需要进行后面的动作； 做完所有的动作与条件类型； 若还有后续的‘行’的数据，则重复上面1-3的步骤，直到所有的数据都读完为止。 变量名字及意义： NF：每一行（$0）拥有的字段总数 NR：目前awk所处理的是‘第几行’数据…

  • 排序命令：sort，wc，uniq

    sort 用法：sort [-fbMnrtuk] [file or stdin] 参数： -f：忽略大小写 -b：忽略最前面的空格部分 -M：以月份名字来排序，如JAN，DEC等的排序方法 -n：使用‘纯数字’排序（默认是以文字类型来排序的） -r：反向排序 -u：就是uniq，相同的数据中，仅出现一行代表 -t：分隔符，默认是用tab键来分割 -k：以那个区间（field）来进行排序 wc 用法：wc [-lwm] 参数： -l：仅列出行 -w：仅列出多少字 -m：多少字符 uniq 用法：uniq [-ic] 参数： -i：忽略大小写 -c：进行计数 摘自：《鸟哥的私房菜》第三版 基础学习篇

  • 选取命令：cut，grep

    cut cut -d ‘分隔字符’ -f fields cut -c 字符范围

  • 通配符与特殊符号

    通配符是bash操作环境中一个非常有用的功能，利用它我们处理数据就更加方便。 *：代表0个到无穷多个任意字符 ?: 代表一定有一个任意字符 []: 同样代表一定有一个在中括号内的字符（非任意字符） [-]: 若有减号在中括号内时，代表『在编码顺序内的所有字符』。例如 [0-9] 代表 0 到 9 之间的所有数字，因为数字的语系编码是连续的 [^]: 若中括号内的第一个字符为指数符号 (^) ，那表示『反向选择』，例如 [^abc] 代表 一定有一个字符，只要是非 a, b, c 的其他字符就接受的意思 特殊字符 #: 注释，这个最常用在script中，视为说明.其后的数据均不执行 \: 转义符号，将特殊字符或通配符还原成一般字符 |: 分隔两个管线命令的界定 ;: 连续性命令的界定（注意，与管线命令并不相同） ~: 用户的主文件夹 $: 使用变量前导符 &: 将指令变成在背景下工作 !: 逻辑运算中的“非” /: 路径分隔符号 >,>>: 数据流重定向，输出导向，代表替换和累加 <,<<: 数据流重定向，输入导向 ”: 单引号，不具有变量置换的功能 “”: 具有变量置换的功能 “: 两个“`”中间为可以先执行的指令 (): 中间为子shell的起始与结束 {}:…

  • shell的变量功能（二）

    变量键盘读取，数组与声明：read，array，declare 1. read 读取来自键盘输入的变量，常被用在 shell script 的撰写当中。 用法： 2. declare / typeset declare 或 typeset 是一样的功能，就是在声明变量的类型。如果使用 declare 后面并没有接任何参数，那么 bash 就会主动的将所有的变量名称与内容通通叫出来， 就好像使用 set 一样。 用法：declare [-aixr] variable 参数： -a ： 将后面名为 variable 的变量定义成为数组 (array) 类型 -i ： 将后面名为 variable 的变量定义成为整数数字 (integer) 类型 -x ： 用法与 export 一样， 就是将后面的 variable 变成环境变量； -r ： 将变量设置成为 readonly 类型， 该变量不可被更改内容，…

  • 我的Python之路

    大概去年六七月份，那时还不懂什么是生物信息学，什么是编程，Python更是听都没听说过，稀里糊涂的就在老师的安排下跟着别的院一个师兄学习生信， 刚见面师兄就跟我讲解了什么是生物信息学，学生信的种种好处，什么不用做实验就能发文章啊，做的好的话读博能去一个不错的实验室啊之类的，当时我是处于懵逼状态，心想“这是一种怎样的操作？”  在师兄的推荐下，我掏了300大洋买了《DNA和蛋白质序列数据分析工具》《鸟哥的私房菜上下》三本书，用来了解和入门什么是生物信息学。说实话，到现在为止这三本书我都没怎么翻过，在我带入门的人里，我也不会推荐《鸟哥的私房菜》这种，这只会让他们望而却步，只要他们想学，我手里也有他们学不完的资源。好的是，师兄手把手的教了我一段时间。记得师兄说过，要想学生信，就必须学会一门编程语言，否则出门千万别说自己是搞生信的，丢人！所以我就又入了编程的坑，在师兄的强烈推荐下，我选择学习Python，师兄也帮我装上Python.2.7和编辑器，并装上Biopython包，然后扔给我一本全英的《Biopython》和一个脚本，说你要是一周之内不能把这个脚本弄懂，就不要学生信了，不适合。当时我差点就一口老血喷出来，心想“老哥，咱先不谈其他，你好歹也给个中文版本的吧，谁跟你这么强，硕士就到英国留学，毕业论文搞个全英的？”无奈，有总比没有强，还好后来我在网上搜到了这本书的中文版。 因此，在我对Python一无所知的情况下，首先学习了biopython，然后买了《python基础教程》《python核心编程》这两本书更进一步的学习。一但入了此坑就很难回头，尽管现在我主要跟着学习生信的老师和一些小伙伴都是Perl大神，也很难把我从这个坑里拉出来。可喜的是，现在python在生信上的应用越来越广，在机器学习方向，python也是处于领先地位，这更加给了我学下去的理由。当然，主学Python之余，R和Perl也是要懂一点的，用来做图和单行操作还是很必要的！ 人生接触的第一个脚本： import re import sys, getopt import operator from Bio import SeqIO from Bio.Seq import Seq from Bio.SeqRecord import SeqRecord from Bio.Alphabet import generic_nucleotide import re import sys, getopt import operator opts, args = getopt.getopt(sys.argv[1:], “c:i:o:”) blast_info = “” out_file = “” for op, value in opts: if op ==”-o”: #ARG_pattern_MIN_7030.fasta…

Got any book recommendations?

Get In Touch