十一 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 →

十一 25

最终还是把博客的全半角中英文标点符合的bug解决了

Posted on by

已经有非常多的小伙伴跟我反映了直接拷贝我的代码无法运行的问题，其实报错的原因很简单，就是中英文标点的bug而已。所以我给大家的理由是不用那么懒，拷贝我的代码，我就是希望你们能手动敲打每一个命令，来熟练记忆使用。
其实，我没那么好心，我就是懒而已。因为这个博客是host在阿里云的免费服务器上面的，各种IP密码我懒得去记忆，就差不多忘记了。当初弄好了我就懒得管了，正好现在博客免费时期快到了，也就想把这个问题解决掉。

我很简单搜索了一下，http://shiyun1013.blog.163.com/blog/static/10774036201301824446708/ 需要连接我博客的ftp，去修改博客里面的文件，

<?php remove_filter('the_content','wptexturize'); ?>  正好看看这个标点符号被改变了吗？
好像还不错，以后大家就可以直接拷贝我的代码去运行啦！
下面是我登陆了ftp，发现以前用rmarkdown写的几个教程，感兴趣的小伙伴可以随便看看！

Continue reading →

十一 25

kegg在线链接图的颜色设置

Posted on by
一般来说， 有了kegg的ID，就可以直接去官网查看具体的通路图片，但是需要把差异基因给标注上去，就有点麻烦了，我以前做过类似的工作，结果没有做笔记，这次相当于重新造了个轮子，好惨！
简单的KEGG图片，看下面的url：
如果要做下面的这个，上调基因用红色表示，下调基因用绿色表示：

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

十一 23

关于multiple mapping我想说的

Posted on by

很多时候，我们都要选取unique mapped的reads，尤其是在RNA-seq和CHIP-seq的时候，但是如何保留，各种教程都不一致，我稍微总结了一下，是因为使用的比对工具不一样导致的！但是主要都反应在sam文件的一系列tag里面~

首先对bwa来说，如果它遇到一个reads可以比对到参考基因在的多个序列，只会随机的选取一个位置来输出到sam文件，但是会加上一个tag是XS:I:<N>来告诉我们第二好的比对情况的比对得分是多少，bowtie也是一样。但是它们都有参数来决定是否只对每个reads输出一条信息，还是输出全部的信息，在bwa是-a的参数，在bowtie里面是-m参数。

但是bowtie2里面取消了这个参数，它们都必须用XS:I:<N>这个tag来挑选unique mapped的reads

但是如果是用hisat来比对的话，决定是否是唯一比对的却是NH这个tag信息。默认情况下一条reads可以输出多条比对结果。

我想起了再补充吧，其实应该找几个例子用IGV看看，就明白了，可是我暂时没有时间了，只是觉得这个很重要，就提一下。

 

十一 23

quantile normalization到底对数据做了什么？

Posted on by

提到normalization很多人都烦了，几十种方法，而对于芯片或者其它表达数据来说，最常见的莫过于quantile normalization啦。那么它到底对我们的表达数据做了什么呢？首先要么要清楚一个概念，表达矩阵的每一列都是一个样本，每一行都是一个基因或者探针，值就是表达量咯。quantile normalization 就是对每列单独进行排序，排好序的矩阵求平均值，得到平均值向量，然后根据原矩阵的排序情况替换对应的平均值，所以normalization之后的值只有平均值了。具体看下面的图： Continue reading →

十一 23

用R的bioconductor里面的stringDB包来做PPI分析

Posted on by

PPI本质上是根据一系列感兴趣的蛋白质或者基因（可以是几百个甚至上千个）来去PPI数据库里面找到跟这系列蛋白质或者基因的相互作用关系！

本次的主角是stringDB，顾名思义用得是大名鼎鼎的string数据库，
本来还以为需要自己上传自己的基因给这个数据库去做分析，没想到他们也开发了R包，主页见： http://www.bioconductor.org/packages/release/bioc/html/STRINGdb.html 而我比较喜欢用编程来解决问题，所以就学了一下这个包，非常好用！
它只需要一个3列的data.frame，分别是logFC,p.value,gene ID,就是标准的差异分析的结果。
然后用string_db$map函数给它加上一列是 string 数据库的蛋白ID，然后用string_db$add_diff_exp_color函数给它加上一列是color。
用string_db$plot_network函数画网络图，只需要 string 数据库的蛋白ID，如果需要给蛋白标记不同的颜色，需要用string_db$post_payload来把color对应到每个蛋白，然后再画网络图。
也可以直接用get_interactions函数得到所有的PPI数据，然后写入到本地，再导入到cytoscape进行画图

Continue reading →

十一 23

java环境变量的问题

Posted on by
有篇文章提到了cytoscape，想着一直没用过这个神器对不起我生信大神的称号呀，就下载了准备安装，居然报错了，简直不可思议，因为一直以为它是java软件，一般不需要安装，结果是exe的，只是依赖于java，报错是EXE4J_JAVA_HOME, No JVM could be found on your system，这是个很常见的错误，我 简单搜索了解决方案https://wincrunch.com/exe4j-java-home-no-jvm-could-be-found-on-your-system/ 居然无效，但是里面有句话引起了我的注意，通常64位的window电脑的java是安装在Program Files 而不是Program Files (x86)，这才是问题所在，我当初图简单，直接用了JDK来安装JRE，所以导致软件安装目录错误。有非常多的生物信息学软件都依赖与java，比如IGV,GSEA,cytoscape，一般来说window电脑安装好了java之后这些软件都挺好用的。那么关于java问题，我整理了3个：

Continue reading →

十一 23

【直播】我的基因组（八）:原始测序数据质量报告

Posted on by

由于我是分期付款，所以我先拿到了我的测序数据的质控结果和比对情况分析报告，需要补齐全款后才能拿到原始测序数据！（中间还出了个小意外，打款的时候不小心多打了30块钱！(⊙o⊙)…不过多打的30块钱想拿回来估计不太可能了，需要填写书面申请表格并且自费快递到公司，这边跨境快递费都不止这个数了） Continue reading →

十一 23

【直播】我的基因组（七）:从整体理解全基因组测序数据的变异位点

Posted on by
首先记住一个很重要的知识点，变异是相对的！

简单说一下什么是找变异，变异跟突变有什么区别呢？举个栗子：有国际组织规定了人类的参考基因组（如UCSC,ENSEMBL,NCBI等，前面帖子都有讲)，就是 AAAAA(这里简化一下，就5个碱基，其实人类基因组多达30亿个)  。现在通过给自己测序得知，我与之对应的是AGCAA，那么我相比国际基因组来说，就是2个变异位点，位于基因组的坐标2和3，但是它们还不能说就是突变。 Continue reading →

十一 23

【直播】我的基因组（六）:变异位点注释数据库的准备

Posted on by

通常一个人的全基因组测序数据可以挖掘到四百万个SNVs(跟参考基因组不一样的单碱基位点)，还有五十万的indels(insertions or deletions),但是得到的数据通常是以vcf文件格式给出的(自行搜索什么是vcf格式)，比如下面：

很明显，正常人是看不懂这些变异位点有啥子一样的，只知道第20条染色体的1230237坐标上面本来是一个T碱基的，但是突变成了G，那么我们必然还想知道，这个位点是在某个基因上面吗？如果是，在基因的外显子还是内含子？它的突变有没有改变该基因的功能呢？有没有影响它的转录和翻译呢？还有世界上有没有其他正常人也是这个位点变异呢？如果有，是哪些人种呢？有没有癌症病人也发现了这个变异呢？如果有，是什么癌症呢？所以我们必须下载一系列的变异位点注释数据库，来全方位的解释我们自己找到那四百万个SNVs和五十万的indels。下面我们一起进行数据库准备。

Continue reading →

十一 15

htseq-counts跟bedtools的区别

Posted on by

我以前写过bedtools和htseq-counts的教程，它们都可以用来对比对好的bam文件进行计数，正好群里有小伙伴问我它们的区别，我就简单做了一个比较，大家可以先看看我以前写的软件教程。写的有的挫：

使用Bedtools对RNA-seq进行基因计数 ，

转录组HTseq对基因表达量进行计数

言归正传，我这里没精力去探究它们的具体原理，只是看看它们数一个read是否属于某个基因的时候，区别在哪里，大家看下图： Continue reading →

十一 14

TPM值就是RPKM的百分比嘛！

Posted on by

很久以前就有人问过这个问题啦，虽然目前主流还是用RPKM/FPKM来形容一个基因的表达量。但是既然大家都说TPM更好，我也来探究一下吧！

我不喜欢看公式，直接说事情，我有一个基因A，它在这个样本的转录组数据中被测序而且mapping到基因组了 5000个的reads，而这个基因A长度是10K，我们总测序文库是50M，所以这个基因A的RPKM值是 5000除以10，再除以50，为10. 就是把基因的reads数量根据基因长度和样本测序文库来normalization 。 Continue reading →

十一 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 →

十一 11

数据库批量注释不可盲目-annovar数据库错误

Posted on by

我对H3F3A这个基因做了两个突变的cellline，分别是G34V和K27M，现在知道这个基因在hg38上面的坐标是：

Genomic Location for H3F3A Gene
Chromosome:  1
Start:226,061,851 bp from pter  End:226,072,002 bp from pter
Size:10,152 bases    Orientation:Plus strand

然后我用samtools结合bcftools把该基因区域的snp位点call出来：

samtools mpileup -r chr1:226061851-226072001 -t "DP4" -ugf ~/reference/genome/hg38/hg38.fa  *sorted.bam | bcftools call -vmO z -o  H3F3A.vcf.gz

Continue reading →

十一 10

一个基因在同一套基因组上面竟然有两个定位！

Posted on by

查了好久的bug，终于搞清楚问题所在了！因为要对基因进行reads计数，所以要拿到基因在基因组上面的染色体起始终止坐标，结果发现了个十分诡异的现象，很多基因有多个坐标，比如下面这个PTPRS 在hg38这个基因组版本，居然有两个定位，因为我是写程序格式化得到的坐标，所以我check了我的程序，http://www.biotrainee.com/thread-472-1-1.html  感兴趣的同学可以点开看看我的代码！

tmp
Continue reading →