十二 09

【直播】我的基因组(十一):测序数据的比对

上一次直播中,我们对拿到手的测序数据进行了质控,测序数据的质量已经得到了保证。那么接下来就可以把它拿来与参考基因组比对了,这里我们先用参考基因组hg19,大家可以参照【直播】我的基因组(五):测试数据及参考基因组的准备来下载参考基因组hg19,我这里选择的是UCSC提供的hg19。然后安装bwa软件进行比对,可以参考【直播】我的基因组(四):计算资源的准备来安装,以及对hg19建立索引。

我们首先简单讲一下为什么要进行比对以及比对过程是怎样的?

可以看到我们到手的测序数据格式是fastq,每条reads都是150个碱基组成,如果只看这fastq,我们没办法得知它的意义,参考基因组那么大(人类约30亿个碱基),这个reads在我们基因组的哪里呢?

简单解释一下,假设人类基因组是123456789,如果我们的测序得到的reads是123,那么我们很明显知道这条reads在基因组的第一个位置,跨越了3个长度,如果我们的reads是567,那么我们也可以看到它在基因组的第5个位置。如果我们看到了一个reads是235567,同样我们也很容易看到它在基因组第2位置,但是在第3个位置参考是4,它却是5,这里可能是测序错误,也可能是这个reads真的变异了!

但是我们的参考基因组远远没有那么简单,30亿个碱基的庞大数目,测序的一条reads也有150个碱基,仅仅用肉眼观察基本不可能判断出它到底在哪里。但并非一定观察不到,如果你有多的不可计的时间及精力的话,手工比对穷极一生来搞定一条reads的比对就很不容易了(当然肯定不会有人这么傻,这里只是说数据量真的很大而已)。然而在我们手上可是有8.9亿条reads,所以我们需要借助计算机来进行比对,现在比较流行的基因组比对工具是bwa和bowtie,它俩的算法不一样,但是我们不需要了解那么具体,只需要知道它可以把我们的fastq测序文件通过与参考基因组的比对生成sam格式(自行搜索了解该格式)的比对结果文件(如下),从sam文件中,我们可以看到每条reads在参考基因组的位置,这条reads是在哪一条染色体,又是在这条染色体的哪个位置就可以一目了然。

对于比对的结果,我们可以用IGV可视化查看,还可以手动查看每个基因的比对情况:

下面我简单讲一下代码

1,下载hg19基因组

cd ~/reference

mkdir -p genome/hg19  && cd genome/hg19

nohup wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz &

tar zvfx chromFa.tar.gz

cat *.fa > hg19.fa

rm chr*.fa

首先要理解linux基础命令,在我们的服务器上面新建好目录,找到hg19的下载链接,用linux自带的wget下载,因为文件太大,所以我们用nohup放在后台下载。下载后是压缩文件 chromFa.tar.gz,在linux里面需要用tar zvfx 来解压tar.gz文件即可。解压开后是一个个文件,需要用cat合并!最终效果如下:

2,安装bwa软件

## Download and install BWA

cd ~/biosoft

mkdir bwa &&  cd bwa

#http://sourceforge.net/projects/bio-bwa/files/

wget https://sourceforge.net/projects/bio-bwa/files/bwa-0.7.15.tar.bz2

tar xvfj bwa-0.7.15.tar.bz2 # x extracts, v is verbose (details of what it is doing), f skips prompting for each individual file, and j tells it to unzip .bz2 files

cd bwa-0.7.15

make

~/biosoft/bwa/bwa-0.7.15/bwa

我所有的软件都安装在自己的home目录下面的biosoft文件夹。同样,也是找的bwa的下载地址,然后解压,然后直接make即可。很多人的服务器会报错zlib.h缺少的问题,看我以前的教程:http://www.bio-info-trainee.com/518.html ,缺少什么你就安装什么,但是缺少的东西需要安装到系统环境变量,但是我的bwa是直接安装到自己的目录,所以我用全路径在调用该软件。如果你的这个命令~/biosoft/bwa/bwa-0.7.15/bwa  能够显示下面的help文档,说明你已经安装成功啦~

3,对hg19参考基因组用bwa构建索引

cd ~/reference

mkdir -p index/bwa && cd index/bwa

nohup time ~/biosoft/bwa/bwa-0.7.15/bwa index   -a bwtsw   -p ~/reference/index/bwa/hg19  ~/reference/genome/hg19/hg19.fa 1>hg19.bwa_index.log 2>&1   &

代码很简单,就是新建好一个文件夹来存放我们的参考基因组的索引,我这里选择的是我的home目录下面的reference/index/bwa/ 文件夹,可以看到如下内容:

我还是用了nohup把这个命令挂在后台,防止掉线,因为要运行2个小时左右,我加上time命令可以看到运行时间,我用了bwa的index模式来索引参考基因在,具体bwa用法可以自己看文档,但是我们只需要学会索引及比对就好了。有点类似于window下面的软件有一个个菜单栏一样,需要自己的鼠标点击来实现一个个功能,在linux下面就是把命令准备好,然后运行。

4.把fastq文件比对到参考基因组

for i in $(seq 1 6) ;do (nohup ~/biosoft/bwa/bwa-0.7.15/bwa  mem -t 5 -M ~/reference/index/bwa/hg19  KPGP-00001_L${i}_R1.fq.gz KPGP-00001_L${i}_R2.fq.gz 1>KPGP-00001_L${i}.sam 2>KPGP-00001_L${i}.bwa.align.log &);done

这个命令就一句话,但是里面的信息量非常大, 需要熟练掌握linux命令以及shell脚本的语法,但是解析起来也很简单,就是因为我们的fastq文件命名是有规律的,根据规律我构造出一个循环命令,里面的i这个变量会自动扩展成1,2,3,4,5,6依次来用bwa  mem 模式来比对,因为是PE150测序,所以选择这个模式,-M就是选择我们上一步构建好的参考基因组,最后面的 1> 和2>是把软件运行结果输出来,分别是标准输出和标准错误输出,大家可以自行搜索。如果fastq文件的命名发生变化,这个shell脚本是运行不了的,需要临时构建,自己得掌握脚本编写,不然就一个个的比对,手动。

大家可以去看【直播】我的基因组(七):从整体理解全基因组测序数据的变异位点,来了解这个命令的运行结果。

请扫描以下二维码关注我们,获取直播系列的所有帖子!

1

十二 09

【直播】我的基因组(九):拿到数据后要做的事情

时隔好几个月,因为各种各样的原因数据终于拿到了自己的手上,真是不容易啊!

拿到数据后,第一件要做的事情就是检查数据传输的完整性,然后备份!我拿到的数据如下:

可以看到,公司给了我测序仪的下机数据(raw data)和他们质控后的clean data,这个过程减少了6G的数据量,对应着约90亿bp的碱基,相当于减少了3个人的全基因组数据。具体推算公式见前面的系列直播贴!

首先我把数据拷贝到了我上上周买的2T移动硬盘里面,再拷贝到我工作电脑一份,服务器一份,私人电脑一份,另外一个移动硬盘一份。然后删除了公司寄给我的硬盘里面的数据,再把硬盘寄回给公司,然后监督他们删除我所有的数据。(做这么多就是为了保护隐私,当然这个大前提是我已经确定数据没有问题了。)

检查数据传输的完整性就是md5校验,看看数据在拷贝过程中有没有意外的损坏(这个在之前下载数据的时候我也说过)!一般传输数据之前,会用md5命令来生成各个文件的md5值,就是下面的MD5.txt文件里面的内容,然后传输数据之后,需要自行用md5sum -c MD5.txt 来校验文件里面记录的文件的完整性,如果显示都是OK,说明文件拷贝传输过程是没有问题的!但这个过程会耗费大量的磁盘读写,磁盘读写能力是有限的,所以开多个进程并不能加快这一过程。

然后我把公司处理好的bam文件上传到服务器做下游分析,我用的winscp软件把文件传到服务器上的!

从明天起,我们就开始正式对基因组进行分析啦!欢迎围观!

请扫描以下二维码关注我们,获取直播系列的所有帖子!

1

十二 01

GSEA的统计学原理试讲

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

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

Continue reading

十二 01

吐血推荐snpedia数据库,非常丰富的snp信息记录

正好,我拿到了自己的全基因组测序数据,而前些天看到朋友圈推送的文章提到有研究表明STAT4上的rs7574865和HLA-DQ的 rs9275319是国人群中乙型肝炎病毒(HBV)相关肝细胞癌(HCC)遗传易感基因,我就想顺便看看自己在这两个位点的变异情况。一般的流程是先找完变异位点,然后用vep/snpEFF对变异位点进行注释,然后看看有没有这两个位点。但我仅仅是想查看这两个位点,所以我会根据它的rsID来找到它的基因组坐标,再直接call这个位置的变异情况。以前我都是用dnSNP来查看rsID的基因组坐标的,
mkdir -p ~/annotation/variation/human/dbSNP
cd ~/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

Continue reading

十一 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的晚上十一点开讲!
讲义的草稿如下,如果你能看懂草稿,能自己学会,就不用参加本次课程啦。
如果需要问我问题,就赶快找我申请加入交流群,提供本次培训的全套视频和代码!!!

Continue reading

十一 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).

Continue reading

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

Continue reading

十一 25

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

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

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

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

Continue reading

十一 25

kegg在线链接图的颜色设置

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

Continue reading

十一 24

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

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

Continue reading

十一 24

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

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

Continue reading

十一 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中心就可以点击安装!

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

下面的链接选择性观看:

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

十一 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进行画图

Continue reading

十一 23

java环境变量的问题

有篇文章提到了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