Category Archives: 基础软件

blast++套件，blat，fastqc，flash，solexaQA，NGS-QC-toolkit，SRA-toolkit，fastx-toolkit 待续

二 06

scalpel软件找indel

Posted on by

Scalpel is available here: http://scalpel.sourceforge.net/
文章是： http://www.nature.com/nmeth/journal/v11/n10/full/nmeth.3069.html
很赞的工具！
软件说明书写的也比较详细：http://scalpel.sourceforge.net/manual.html
他提供了3种情况的找INDELs变异，我目前需要用的就是对我的全基因组测序数据来找，所以用single模式：
为了节省对计算资源的消耗，作者建议我单独对每条染色体分别处理。 Continue reading →

十二 15

制作自己的gene set文件给gsea软件

Posted on by

熟悉GSEA软件的都知道，它只需要GCT,CLS和GMT文件，其中GMT文件，GSEA的作者已经给出了一大堆！就是记录broad的Molecular Signatures Database (MSigDB) 已经收到了18026个geneset，但是我奇怪的是里面竟然没有包括cancer testis的gene set，MSigDB的确是多，但未必全，其实里面还有很多重复。而且有不少几乎没有意义的gene set。那我想做自己的gene set来用gsea软件做分析，就需要自己制造gmt格式的数据。因为即使下载了MSigDB的gene set，本质上就是gmt格式的数据而已：http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats#GMT:_Gene_Matrix_Transposed_file_format_.28.2A.gmt.29 Continue reading →

十二 01

GSEA的统计学原理试讲

Posted on by

GSEA这个java软件使用非常方便，只需要根据要求做好GCT/CLS格式的input文件就好了。我以前也写个用法教程：

但说到统计学原理，就有点麻烦了，我试着用自己的思路阐释一下：
假设芯片或者其它测量方法测到了2万个基因，那么这两万个基因在case和control组的差异度量(六种差异度量，默认是signal 2 noise，GSEA官网有提供公式，也可以选择大家熟悉的foldchange)肯定不一样,那么根据它们的差异度量，就可以对它们进行排序，并且Z-score标准化，在下图的最底端展示的就是

Continue reading →

十一 25

用BioNet这个bioconductor包来找 maximal-scoring subgraph

Posted on by
## 此包是为了解决一个难题： maximal-scoring subgraph (MSS) problem ，在一个巨大的复杂网络里面找到significantly differentially expressed subnetworks，就是说，得到了几百个差异基因，去PPI数据库做网络图的时候，发现还是巨大无比，所以需要用这个包来精简我们的网络图。
heuristically的中文意思：启发性地
## 而这个R包可以整合多种数据结果来给一个网络打分，
它整合了PPI网络分析和寻找功能模块的需求。
重点就是根据一个"igraph" or "graphNEL"对象和打分来找最大的MSS
subnet <- subNetwork(dataLym$label, interactome)
module <- runFastHeinz(subnet, scores)
plotModule(module, scores=scores, diff.expr=logFC) #这个就是精简后的我们的网络图。
其实另外一个函数也有类似的功能，dNetFind https://rdrr.io/cran/dnet/man/dNetFind.html

Continue reading →

十一 24

cytoscape五步曲之二：在cytoscape里面生成网络图

Posted on by
通过上一讲大家应该明白了，网络图是为了展现分子之间的连接关系的，并不是一定要用cytoscape来做，只需要根据连接关系给我们的所有点安排一个坐标，然后把相应的线连接起来即可！那么既然我们要学习cytoscape，肯定是要用cytoscape做好第一步，就是根据输入数据来做网络图。
可以先了解一下cytoscape定义好的输入数据，
http://wiki.cytoscape.org/Cytoscape_User_Manual/Network_Formats 当然，其实木有意义！因为我们不可能拿到cytoscape的输入文件（cys格式的），除非是你朋友传给你的。我们肯定是根据txt.csv等分割的文本文件来做网络图。

Continue reading →

十一 24

cytoscape五步曲之一：明白什么是网络图

Posted on by
想了想还是写一个系列教程吧，问的朋友也太多了，主要是因为cytoscape跟python一样，经历了从2到3的进化阵痛过程，而且进化的面目全非了！！！很多人拿着2.x的说明书教程，视频，然后下载的却是3.x版本的cytoscape，真可怕！！！
已经从两万个芯片探测到的基因里面找到了近千个差异基因了，对它们做了GO/KEGG分析还是抓不住重点，看到文献说可以用PPI数据库做network analysis之后找hub基因，也也许可以说明一些问题！
提到 network analysis ，我想起来我以前总结过 R语言画网络图的三部曲，里面讲到过网络分析的基本原理！

Continue reading →

十一 24

cytoscape五步曲之三：安装各种插件

Posted on by

软件安装我就不多说了，直接去官网下载即可，请务必下载3.x版本，我讲的是 最新版教程！

本次讲解如何给cytoscape安装插件，cytoscape本身是一个平台，学者可以在上面开发各种各样功能的插件实现不同的分析需求，类似于R语言这个平台，人们在上面安装包一样。R里面如何安装包我博客讲了4次，基本上看完的人都会懂。而cytoscape不一样，它的插件安装非常简单！非常简单！非常简单！

你只需要去cytoscape的APP中心找到包，如果你打开了cytoscape的界面，那么网页就会有install的字样，非常显眼，点击就自动安装了，这个时候会安装到

C:\Users\jimmy1314\CytoscapeConfiguration\3\apps\installed 这个目录！！~ 在你的电脑里面 jimmy1314 不一样

如果你这个时候并没有打开cytoscape的界面，那么网页就会有download的字样，也是非常显眼，点击就可以下载， 下载之后你需要自己把下载的jar文件放到cytoscape的安装路径，一般默认是

C:\Program Files\Cytoscape_v3.3.0\apps

最后，cytoscape提供了APP中心，就跟苹果手机安卓手机安卓软件一样，直接在cytoscape软件的菜单栏app中心就可以点击安装！

我要说的就是这么多了，我安装了十几个插件了，都没有什么问题，如果大家有遇到安装不了的，随时报告我，我来更新教程！联系jmzeng1314@163.com 

下面的链接选择性观看：

http://wiki.cytoscape.org/Cytoscape_3/UserManual

十一 14

仅仅对感兴趣的基因call variation

Posted on by

有这个需求，是因为我们经常对某些细胞系进行有针对性的设计变异，比如BAF155的R1064K呀，H3F3A的K27呀，那我我们拿到高通量测序数据的时候，就肯定希望可以快速的看看这个基因是否被突变成功了。现在比对几乎不耗费什么时间了，但是得到的sam要sort的时候还是蛮耗费时间的。假设，我们已经得到了所有样本的sort好的bam文件，想看看自己设计的基因突变是否成功了，可以有针对性的只call 某个基因的突变！

1

Continue reading →

十一 12

仔细探究picard的MarkDuplicates 是如何行使去除PCR重复reads功能的

Posted on by

本帖紧跟前面的仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的

同样的我们也是分单端和双端测序来看结果，并且比较两个工具的区别！

首先对于那个单端数据，samtools给出的结果是：[bam_rmdupse_core] 25 / 53 = 0.4717 in library Continue reading →

十一 12

仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的

Posted on by
在做这个去除PCR重复reads时候必须要明白为什么要做这个呢？WGS？WES？RNA-SEQ?CHIP-SEQ？都需要吗？随机打断测序才需要？特异性捕获不需要？
搞明白了，我们就开始做，首先拿一个小的单端测序数据比对结果来做测试！
samtools rmdup -s tmp.sorted.bam tmp.rmdup.bam
[bam_rmdupse_core] 25 / 53 = 0.4717 in library
我们的测试数据里面有53条records根据软件算出了25条reads都是PCR的duplicate，所以去除了！

Continue reading →

十一 07

mysql的table居然有最大列限制

Posted on by
想着把TCGA的RPKM值矩阵表格写入到mysql，然后做一个查询网页给生物学家，我下载的是所有TCGA收集的mRNA表达数据集数据集-GSE62944 ，共9264个癌症样本，和741个正常组织的表达数据。当我想写入癌症表达矩阵的时候，报错了：
Error in .local(conn, statement, ...) :
could not run statement: Too many columns
简单搜索了一下，发现是mysql有最大列的限制，但是我不是很懂计算机，所以没太看明白该如何调整参数使得mysql列限制扩充：http://dev.mysql.com/doc/refman/5.7/en/column-count-limit.html 所以就把癌症表达矩阵根据癌症拆分了，癌种数量如下：

Continue reading →

十一 03

热图最佳实践-pheatmap

Posted on by

用pheatmap来绘图首先要安装这个包，它就一个功能，画出热图即可，号称是pretty heatmap，的确比其它的好用很多。我以前写过《一步一步学习heatmap.2》的教程，很简单的那种，所以就没有公布在博客上面，结果发现很多其它博客居然能先我一步发出。其实包括本次的pheatmap指南，都没什么好发，的在R里面也是傻瓜式出图，无法就是自己熟练一下参数而已，又不是开发一个包，没什么技术含量。我这里单独提一下pheatmap是因为它的确非常好用，将会是我画热图的不二之选。比如下面这个，是我最喜欢的： Continue reading →

十 31

用samtools idxstats来对de novo的转录组数据计算表达量

Posted on by

de novo的转录组数据，比对的时候一般用的是自己组装好的trinity.fasta序列(挑选最长蛋白的转录本序列)来做参考，用bowtie2等工具直接将原始序列比对即可。所以比对 sam/bam文件本身就包含了参考序列的每一条转录本序列ID，直接对 sam/bam文件进行counts就知道每一个基因的表达量啦！

本来我是准备自己写脚本对sam文件进行counts就好，但是发现了samtools自带这样的工具：http://www.htslib.org/doc/samtools.html 

如果是针对基因组序列，那么这个功能用处不大，但是针对转录本序列，统计出来的就是我们想要的转录本表达量。 Continue reading →

十 22

使用trimmomatic对illumina数据做质控-去接头还有去除低质量碱基

Posted on by
因为一直拿到的是公司给的特别好的数据，所以没太关注质控这个问题，最近拿到了raw data，才发现其实里面的门道挺多的。前面都是用cutadapt这个python软件来去除接头的，但是它有一个弊端，需要自己指定接头文件。正好朋友推荐了trimmomatic，是java软件，所以直接Google找到其官网，然后下载二进制版本解压即可使用！
反正对我的illumina测序数据来说，直接用它就可以把raw data 变成 clean data啦！
1

Continue reading →

十 15

R一大利器之对象的操作函数查询

Posted on by

对于生物出身的部分生物信息学工程师来说，很多计算机概念让人很头疼，尤其是计算机语言里面的高级对象。我以前学编程的时候，给我一个变量，一个数据，一个hash，我就心满意足了，可以解决大部分我数据处理问题，可事情远比想象之中复杂。因为很多高手喜欢用封装，代码复用，喜欢用高级对象。在R的bioconductor里面尤其是如此，经常会遇到各种包装好的S3，S4对象，看过说明书，倒是知道一些对象里面有什么，可以去如何处理那些对象，提取我们想要的信息，比如我就写过一系列的帖子：

Continue reading →