使用Cytoscape分析和可视化网络(1)

获得网络数据
获得网络数据的方法有三种,第一种是通过检索蛋白相互作用数据库构建网络。
Cytoscape欢迎界面:
File-Welcome
如果启动时,没有显示欢迎界面,可以使用菜单”help” 》 “show welcome screen”打开。在欢迎界面,选择”From Network Database”,将弹出如下界面:
Cytoscape_from_database
在搜索框中输入基因或者蛋白的名称,点击”Search”, 在结果框中选择有记录的数据库,”Automatic Network Merge”(自动合并网络)可以选择,也可以手动合并网络,最后点击”Import”。下图是其中得到的一个网络:
Cytoscape_KCTD10

第二种得到网络的方法是通过Agilent Literature Search插件进行文本挖掘,通俗地来讲就是让机器帮你读文献。
安装Agilent Literature Search插件, 菜单”Apps” 》 “APP Manager”, 选择”Agilent Literature Search”, 点击”Install”进行安装,如下图所示:
Cytoscape_plugin_install
安装成功后,会在菜单”Apps”下,看到”Agilent Literature Search”,点击使用它,如下图所示:
Cytoscape_text_mining
“Terms”框这里依然输入”KCTD10″,”Max Engine Matches(最多读多少篇文献)”这里设置为100,其他参数这里不做修改,你可以根据需要自行设置,点击绿三角图标进行文本挖掘。 下图是通过Agilent Literature Search插件文本挖掘得到的网络:
KCTD10_text_mining

得到网络的第三种方法是从网络文件中导入,这里先介绍一下Cytoscape支持的网络文件格式。
.sif(Simple Interaction Format,简单相互作用格式)

nodeA interactionType1 nodeB
nodeA interactionType2 nodeC
NodeF

.sif是最直接的格式,用户可以通过文本编辑器定义网络。.sif的每一行包含三个标识,第一个标识是源节点,第二个标识是相互作用类型,第三个标识是目标节点。上面的示例文件表示,nodeA和nodeB通过interactionType1相互作用有一条边,nodeA和nodeC通过interactionType2相互作用有一条边,NodeF定义了一个节点,但是没有边。

.noa(节点文件)

AttributeName
nodeA = value1
nodeB = value2
nodeC = value1

第一行是属性名,例如’SubcellularLocation’,接下来的每一行包含一个节点名、等号、节点属性值,属性值可以是数值型也可以是字符型,例如’4.12’或者’nucleus’。

.eda(边文件)

AttributeName
nodeA (interactionType) nodeB = 0.56
nodeB (interactionType) nodeC = 0.918
nodeB (interactionType) nodeA = 0.3412

边文件跟节点文件类似,这里的边名由源节点、相互作用类型、目标节点组成。

表达值文件格式(Expression data file format)
需要将扩展名改为.mrna或者.pvals才能被Cytoscape识别,改为哪一个扩展名没有区别。

Gene Label Experiment1 Experiment2 Experiment1 Experiment2
geneA labelA valueA_1 valueA_2 pvalueA_1 pvalueA_2
geneB labelB valueB_1 valueB_2 pvalueB_1 pvalueB_2
geneC labelC valueC_1 valueC_2 pvalueC_1 pvalueC_2

表达值数据是一个矩阵,每一行代表一个基因或者蛋白在不同实验条件下的表达值。第一行是列标签,第一列是基因名,第二列是基因的描述,接下来的列是表达值,每一个实验一列。如果数据包含P-value或者其他衡量显著性的值,每一个实验则对应两列,第一列表达值,第二列P-value,这两列的列标签必须相同。使用jActiveModules插件分析功能模块时,表达值数据必须包含0(最显著)到1(最不显著)的P-value。

从文件导入网络:
欢迎界面中,选择”From Network File”;
或者在菜单中,”File” 》 “Import” 》 “Network” 》 “File”。
这里导入示例文件galFiltered.sif,这个文件包含在Cytoscape Nature Portocol的补充信息中。

使用Cytoscape分析和可视化网络(简介)

Cytoscape是一款用于分析和图形化展示网络的软件。Cytoscape最基本的组织形式是网络图,分子(基因或者蛋白)表示一个节点,分子间的相互作用表示一个边。Cytocape软件本身只提供一些基本的功能,更多深入分析的功能需要依靠插件。而Cytocape最具魅力的地方也在于它拥有丰富的插件,例如,用于GO注释的BiNGO,用于功能模块分析的jActiveModules,用于文本挖掘构建网络的Agilent Literature Search……

Cytoscape开源、免费、跨平台,可以从它的官网(http://www.cytoscape.org/)下载,并免费使用。当前最新版本为Cytoscape 3.3.0,需要先安装Java 8。
Cytoscape用户指南:http://opentutorials.cgl.ucsf.edu/index.php/Portal:Cytoscape3

Cytoscape相关文献:
1. Shannon, P., Markiel, A., Ozier, O., Baliga, N.S., Wang, J.T., Ramage, D., Amin, N., Schwikowski, B. and Ideker, T. (2003) Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome research, 13, 2498-2504.
2. Cline, M.S., Smoot, M., Cerami, E., Kuchinsky, A., Landys, N., Workman, C., Christmas, R., Avila-Campilo, I., Creech, M., Gross, B. et al. (2007) Integration of biological networks and gene expression data using Cytoscape. Nature protocols, 2, 2366-2382.
3. Smoot, M.E., Ono, K., Ruscheinski, J., Wang, P.L. and Ideker, T. (2011) Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics, 27, 431-432.
4. Saito, R., Smoot, M.E., Ono, K., Ruscheinski, J., Wang, P.L., Lotia, S., Pico, A.R., Bader, G.D. and Ideker, T. (2012) A travel guide to Cytoscape plugins. Nature methods, 9, 1069-1076.

本教程依据07年和12年的两篇文献来写的,相当于一个中文简略版,建议大家直接读英文原文。

Cytoscape分析和可视化网络,大概可以分为五个步骤:
1. 获得网络数据(Obtain network data)
2. 探索网络并输出(Explore network and generate layout)
3. 用表达信息或者其他功能属性注释(Annotate with attribute and expression data)
4. 分析网络特征(Analyze network features)
5. 基因功能富集(Detect enriched gene functions)

使用VennDiagram包绘制韦恩图

维恩图,也叫文氏图,用于显示元素集合重叠区域的图示。最简单的方法是用PPT画几个圆圈,然后标上数字。R中的VennDiagram包最多可以画5组数据的韦恩图,且图的质量可以达到发表的要求。
VennDiagram文章发表在BMC Bioinformatics杂志上,下图就是文章中的插图。
VennDiagram

安装VennDiagram包

install.packages("VennDiagram")

载入VennDiagram包

library(VennDiagram)

绘制韦恩图

list1 <- c(1,2,3,4,5,6,7,8,9,10)
list2 <- c(6,7,8,9,10,11,12,13,14,15)
list3 <- c(1,2,3,4,5,11,12,13,14,15)
venn.diagram(list(A=list1,B=list2,C=list3), "figure.tiff")

VennDiagram_figure

venn.diagram()函数最基本的两个参数就是输入和输出,输入是一个包含每组数据的list,输出指定文件名即可。venn.diagram()函数是输入是每组数据的详单,有时候我们并不知道每组数据具体元素,只是知道各组数据之间有多少交集,那该怎么办呢?

VennDiagram包提供了draw.single.venn()、draw.pairwise.venn()、draw.triple.venn()、draw.quad.venn()、draw.quintuple.venn()函数,分别绘制1-5组数据的韦恩图。
还是以上面的例子,用draw.triple.venn()作图。

draw.triple.venn(area1=10,area2=10,area3=10,n12=5,n23=5,n13=5,n123=0,category=c("A","B","C"))

area1,area2,area3表示三组数据中元素的数量,n12,n23,n13,n123表示各组数据之间交集的数量,category指定每组数据的名称。

如果要把图画得好看点,可以添加其他的修饰参数。?VennDiagram查看帮助文档。

一个可以在线绘制韦恩图的网址:
http://bioinformatics.psb.ugent.be/webtools/Venn/

RNA-Seq数据分析流程

RNA-Seq分析流程用的最多的是Tuexo workflow,如下图所示,
tuxedo_workflow

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

这套流程最大的优势在于,软件衔接非常好,Tophat比对、cufflinks组装、cuffdiff差异表达,Nature protocol还发表了一篇详细的操作手册。

但是,很多人都说cufflinks推出的FPKM算表达量并不好。cufflinks这帮大牛后面也意识到了问题,推出了cuffquant、cuffnorm、cuffdiff2,做出了一些修正,效果还是差强人意。cufflinks这帮大牛又开发出HISAT来替代Tophat,StringTie替代cufflinks。

HISAT2比对,StringTie组装,使用HTseq提供的脚本重新计算转录本对应的reads数目,用DEseq2进行差异表达分析,是最靠谱的选择。

Python分割fasta文件

CPC可以从转录组得到的fasta序列中分析出哪些是编码RNA,哪些是lncRNA,但是它的速度实在是太慢了。我写了一个脚本,将fasta文件分割成更小的文件,在不同计算节点上面同时运行CPC以提高速度。CPC在balst比对这一步是支持多线程的,但是并不支持跨节点。对于某些不直接支持多线程的程序,将输入文件分割,再同时运行程序,是一种常用思路。

用法:python3 split_fasta.py -i input.fasta -o prefix -x split_number
-i 指定输入的fasta文件
-o 指定输出文件的前缀,默认是split_
-x 指定分割文件中fasta序列的数目,最后一个文件小于等于指定的数目,默认是1

命令:grep “>” file.fasta |wc -l
可以用来查看fasta文件中,有多少条fasta序列

代码如下:
split_fasta.py

# -*- coding: utf-8 -*-
 
import sys, getopt
 
# get parameter
opts, args = getopt.getopt(sys.argv[1:], "hi:o:x:")
X = 1
input_file = ""
prefix = "split_"
for op, value in opts:
    if op == "-i":
        input_file = value
    elif op == "-o":
         prefix = value
    elif op == "-x":
        X = int(value)
    elif op == "-h":
        print("Usage: python3 split_fasta.py -i input.fasta -o prefix -x split_number")
        print("default prefix = split_")
        print("default split_number = 1")
        sys.exit()
 
FA_in_file = open(input_file, "r")
 
# read fasta file to a list
fa_Info = []
fa_Seq = []
fa_Num = -1
for Y in FA_in_file.readlines():
    Y = Y.rstrip()
    if Y[0] == ">":
        fa_Info.append(Y)
        fa_Num = fa_Num + 1
        fa_Seq.append("")
    else:
        fa_Seq[fa_Num] = fa_Seq[fa_Num] + Y
 
# split the fasta list to multipe files
file_Num = (fa_Num + 1)//X + 1
for i in range(file_Num):
    exec(prefix + str(i + 1) + ' = open("' + prefix + str(i + 1) + '.fasta"' + ', "w")')
    start = i * X
    end = (i + 1) * X
    if end > fa_Num + 1:
        end = fa_Num + 1
    for j in range(start, end, 1):
        exec(prefix + str(i + 1) + '.write(fa_Info[j] + "\n")')
        while len(fa_Seq[j]) > 60:
            exec(prefix + str(i + 1) + '.write(fa_Seq[j][:60] + "\n")')
            fa_Seq[j] = fa_Seq[j][60:]
        else:
            exec(prefix + str(i + 1) + '.write(fa_Seq[j] + "\n")')
    exec(prefix + str(i + 1) + '.close()')
 
FA_in_file.close()

转录组差异表达分析工具Ballgown

Ballgown是分析转录组差异表达的R包。

软件安装:
运行R,
source(“http://bioconductor.org/biocLite.R”)
biocLite(“ballgown”)
R会自动安装Ballgown,及相应的依赖包。

Ballgown的输入文件
StringTie使用-B参数直接生成Ballgown的输入文件,Cufflinks的输出结果需要使用Tablemaker生成Ballgown的输入文件。
一个有5个输入文件,分别是:
e_data.ctab,外显子水平表达量;
i_data.ctab,内含子水平表达量;
t_data.ctab,转录本水平表达量;
e2t.ctab,外显子与转录本的对应关系;
i2t.ctab,内含子与转录本的对应关系。

Tablemaker
tablemaker -p 4 -q -W -G merged.gtf -o sample01_output sample_01/accepted_hits.bam
-p 指定线程数
-q 去冗余
-W 运行模式是tablemaker,而不是Cufflinks模式
-G 指定组装好的GTF文件,这个文件由cuffmerge生成
BAM文件是最后一个参数,这个文件由Tophat生产

运行Ballgown
运行R
载入Ballgown包
library(ballgown)
载入数据,并创建一个ballgown项目
bg = ballgown(dataDir=’D:\extdata’, samplePattern=’sample’, meas=’all’)
指定分组,及重复样本数目:
pData(bg) = data.frame(id=sampleNames(bg), group=rep(c(1,0), each=3))
差异表达分析:
stat_results = stattest(bg, feature=’transcript’, meas=’FPKM’, covariate=’group’)

有些地方我写得不是很清楚,有一个英文的Tutorial写得非常好。
https://github.com/alyssafrazee/ballgown

转录组组装工具StringTie

StringTie由约翰霍普金斯大学联合德州大学西南医学中心开发,能够组装转录本并预计表达水平。它应用网络流算法和可选的denovo组装,将复杂的数据集组装成转录本。 StringTie-0

 

 

 

 

 

 

 

 

 

相对于其他拼接软件(Cufflinks, IsoLasso, Scripture,Traph等),StringTie能够拼接出更完整、更准确的基因,并且StringTie采用拼接和定量同步进行,相对于其他方法,其定量结果更加准确。 StringTie-1

文中指出,对于从人类血液中获得的reads,StringTie正确组装了10,990个转录本,而Cufflinks只组装了7,187个。而对于模拟的数据集,StringTie正确组装了7,559个转录本,比Cufflinks的6,310个提高了20%。此外,它的运行速度也比其他组装软件更快。

软件安装
下载StringTie并解压:
tar zxvf stringtie-1.1.2.Linux_x86_64.tar.gz
将HISAT2目录添加到环境变量:
vi ~/.bashrc
在文件末位添加:
export PATH=/home/chenwen/bin/stringtie-1.1.2.Linux_x86_64:$PATH
保存退出
source ~/.bashrc

StringTie的输入文件
HISAT的输出文件为SAM格式,需要经过两步转换成StringTie可以使用的BAM格式:
1. SAM转BAM,并排序:
samtools view -S -b input.sam | samtools sort – input.sorted
2.修改HI标签:
samtools view -h input.bam | perl -ne ‘if(/HI:i:(d+)/) {$m=$1-1; $_ =~ s/HI:i:(d+)/HI:i:$m/} print $_;’| samtools view -bS – > input.correct.bam

运行StringTie
stringtie SRR534294.correct.bam -p 16 -G genes.gtf -B -o stout/transcripts.gtf
-p 指定线程数,默认1
-G 指定参考的转录组注释文件
-B 生成用于Ballgown 分析的文件
-o 指定输出文件
更多参数请查看HISAT2的操作手册:
https://ccb.jhu.edu/software/stringtie/index.shtml?t=manual