学术咖

Category: Bioinformatics

  • WGCNA（my project）

    setwd(“D:/Ma_transcription/WGCNA/”) database <- read.table(file = “COS_deseq_counts_normalized.txt”, sep = “\t”, header = T, row.names = 1, stringsAsFactors = F) library(reshape2) library(‘WGCNA’) enableWGCNAThreads()#打开多线程 WGCNA_matrix = t(database[order(apply(database,1,mad), decreasing = T)[1:10000],]) subname=sapply(colnames(database),function(x) strsplit(x,”_”)[[1]][1]) datTraits = data.frame(gsm=names(database), subtype=subname) rownames(datTraits)=datTraits[,1] head(datTraits) # Choose a set of soft-thresholding powers powers = c(c(1:10), seq(from = 12, to=20, by=2)) datExpr <- WGCNA_matrix # Call…

  • DESeq2使用流程

    library(“DESeq2”) database

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

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

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

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

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