质控之前我们在直播八的时候分析过，公司也给了我质控后的的数据，但是毕竟是别人做的，我们做为一个数据分析师，自己动手来验证一下公司给出的报告也是再好不过的了。大家可以跟着我先将下载数据进行一下质控。 Continue reading
十二
09
十二
09
【直播】我的基因组（十一）:测序数据的比对
上一次直播中，我们对拿到手的测序数据进行了质控，测序数据的质量已经得到了保证。那么接下来就可以把它拿来与参考基因组比对了，这里我们先用参考基因组hg19，大家可以参照【直播】我的基因组（五）:测试数据及参考基因组的准备来下载参考基因组hg19，我这里选择的是UCSC提供的hg19。然后安装bwa软件进行比对，可以参考【直播】我的基因组（四）：计算资源的准备来安装，以及对hg19建立索引。 Continue reading
十二
09
【直播】我的基因组（九）:拿到数据后要做的事情
时隔好几个月，因为各种各样的原因数据终于拿到了自己的手上，真是不容易啊！
拿到数据后，第一件要做的事情就是检查数据传输的完整性，然后备份！我拿到的数据如下：
可以看到，公司给了我测序仪的下机数据(raw data)和他们质控后的clean data，这个过程减少了6G的数据量，对应着约90亿bp的碱基，相当于减少了3个人的全基因组数据。具体推算公式见前面的系列直播贴！
首先我把数据拷贝到了我上上周买的2T移动硬盘里面，再拷贝到我工作电脑一份，服务器一份，私人电脑一份，另外一个移动硬盘一份。然后删除了公司寄给我的硬盘里面的数据，再把硬盘寄回给公司，然后监督他们删除我所有的数据。(做这么多就是为了保护隐私，当然这个大前提是我已经确定数据没有问题了。)
检查数据传输的完整性就是md5校验，看看数据在拷贝过程中有没有意外的损坏（这个在之前下载数据的时候我也说过）！一般传输数据之前，会用md5命令来生成各个文件的md5值，就是下面的MD5.txt文件里面的内容，然后传输数据之后，需要自行用md5sum -c MD5.txt 来校验文件里面记录的文件的完整性，如果显示都是OK，说明文件拷贝传输过程是没有问题的！但这个过程会耗费大量的磁盘读写，磁盘读写能力是有限的，所以开多个进程并不能加快这一过程。
然后我把公司处理好的bam文件上传到服务器做下游分析,我用的winscp软件把文件传到服务器上的！
从明天起，我们就开始正式对基因组进行分析啦！欢迎围观！
请扫描以下二维码关注我们，获取直播系列的所有帖子！
十二
01
java版本GSEA软件的ES score图片的修改
首先要明白这个ES score图片里面的数据是什么，这样才能修改它，因为java是一个封闭打包好的软件，所以我们没办法在里面修改它没有提供的参数，运行完GSEA，默认输出的图就是下面这样： Continue reading
十二
01
GSEA的统计学原理试讲
GSEA这个java软件使用非常方便，只需要根据要求做好GCT/CLS格式的input文件就好了。我以前也写个用法教程：
但说到统计学原理，就有点麻烦了，我试着用自己的思路阐释一下：
假设芯片或者其它测量方法测到了2万个基因，那么这两万个基因在case和control组的差异度量(六种差异度量，默认是signal 2 noise，GSEA官网有提供公式，也可以选择大家熟悉的foldchange)肯定不一样,那么根据它们的差异度量，就可以对它们进行排序，并且Z-score标准化，在下图的最底端展示的就是
十二
01
吐血推荐snpedia数据库，非常丰富的snp信息记录
正好，我拿到了自己的全基因组测序数据，而前些天看到朋友圈推送的文章提到有研究表明STAT4上的rs7574865和HLA-DQ的 rs9275319是国人群中乙型肝炎病毒（HBV）相关肝细胞癌（HCC）遗传易感基因，我就想顺便看看自己在这两个位点的变异情况。一般的流程是先找完变异位点，然后用vep/snpEFF对变异位点进行注释，然后看看有没有这两个位点。但我仅仅是想查看这两个位点，所以我会根据它的rsID来找到它的基因组坐标，再直接call这个位置的变异情况。以前我都是用dnSNP来查看rsID的基因组坐标的，
mkdir -p ~/annotation/variation/human/dbSNPcd ~/annotation/variation/human/dbSNP## https://www.ncbi.nlm.nih.gov/projects/SNP/## ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh38p2/## ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/nohup wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/VCF/All_20160601.vcf.gz &wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/VCF/All_20160601.vcf.gz.tbi
十一
29
如何安装别人开发的未发表的包
我以为我写完了R包终极解决方案！ 之后，应该不会再有任何关于R包安装的问题产生了，但仔细回过头来看才发现，我介绍的都是如何从CRAN或者bioconductor里面安装正规发布的包，但是有很多人开发的是自己私人的包，而我们有的确非常需要用怎么办？？这个时候就需要下载别人开发的包来安装了。比如我R包地址见github：https://github.com/jmzeng1314/humanid Continue reading
十一
29
如何开发自己的R包
随着R语言的流行度的提高，开发一个R包已经不再是专业程序猿才有的技能了。我这里讲的不是如何写一个包含了复杂统计公式或者发表一篇SCI文章的包，而是简简单单的用Rstudio自带的创建包的功能把自己的几个函数和数据打包！！！我R包地址见github：https://github.com/jmzeng1314/humanid Continue reading
十一
29
R来完成表达芯片分析全流程
包括如何从GEO下载数据，如何分组，两组直接如何找差异，差异基因如何去注释，包括GO/KEGG注释，还有特殊数据库，自定义数据库的注释，比如oncogene或者tumor suppress genes，TF的gene注释，还有GSEA软件的分析。
然后是对选择好的差异基因去string等PPI数据库拿到网络数据，在R或者cytoscape里面画网络图，然后是用MCODE插件和bioNet包来对网络找sub-network或者module，和hub genes。
就拿GSE42872 这个数据来做例子吧，希望听众具有基础R知识，了解什么是bioconductor，然后具有基础生物学知识，知道什么是基因，什么是表达，什么是通路，什么是富集，什么是注释。
总共10讲，每次半小时，每周3,4,6的晚上十一点开讲！
讲义的草稿如下，如果你能看懂草稿，能自己学会，就不用参加本次课程啦。
如果需要问我问题，就赶快找我申请加入交流群，提供本次培训的全套视频和代码！！！
十一
28
R语言画网络图三部曲之igraph
经过热心的小伙伴的提醒，我才知道我以前写的R语言画网络图三部曲竟然漏掉了最基础的一个包，就是igraph，不了解这个，后面的两个也是无源之水。
R语言画网络图三部曲之networkD3
R语言画网络图三部曲之sna
其实包括了3个包：igraph/RBGL/Rgraphviz
用到了一个测试数据，是构建好的PPI网络对象：We will first analyse a curated data set of protein-protein interactions in the yeast Saccharomyces cerevisiae extracted from published papers. This data set comes from with an R package called “yeastExpData”, which calls the data set “litG”. This data was first described in a paper by Ge et al (2001) in Nature Genetics (http://www.nature.com/ng/journal/v29/n4/full/ng776.html).
十一
25
hisat2+stringtie+ballgown
早在去年九月，我就写个博文说 RNA-seq流程需要进化啦！ http://www.bio-info-trainee.com/1022.html ，主要就是进化成hisat2+stringtie+ballgown的流程，但是我一直没有系统性的讲这个流程，因为我觉真心木有用。我只用了里面的hisat来做比对而已！但是群里的小伙伴问得特别多，我还是勉为其难的写一个教程吧，你们之间拷贝我的代码就可以安装这些软件的！然后自己找一个测试数据，我的脚本很容易用的！ Continue reading
十一
25
用BioNet这个bioconductor包来找 maximal-scoring subgraph
## 此包是为了解决一个难题： 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
十一
25
最终还是把博客的全半角中英文标点符合的bug解决了
已经有非常多的小伙伴跟我反映了直接拷贝我的代码无法运行的问题，其实报错的原因很简单，就是中英文标点的bug而已。所以我给大家的理由是不用那么懒，拷贝我的代码，我就是希望你们能手动敲打每一个命令，来熟练记忆使用。
其实，我没那么好心，我就是懒而已。因为这个博客是host在阿里云的免费服务器上面的，各种IP密码我懒得去记忆，就差不多忘记了。当初弄好了我就懒得管了，正好现在博客免费时期快到了，也就想把这个问题解决掉。
我很简单搜索了一下，http://shiyun1013.blog.163.com/blog/static/10774036201301824446708/ 需要连接我博客的ftp，去修改博客里面的文件，
<?php remove_filter('the_content','wptexturize'); ?> 正好看看这个标点符号被改变了吗？
好像还不错，以后大家就可以直接拷贝我的代码去运行啦！
下面是我登陆了ftp，发现以前用rmarkdown写的几个教程，感兴趣的小伙伴可以随便看看！
十一
25
kegg在线链接图的颜色设置
一般来说， 有了kegg的ID，就可以直接去官网查看具体的通路图片，但是需要把差异基因给标注上去，就有点麻烦了，我以前做过类似的工作，结果没有做笔记，这次相当于重新造了个轮子，好惨！
简单的KEGG图片，看下面的url：
可以看出来只需要变化KEGG的ID即可。
如果要做下面的这个，上调基因用红色表示，下调基因用绿色表示：
十一
24
cytoscape五步曲之二：在cytoscape里面生成网络图
通过上一讲大家应该明白了，网络图是为了展现分子之间的连接关系的，并不是一定要用cytoscape来做，只需要根据连接关系给我们的所有点安排一个坐标，然后把相应的线连接起来即可！那么既然我们要学习cytoscape，肯定是要用cytoscape做好第一步，就是根据输入数据来做网络图。
可以先了解一下cytoscape定义好的输入数据，
http://wiki.cytoscape.org/Cytoscape_User_Manual/Network_Formats 当然，其实木有意义！因为我们不可能拿到cytoscape的输入文件（cys格式的），除非是你朋友传给你的。我们肯定是根据txt.csv等分割的文本文件来做网络图。
十一
24
cytoscape五步曲之一：明白什么是网络图
想了想还是写一个系列教程吧，问的朋友也太多了，主要是因为cytoscape跟python一样，经历了从2到3的进化阵痛过程，而且进化的面目全非了！！！很多人拿着2.x的说明书教程，视频，然后下载的却是3.x版本的cytoscape，真可怕！！！
已经从两万个芯片探测到的基因里面找到了近千个差异基因了，对它们做了GO/KEGG分析还是抓不住重点，看到文献说可以用PPI数据库做network analysis之后找hub基因，也也许可以说明一些问题！
提到 network analysis ，我想起来我以前总结过 R语言画网络图的三部曲，里面讲到过网络分析的基本原理！
十一
24
cytoscape五步曲之三：安装各种插件
软件安装我就不多说了，直接去官网下载即可，请务必下载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中心就可以点击安装！
一般来说每一个插件都会对应的文章发表，而且一定在cytoscape的APP中心有展示，如果没有，可能是冒牌的哦！比如下面的插件
我要说的就是这么多了，我安装了十几个插件了，都没有什么问题，如果大家有遇到安装不了的，随时报告我，我来更新教程！联系jmzeng1314@163.com
下面的链接选择性观看：
十一
23
关于multiple mapping我想说的
很多时候，我们都要选取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到底对数据做了什么？
提到normalization很多人都烦了，几十种方法，而对于芯片或者其它表达数据来说，最常见的莫过于quantile normalization啦。那么它到底对我们的表达数据做了什么呢？首先要么要清楚一个概念，表达矩阵的每一列都是一个样本，每一行都是一个基因或者探针，值就是表达量咯。quantile normalization 就是对每列单独进行排序，排好序的矩阵求平均值，得到平均值向量，然后根据原矩阵的排序情况替换对应的平均值，所以normalization之后的值只有平均值了。具体看下面的图： Continue reading
十一
23
用R的bioconductor里面的stringDB包来做PPI分析
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进行画图