昨天，我们了解了一下SAM格式的比对结果，不知道大家理解的怎么样。但是全基因组测序数据实在是太大了，即使比对后把sam文件压缩成二进制的bam文件也还有55G（如何压缩转换可查看直播十二），如果完整的导入IGV查看会略微考验计算机配置。
十二
19
十二
19
【直播】我的基因组（13）:了解sam格式比对结果
比对工具到现在已经多如牛毛了，见列表： https://en.wikipedia.org/wiki/List_of_sequence_alignment_software 。但是能被大多数人熟知的，就是bowtie和bwa（我们在十一讲中用的才是bwa），它们把测序数据比对到参考基因组之后，都会生成一个sam格式的文件。随后的大部分分析都是基于sam格式进行的分析，虽然Jimmy多次强调这些基础知识的重要性需要大家私下自学。但是由于这个sam文件实在是太重要了！！！所以，不得不亲自抽出一讲来说说它，后面也会基于此写十多篇文章： Continue reading
十二
15
制作自己的gene set文件给gsea软件
熟悉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
十二
15
CpG Islands记录文件下载的4种方式
这个也是读者来信最多的，关于基因组某些区域的起始终止坐标的下载问题，genomic feature的问题，一般是gtf文件或者bed文件，比如人类hg19上面的所有外显子的坐标记录文件，所有基因的坐标记录文件，所有lncRNA，rRNA等等，我这里拿CpG Islands记录文件下载的4种方式举例子给大家说明一下： Continue reading
十二
15
deeptools辅助CHIP-seq数据分析-可视化
有很多读者来信，CHIP-seq数据比对后的bam文件如果根据基因组的所有基因来画热图，profile图呢？
十二
11
gene symbol 中的奇怪开头基因
这本是我为论坛的基础板块写的一个基础知识点，但是浏览量实在有限，不忍它蒙尘，特在博客重新发布一次！原帖见：http://www.biotrainee.com/thread-511-1-1.html
gene symbol 是非常官方的，由HUGO 组织负责维护，有专门的数据库HGNC database of human gene names | HUGO
以前分析数据的时候，有一些基因的symbol很奇怪，让我百思不得其解，比如
C orf 系列基因，
HS.系列基因，
KRTAP系列基因，
LOC系列基因，
MIR系列基因，
LINC系列基因
它们往往一个系列，就有好几百个基因；
C12orf44; Chromosome 12 Open Reading Frame 44; 这个是C orf系列基因的意思
MIR系列基因应该是 miRNA相关的基因
LINC系列基因应该就是long intergenic non-protein coding RNA
LOC系列基因，是非正式的，推定的，日后可能被更合适的名字替代
我这里做好了所有的基因对应关系，去生信菜鸟团QQ群里下载吧，共47938个基因的symbol和entrez gene id还有name，还有alias的对应!
还有一些RNA基因，根本就没有symbol，比如：CTA/B/C/D系列的
Aliases for ENSG00000271971 Gene
Quality Score for this RNA gene is 1
Aliases for ENSG00000271971 Gene
CTD-2006H14.2 5
External Ids for ENSG00000271971 Gene
Ensembl: ENSG00000271971
还有，如果你看到HS.开头的基因，它是unigene的ID了，已经不再是symbol啦。
十二
11
用R获取芯片探针与基因的对应关系三部曲-NCBI下载对应关系
这是系列文章，请先看：
用R获取芯片探针与基因的对应关系三部曲-bioconductor
ncbi现有的GPL已经过万了，但是bioconductor的芯片注释包不到一千，虽然bioconductor可以解决我们大部分的需要，比如affymetrix的95,133系列，深圳1.0st系列，HTA2.0系列，但是如果碰到比较生僻的芯片，bioconductor也不会刻意为之制作一个bioconductor的包，这时候就需要自行下载NCBI的GPL信息了，也可以通过R来解决：
##本质上是下载一个文件，读进R里面，然后解析行列式，得到芯片探针与基因的对应关系，看下面的代码，你就能理解了。 Continue reading
十二
10
解决阿里云博客的虚拟主机升级问题
首先感谢生信菜鸟团的各个小伙伴的大力支持，在阿里云的2年免费虚拟主机到期之后，我成功了续费了，但是坑爹的阿里云居然把我的IP地址和用户名都给替换了，导致了一些莫名其妙的bug。
虽然只是warning，不影响网站访问，但实在是影响界面美观,如下：
十二
09
【直播】我的基因组（十二）:先粗略看看几个基因吧
昨天我们说到，测序得到的fastq文件map到基因组之后，我们通常会得到一个sam或者bam为扩展名的文件。SAM的全称是sequence alignment/map format，而BAM就是SAM的二进制文件。通常sam文件太大，我们会生成bam文件来节省空间。sam文件和bam文件的转换用samtools这个软件就可以完成。 Continue reading
十二
09
【直播】我的基因组（十）:测序数据质量控制
质控之前我们在直播八的时候分析过，公司也给了我质控后的的数据，但是毕竟是别人做的，我们做为一个数据分析师，自己动手来验证一下公司给出的报告也是再好不过的了。大家可以跟着我先将下载数据进行一下质控。 Continue reading
十二
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