三 23

我的优酷视频网盘分享

Posted on by

因为视频教程架构出了点问题，所有不会更新了。现在公布所有下载链接，大家不需要去优酷观看了，太不清晰了，而且还有广告！

我的自频道-优酷视频 

视频百度云盘链接：http://pan.baidu.com/s/1jIvwRD8

下面是以前写的大纲，虽然我不做视频了，但是大家可以按照下面的大纲来自学，希望能对你有所帮助！！！

Continue reading →

三 15

perl程序技巧-检验系统环境或模块安装

Posted on by

这个程序是我在VirusFinder里面发现的，大家可以自行搜索它！

非常好用，建议大家写程序都可以加上这个！

print "\nChecking Java version...\n\n";
my $ret = `java -version 2>&1`;
print "$ret\n";
if (index($ret, '1.6') == -1) {
    printf "Warning: The tool Trinity of the Broad Institute may require Java 1.6.\n\n";
}
print "\nChecking SAMtools...\n\n";
$ret = `which samtools 2>&1`;
if (index($ret, 'no samtools') == -1) {
    printf "%-30s\tOK\n\n", 'SAMtools';
}else{
    printf "%-30s\tnot found\n\n", 'SAMtools';
}
my @required_modules = ("Bio::DB::Sam",
                        "Bio::DB::Sam::Constants",
                        "Bio::SeqIO",
                        "Bio::SearchIO",
                        "Carp",
                        "Config::General",
                        "Cwd",
                        "Data::Dumper",
                        "English",
                        "File::Basename",
                        "File::Copy",
                        "File::Path",
                        "File::Spec",
                        "File::Temp",
                        "FindBin",
                        "Getopt::Std",
                        "Getopt::Long",
                        "IO::Handle",
                        "List::MoreUtils",
                        "Pod::Usage",
                        "threads");
print "\nChecking CPAN modules required by VirusFinder...\n\n";
my $count = 0;
for my $module (@required_modules){
eval("use $module");
if ($@) {
printf "%-30s\tFailed\n", $module;
                $count++;
}
else {
printf "%-30s\tOK\n",     $module;
}
}
if ($count==1){
print "\n\nOne module may not be installed properly.\n\n";
}elsif ($count > 1){
print "\n\n$count modules may not be installed properly.\n\n";
}else{
print "\n\nAll CPAN modules checked!\n\n";
}

 

三 15

perl模块终极解决方案-下

Posted on by

其实可以手动下载local::lib, 这个perl模块，然后自己安装在指定目录，也是能解决模块的问题！

下载之后解压，进入：

 $ perl Makefile.PL --bootstrap=~/.perl  ##这里设置你想把模块放置的目录
 $ make test && make install
 $ echo 'eval $(perl -I$HOME/.perl/lib/perl5 -Mlocal::lib=$HOME/.perl)' >> ~/.bashrc

等待几个小时即可！！！

添加好环境变量之后，就可以用

perl -MCPAN -Mlocal::lib -e 'CPAN::install(LWP)'
这样的模式下载模块了，所有的模块都会存储在$HOME/.perl/lib/perl5 里面！！！
如果是新写的perl程序，需要在开头加入 use local::lib;   
# sets up a local lib at ~/perl5才能使用该模块！  非常重要，非常重要，非常重要！！！

 

其实你也可以直接打开 ~/.bashrc，然后写入下面的内容

PERL5LIB=$PERL5LIB:/PATH_WHERE_YOU_PUT_THE_PACKAGE/source/bin/perl_module; export PERL5LIB
可以把perl模块安装在任何地方，然后通过这种方式去把模块加载到你的perl程序！

 

三 15

perl模块终极解决方案-上

Posted on by

不管别人怎么说，反正我是非常喜欢perl语言的！

也会继续学习，以前写过不少perl模块的博客，发现有点乱，正好最近看到了关于local::lib这个模块。

居然是用来解决没有root权限的用户安装，perl模块问题的！

首先说一下，如果是root用户，模块其实没有问题，直接用cpan下载器，几乎能解决所有的模块下载安装问题！

但是如果是非root用户，那么就麻烦了，很难用自动的cpan下载器，这样只能下载模块源码，然后编译，但是编译有个问题，很多模块居然是依赖于其它模块的，你的不停地下载其它依赖模块，最后才能解决，特别麻烦

但是，只需要运行下面的代码：

wget -O- http://cpanmin.us | perl - -l ~/perl5 App::cpanminus local::lib
eval `perl -I ~/perl5/lib/perl5 -Mlocal::lib`
echo 'eval `perl -I ~/perl5/lib/perl5 -Mlocal::lib`' >> ~/.profile
echo 'export MANPATH=$HOME/perl5/man:$MANPATH' >> ~/.profile

就能拥有一个私人的cpan下载器，~/.profile可能需要更改为.bash_profile, .bashrc, etc等等，取决于你的linux系统！

然后你直接运行cpanm Module::Name，就跟root用户一样的可以下载模块啦！

或者用下面的方式在shell里面安装模块，其中ext是模块的安装目录，可以修改

 perl -MTime::HiRes -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext Time::HiRes;
perl -MFile::Path -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext File::Path;
perl -MFile::Basename -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext File::Basename;
perl -MFile::Copy -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext File::Copy;
perl -MIO::Handle -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext IO::Handle;
perl -MYAML::XS -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext YAML::XS;
perl -MYAML -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext YAML;
perl -MXML::Simple -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext XML::Simple;
perl -MStorable -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext Storable;
perl -MStatistics::Descriptive -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext Statistics::Descriptive;
perl -MTie::IxHash -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext Tie::IxHash;
perl -MAlgorithm::Combinatorics -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext Algorithm::Combinatorics;
perl -MDevel::Size -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext Devel::Size;
perl -MSort::Key::Radix -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext Sort::Key::Radix;
perl -MSort::Key -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext Sort::Key;
perl -MBit::Vector -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext Bit::Vector;
perl -M"feature 'switch'" -e 1 > /dev/null 2>&1 || cpanm -v --notest -l ext feature;

下面是解释为什么这样可以解决问题！！！

 

What follows is a brief explanation of what the commands above do.

wget -O- http://cpanmin.us fetches the latest version of cpanm and prints it to STDOUT which is then piped to perl - -l ~/perl5 App::cpanminus local::lib. The first - tells perl to expect the program to come in on STDIN, this makes perl run the version of cpanm we just downloaded.perl passes the rest of the arguments to cpanm. The -l ~/perl5 argument tells cpanm where to install Perl modules, and the other two arguments are two modules to install. [App::cpanmins]1 is the package that installs cpanm. local::lib is a helper module that manages the environment variables needed to run modules in local directory.

After those modules are installed we run

eval `perl -I ~/perl5/lib/perl5 -Mlocal::lib`

to set the environment variables needed to use the local modules and then

echo 'eval `perl -I ~/perl5/lib/perl5 -Mlocal::lib`' >> ~/.profile

to ensure we will be able to use them the next time we log in.

echo 'export MANPATH=$HOME/perl5/man:$MANPATH' >> ~/.profile

will hopefully cause man to find the man pages for your local modules.

 

 

 

 

 

 

 

 

三 15

基因组各种版本对应关系

Posted on by 
我是受到了SOAPfuse的启发才想到整理各种基因组版本的对应关系，完整版！！！
以后再也不用担心各种基因组版本混乱了，我还特意把所有的下载链接都找到了，可以下载任意版本基因组的基因fasta文件，gtf注释文件等等！！！
首先是NCBI对应UCSC，对应ENSEMBL数据库：
GRCh36 (hg18): ENSEMBL release_52.
GRCh37 (hg19): ENSEMBL release_59/61/64/68/69/75.
GRCh38 (hg38): ENSEMBL  release_76/77/78/80/81/82.
可以看到ENSEMBL的版本特别复杂！！！很容易搞混！
但是UCSC的版本就简单了，就hg18,19,38, 常用的是hg19，但是我推荐大家都转为hg38
看起来NCBI也是很简单，就GRCh36,37,38，但是里面水也很深！
Feb 13 2014 00:00    Directory April_14_2003
Apr 06 2006 00:00    Directory BUILD.33
Apr 06 2006 00:00    Directory BUILD.34.1
Apr 06 2006 00:00    Directory BUILD.34.2
Apr 06 2006 00:00    Directory BUILD.34.3
Apr 06 2006 00:00    Directory BUILD.35.1
Aug 03 2009 00:00    Directory BUILD.36.1
Aug 03 2009 00:00    Directory BUILD.36.2
Sep 04 2012 00:00    Directory BUILD.36.3
Jun 30 2011 00:00    Directory BUILD.37.1
Sep 07 2011 00:00    Directory BUILD.37.2
Dec 12 2012 00:00    Directory BUILD.37.3
可以看到，有37.1,   37.2，  37.3 等等，不过这种版本一般指的是注释在更新，基因组序列一般不会更新！！！
反正你记住hg19基因组大小是3G，压缩后八九百兆即可！！！
如果要下载GTF注释文件，基因组版本尤为重要！！！
对于ensembl：
变幻中间的release就可以拿到所有版本信息：ftp://ftp.ensembl.org/pub/
对于UCSC，那就有点麻烦了：
需要选择一系列参数：
2. Select the following options:
clade: Mammal
genome: Human
assembly: Feb. 2009 (GRCh37/hg19)
group: Genes and Gene Predictions
track: UCSC Genes
table: knownGene
region: Select "genome" for the entire genome.
output format: GTF - gene transfer format
output file: enter a file name to save your results to a file, or leave blank to display results in the browser
3. Click 'get output'.
 现在重点来了，搞清楚版本关系了，就要下载呀！
UCSC里面下载非常方便，只需要根据基因组简称来拼接url即可：
或者用shell脚本指定下载的染色体号：
for i in $(seq 1 22) X Y M;
do echo $i;
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
## 这里也可以用NCBI的：ftp://ftp.ncbi.nih.gov/genomes/M_musculus/ARCHIVE/MGSCv3_Release3/Assembled_Chromosomes/chr前缀
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg19.fasta;
done
rm -fr chr*.fasta
三 15

融合基因检测软件-soapfusion

Posted on by

开发单位：华大，SOAP系列软件套装！

功能：检测合基因
优点：在现有的各种软件里面表现算是最好的
算法：是hash index，跟其它bwt算法不太一样
其它软件有： FusionSeq [21], deFuse [22], TopHat-Fusion [23], FusionHunter [24], SnowShoes-FTD [25], chimerascan [26] and FusionMap [27]
具体的算法我没看，因为只是有需求，正好有一些RNA-seq数据又想看看样本融合基因情况。所以就测试这个软件，通俗点说，融合基因原理其实很简单，如果有足够多的reads一部分比对到一个基因，另一部分比对到另一个基因，就可以说明它们两个基因发生了融合现象！如果是PE测序，那么更方便，左右两端reads比对情况也可以考虑。我就不多说废话了，直接上教程吧！
一，软件安装
下载压缩包，解压后即可使用！！！
推荐用最新版，然后看作者说明书的时候也要看清楚！
我反正好几次都搞糊涂了，最后联系了作者才搞明白，作者说他想更新到2.0版本，直接用HISAT的比对sam文件来做，但是还在筹备中，我觉得有点悬！
1
解压后是一堆perl程序，都在source目录下，source目录下面还有bin下面附带了几个第三方软件，包括bwa，blast和soap，最后都用得着！
有个很重要的问题，一定要软件自带的perl模块添加到perl的环境变量。不然那些perl程序运行会报错！
配置文件需要修改，就把几个目录放进去即可

二，输入数据准备

这里最重要的就是制作数据库！！！
作者给了非常详细的制作过程，我觉得还是不够清楚，所以再讲一遍！
首先下载5个文件：
6.5K Jun 15  2009 cytoBand.txt.gz
3.0G Oct 12  2012 hg19.fa
2.5M Mar 15 10:30 HGNC_Gene_Family_dataset
38M Feb  8  2014 Homo_sapiens.GRCh37.75.gtf.gz
202 Jan 19 16:07 HumanRef_refseg_symbols_relationship.list

文件下载地址，作者已经给出了！

我把这些文件都放在的当前文件夹下面的raw这个子文件夹，因为我要当前文件夹作为该软件的database文件夹！！！
然后运行命令！
我在SOAPfuse-v1.27文件下面运行：
perl ../SOAPfuse-v1.27/source/SOAPfuse-S00-Generate_SOAPfuse_database.pl  \
-wg raw/hg19.fa  -gtf raw/Homo_sapiens.GRCh37.75.gtf.gz  -cbd raw/cytoBand.txt.gz   -gf raw/HGNC_Gene_Family_dataset \
-rft raw/HumanRef_refseg_symbols_relationship.list \
 -sd ../SOAPfuse-v1.27 -dd ./

这一步耗时很长，4~6小时，创造了transcript.fa和gene.fa，然后还对他们建立bwa和soap的index，所以有点慢！

构建成功会有提示：
Congratulations!
You have constructed SOAPfuse database files successfully.
These database files are all stored in directory you supplied:
/home/jmzeng/biosoft/SOAPfuse/db_v1.27/
They are all generated based on public data files you supplied:
whole_genome_fasta_file:   /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/hg19.fa
gtf_annotation_file:       /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/Homo_sapiens.GRCh37.75.gtf.gz
Chr_Bandregion_file:       /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/cytoBand.txt.gz
HGNC_gene_family_file:     /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/HGNC_Gene_Family_dataset
gtf_segname2refseg_list:   /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/HumanRef_refseg_symbols_relationship.list
这些目录很重要，接下来制作配置文件会用得着！
To use these database files, just set the 'DB_db_dir' in config file as belowed:
DB_db_dir  =   /home/jmzeng/biosoft/SOAPfuse/db_v1.27
配置文件需要修改下面5个
DB_db_dir = /DATABASE_DIR/
PG_pg_dir = /TOOL_DIR/source/bin
PS_ps_dir = /TOOL_DIR/source
PD_all_out = /out_directory/
PA_all_fq_postfix = PostFix
其实你仔细阅读了说明书，你就知道该修改成什么样子了！
最后制作sample list文件
我这里只有一个sample,所以文件就一句话即可
test test test 100
所以我的有下面两个文件，都是为了顺应作者的需求我才搞了test/test/test这么无聊的东西！！！
/home/jmzeng/test_for_soapfuse/test/test/test_1.fq.gz
/home/jmzeng/test_for_soapfuse/test/test/test_2.fq.gz
如果你有多个sample需要一起运行，你就要仔细读作者的readme了，它把这个配置文件搞得特别复杂！！！

三，运行命令

如果文件都准备好了，运行命令非常简单！！
perl SOAPfuse-RUN.pl -c <config_file> -fd <WHOLE_SEQ-DATA_DIR> -l <sample_list> -o <out_directory> [Options]

运行的非常慢！！！

因为需要重新比对，知道

四，数据结果解读

结果，作者已经说的很清楚了，我就不多说了！
三 15

使用virusSeq对NGS数据分析病毒整合位点

Posted on by
开发单位：安德森癌症研究所
功能：对NGS数据进行分析，探测已知病毒在人基因组整合情况
从fastq文件开始，需要借助MOSAIK进行比对
主程序就两个perl程序，不需要安装！

Continue reading →

三 15

新的比对工具MOSAIK

Posted on by
功能：序列比对，类似于BWA，Bowtie
优点：全平台，甚至支持pacbio的三代测序长reads
算法：是hash index，跟其它bwt算法不太一样
作者：WP Lee - ‎2014 - ‎被引用次数：70 - ‎相关文章

Continue reading →

三 12

bioconductor中文社区招募站长

Posted on by 
我已经构建好bioconductor中文社区的雏形，大家可以进去看看！

写在前面

突然发现我的bioconductor.cn这个域名都快要过期了！

哈哈，才想起一年前的计划到现在还没开始实施，实在不像我的风格，可能是水平到了一定程度吧，很多初级工作不像以前那样事无巨细的把关了。正好，借这个机会找几个朋友帮我一起完成这个bioconductor中文社区计划！

 

Continue reading →

三 10

THis is a test for markdown

Posted on by

A wonderful tool to write blog !

Hi,all. I just started to use github and found this magic way to write a blog.
so I do a simpler test here.
Don't be confused, there is nothing new for bioinformatic here.

And it very easy to use Markdown in wordpress,just download a plugin.
Below is some useful links to help you familar with markdown as soon as possible.

  • http://bhttp://mahua.jser.me/log.csdn.net/kaitiren/article/details/38513715
  • https://github.com/guodongxiaren/README/blob/master/README.md
  • http://wowubuntu.com/markdown/
  • http://mahua.jser.me/

Also I want to recommend some useful web-editor for markdown.

  • http://mahua.jser.me/
  • https://www.madoko.net/
  • https://www.zybuluo.com/mdeditor

It's easy to use Markdown, but to get the most out of it, you still need to understand it and keep practising

三 03

用R语言包从EBI的arrayexpress数据库里面下载芯片数据

Posted on by
这个包跟GEOquery区别不是很大，只不过一个是正对NCBI的GEO数据库，一个是针对EBI的arrayexpress数据库，只有对写自动化脚本的人来说才有需求，一般个人分析者都是自己去数据库主页里面查找，然后拿到下载链接，一个个下载。
从EBI的arrayexpress数据库里面下载芯片数据：
update to 2016-3-1 11:41:27
63890 experiments
1912744 assays
40.53 TB of archived data 数据量还是蛮大的
所有的data，都可以在ftp服务器里面下载：ftp://ftp.ebi.ac.uk/pub/databases/arrayexpress/data/experiment/BUGS/
根据ID号很整齐的储存着。
也可以用一个R语言包：ArrayExpress R package
2009年，那个时候R语言用的人很少，这个简单的包都可以发文章，现在看来简直不可思议！
其实大部分数据都是跟GEO数据库对应的：比如https://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-55645/  对应于：GEO - GSE55645
比如对NASH表达数据查找：https://www.ebi.ac.uk/arrayexpress/search.html?query=NASH++expression  30条结果里面只有4条是arrayexpress数据库独有的！
biocLite("ArrayExpress")
library(ArrayExpress)
1
如果用R语言，搜索如下：
可以用sets = queryAE(keywords = "NASH+expression", species = "homo+sapiens")
2
效果是一样的！
下载数据用：
back = getAE("E-MEXP-3291")
下载其实也就是里面存储了链接，直接调用R语言的下载函数即可！
3
一般没必要下载原始测序文件，直接用下面这个函数就可以得到一个数据对象，可以直接得到表达矩阵和实验的metadata

rawset = ArrayExpress("E-MEXP-3291")

二 28

用R的shiny包写一个基因的ID转换小程序

Posted on by
基因的ID现在可能有一千多种，但是我们一般只用NCBI的entrez ID，和HUGO组装规定的gene symbol和gene name
再者有人会用ensembl的ID，至于UCSC的ID很少人用了，更别说剩余的一大堆稀奇古怪的ID了。
ID_example
其实没必要自己写一个软件来转换，因为非常多的现成软件都可以做到!
不过可以拿这个简单的任务来练手！
首先自己做好数据库文件：
我是把数据写到了一个sqlite文件里面了！结构很简单
然后设计软件的界面，俗称UI端，也是很容易，在shiny里面，记住几个控件的函数即可！
UI_1 UI_2
具体代码见我的github代码：https://github.com/jmzeng1314/my-R/tree/master/4-ID-convert-using-shiny
二 28

在R里面操作SQLite

Posted on by
我前面写到过如何把数据库写到mysql，但是发现其实msyql并不方便，需要连接数据库什么的，如果发布一个离线小网页，这时候sqlite的优点就显示出来了！
基础代码很简单：
library(RSQLite)
sqlite    <- dbDriver("SQLite")
con <- dbConnect(sqlite,"hg19_bioconductor.sqlite") # makes a new file
suppressMessages(library(org.Hs.eg.db))
kegg2ID=toTable(org.Hs.egPATH)
#[1] "gene_id" "path_id"
dbWriteTable(con,'keggID2geneID',kegg2ID,row.name=F,overwrite=T)
具体代码，可以看我的github主页：https://github.com/jmzeng1314/my-R/blob/master/3-get-hg19-gene-mapping/get-hg19-gene-mapping-bioconductor.R
做出来的数据，如下，就是几个table存储在文件里面！
 2
最后这些数据都保存在了当前工作目录下的hg19_bioconductor.sqlite文件里面！
在其它程序里面就可以直接调用这个文件，而不需要加载一大堆的包了！
1
library(KEGG.db)
library(GO.db)

library(org.Hs.eg.db )

 

二 28

把bioconductor的gene mapping信息上传到mysql数据库

Posted on by
在R语言里面的bioconductor系列包里面有一个物种注释信息包，其中人类是org.Hs.eg.db
里面有基于人的hg19基因组版本的大部分基因信息之间的转换数据！
包括基因的entrez ID，symbol，name，locus，refseq，kegg pathway，GO ontology对应关系！
但是它们散布在该包的各个数据结构里面！
> ls("package:org.Hs.eg.db")
 [1] "org.Hs.eg"                "org.Hs.eg.db"            
 [3] "org.Hs.eg_dbconn"         "org.Hs.eg_dbfile"        
 [5] "org.Hs.eg_dbInfo"         "org.Hs.eg_dbschema"      
 [7] "org.Hs.egACCNUM"          "org.Hs.egACCNUM2EG"      
 [9] "org.Hs.egALIAS2EG"        "org.Hs.egCHR"            
[11] "org.Hs.egCHRLENGTHS"      "org.Hs.egCHRLOC"         
[13] "org.Hs.egCHRLOCEND"       "org.Hs.egENSEMBL"        
[15] "org.Hs.egENSEMBL2EG"      "org.Hs.egENSEMBLPROT"    
[17] "org.Hs.egENSEMBLPROT2EG"  "org.Hs.egENSEMBLTRANS"   
[19] "org.Hs.egENSEMBLTRANS2EG" "org.Hs.egENZYME"         
[21] "org.Hs.egENZYME2EG"       "org.Hs.egGENENAME"       
[23] "org.Hs.egGO"              "org.Hs.egGO2ALLEGS"      
[25] "org.Hs.egGO2EG"           "org.Hs.egMAP"            
[27] "org.Hs.egMAP2EG"          "org.Hs.egMAPCOUNTS"      
[29] "org.Hs.egOMIM"            "org.Hs.egOMIM2EG"        
[31] "org.Hs.egORGANISM"        "org.Hs.egPATH"           
[33] "org.Hs.egPATH2EG"         "org.Hs.egPFAM"           
[35] "org.Hs.egPMID"            "org.Hs.egPMID2EG"        
[37] "org.Hs.egPROSITE"         "org.Hs.egREFSEQ"         
[39] "org.Hs.egREFSEQ2EG"       "org.Hs.egSYMBOL"         
[41] "org.Hs.egSYMBOL2EG"       "org.Hs.egUCSCKG"         
[43] "org.Hs.egUNIGENE"         "org.Hs.egUNIGENE2EG"     
[45] "org.Hs.egUNIPROT"        
> 
其实比较重要的信息，我们需要把它们统一关联起来，做成一个表格，这样可以上传到数据库里面，做成网页，可视化展现，查找！
suppressMessages(library(org.Hs.eg.db))
all_EG=mappedkeys(org.Hs.egSYMBOL)
###下面的表格最后是基于entrez ID而关联起来的！
tmp=unlist(as.list(org.Hs.egSYMBOL))
EG2Symbol=data.frame(EGID=names(tmp),symbol=as.character(tmp))
 
tmp=unlist(as.list(org.Hs.egENSEMBL))
EG2ENSEMBL=data.frame(EGID=names(tmp),ENSEMBL=as.character(tmp))
 
tmp=unlist(as.list(org.Hs.egGENENAME))
EG2name=data.frame(EGID=names(tmp),name=as.character(tmp))
 
tmp=unlist(as.list(org.Hs.egMAP))
EG2MAP=data.frame(EGID=names(tmp),MAP=as.character(tmp))
     
###EG2GO and    using mySQL
EG2path=as.list(org.Hs.egPATH)
EG2path=lapply(EG2path, function(x) paste(x,collapse = ":"))
tmp=unlist(EG2path)
EG2path=data.frame(EGID=names(tmp),path=as.character(tmp))
 
#as.list(head(org.Hs.egGO2ALLEGS))
GO2allEG=as.list(org.Hs.egGO2ALLEGS)
tmp=unlist(GO2allEG)
names(tmp)=substring(names(tmp),1,10) ## change GO:0000002.IMP to GO:0000002 
EG2allGO <- tapply(tmp,tmp,function(x){names(x)})
EG2allGO=lapply(EG2allGO, function(x) paste(x,collapse = ":"))
tmp=unlist(EG2allGO)
EG2GO=data.frame(EGID=names(tmp),GO=as.character(tmp))
 
tmp=merge(EG2Symbol,EG2MAP,by='EGID',all=TRUE)
tmp=merge(tmp,EG2ENSEMBL,by='EGID',all=TRUE)
tmp=merge(tmp,EG2path,by='EGID',all=TRUE)
tmp=merge(tmp,EG2name,by='EGID',all=TRUE)
my_gene_mapping=merge(tmp,EG2GO,by='EGID',all=TRUE)
##[1] "EGID"    "symbol"  "MAP"     "ENSEMBL" "path"    "name"    "GO"
1
因为我用的merge的all=TRUE，所以最后记录会越来越多
最后的结果就是下面这样，各种ID转换都储存到了my_gene_mapping这个数据对象里面！
可以看到有些基因是没有ensemble数据库对应的ID的，非常多的基因没有被注释到KEGG通路数据库！
我这里如果一个基因对应多个通路，我用冒号连接起来了，成了一个字符串！在后面会有用！

2

现在需要把它上传到数据库里面，这样我就可以用R的shiny可视化网页来查找数据库，做ID转换啦！！！
suppressMessages(library(RMySQL))
con <- dbConnect(MySQL(), host="127.0.0.1", port=3306, user="root", password="11111111")
dbSendQuery(con, "USE test")
dbWriteTable(con,'my_gene_mapping',my_gene_mapping)
dbDisconnect(con)
然后就可以在自己的mysql里面看到它啦！！！ 
其实 "org.Hs.egOMIM"  也很重要，不过我这里只是举个例子，就不深究了！
还有一点，就是那个pathway ID，因为我要以entrez ID为主键，所以把多个pathway给合并了！
suppressMessages(library(org.Hs.eg.db))
kegg2ID=toTable(org.Hs.egPATH)
#[1] "gene_id" "path_id"
dbWriteTable(con,'kegg2ID_bioconductor',kegg2ID,row.name=F,overwrite=T)
go2id=toTable(org.Hs.egGO2ALLEGS)
## gene_id      go_id Evidence Ontology
dbWriteTable(con,'go2id_bioconductor',go2id,row.name=F,overwrite=T)
dbDisconnect(con)
用这个代码，可以在数据库里面多创建两个表，会有用的！
二 28

用crossmap代替liftover做基因组坐标转换

Posted on by
其实国际三大主流生物信息学数据库运营单位都出了自己的基因组坐标转换，它们分别是 (UCSC liftOver, NCBI Remap, Ensembl API)
Ensembl's Assembly Converter.是基于crossmap的，我觉得挺好用的，就介绍给大家！！！

This online tool currently uses CrossMap, which supports a limited number of formats (see our online documentation for details of the individual data formats listed below). CrossMap also discards metadata in files, so track definitions, etc, will be lost on conversion.

Important note: CrossMap converts WIG files to BedGraph internally for efficiency, and also outputs them in BedGraph format.

但是不知道为什么UCSC的liftover最出名，我也写过它的教程，（http://www.bio-info-trainee.com/?p=990）

Continue reading →

二 22

没必要学shell进阶语法

Posted on by

因为大部分生物信息学软件都是linux版本的，所以生物信息学数据分析工作者必备技能就是linux，但是大部分人只是拿他当个中转站，我以前也是，直到接触了大批量的任务，自动化流程，才明白这里面的水太深了，不过无所谓，凭我个人的观点，其实shell的进阶语法真的不必要！

当然，只是我一家之言！
我实在是不想去背诵大括号，小括号，中括号以及双重括号到底区别是什么！
http://www.bio-info-trainee.com/?p=1018  [],[[]],(),(()),{},{{}},以及在前面加上$的区别，以及它们互相杂交组合的区别！！！

我也不想去搞明白操作符两边是否加空格的区别是什么了。

if((i%5==0)) 来判断变量是否被一个数整除
i=$((i+1))来表示变量自增。
这些东西真的很诡异！
如果你有qsub，condor等任务提交系统，那么你只需要熟悉他们就可以了，但大部分散兵游勇的生物信息学家并没有集群，所以压根不会接触任务提交系统，就需要些自动化脚本了！
受限制与机器的cpu以及内存数，需要判断提交了多少任务，等待多久再执行，所以会把一个简单的自动化脚本写的很复杂！
比如下面这个脚本：cat >download_hg38_from_UCSC.sh

for i in $(seq 1 22) X Y M;
do echo $i;
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg38.fa;
done
rm -fr chr*.fa
可以下载hg38基因组的fasta文件，但是是分染色体一个个下载的！
再比如下面这个，批量做GSEA分析的脚本：
while read id
do
echo $id
gene=`echo $id |awk '{print $1}'`
probe=`echo $id |awk '{print $2}'`
echo $i
do_GSEA $probe $gene; ##这里是我自己定义的一个function，就不贴出来了
if((i%5==0))
then
sleep 10  ##重点就在这里，每次提交的任务有限制，所以需要休息，不然机器的cpu负载太高！
fi
i=$((i+1))
done <$1

如果，还有其它功能需要实现，我们可以把脚本写的更负载，纯粹的用shell，需要搜索更多的shell技巧。

但是事实上并没有这个必要，我们现在有了更方便的脚本语言，比如我所擅长的perl
我写一个nohup提交任务的脚本！
## perl nohup.pl   deep_count.sh  0
## perl nohup.pl   deep_count.sh  1
## perl nohup.pl   deep_count.sh  2

[perl]
## perl nohup.pl   deep_count.sh  0
## perl nohup.pl   deep_count.sh  1
## perl nohup.pl   deep_count.sh  2
$i=1;
open FH,$ARGV[0];
while(<FH>){
    chomp;
    next unless $.%3==$ARGV[1];
    $cmd="nohup  $_  &";
    print "$cmd\n";
    system($cmd);
    sleep(10800) if $i%5==4;
    $i++;
    #exit;
}
[/perl]

我尝试过用shell，写了很久，总是报错，但是用perl，一分钟我就写完了，所以，最好是用自己熟悉的一种语法最好！
二 21

画基因结构图！

Posted on by
一个基因有一个甚至多个转录本！根据转录及翻译现象可以把一个基因人为的定义成多种结构！
首先自己查资料搞明白转录本，外显子和内含子，5端UTR区域和CDS区域，还有3端UTR区域的概念。
真核生物的基因含有外显子和内含子，是区别原核生物的特征之一，能转录的就是外显子，不能转录的就是内含子！
内含子（英语：Intron）是一个基因中非编码DNA片段，它分开相邻的外显子。更精确的定义是：内含子是阻断基因线性表达的序列。DNA上的内含子会被转录到前体RNA中，但RNA上的内含子会在RNA离开细胞核进行转译前被剪除。在成熟mRNA被保留下来的基因部分被称为外显子。
通常我们看到基因结构示意图会像下面这样！
gene-structure
但是对于这个图，我之前有非常多的疑问，直到我自己写程序去仔细统计，画图探究后才真正搞明白！
从图中可以看到5端UTR区域后面接着的是第一外显子，但是事实上不是！
外显子和内含子是转录本上面的概念，是基于转录这个行为的定义。
而5端UTR区域和CDS区域，还有3端UTR区域,是基于翻译这种行为的定义！
把它们画成这样，是严重的干扰信息。
实际上，我们需要更正很多概念，我检验了一下，得到下面几个正确的：
1，如果基因有多个转录本，基因的起始坐标，就是该基因所有转录本的第一个外显子的起始坐标的最小值，同理基因的终止坐标，就是该基因的所有转录本的最后一个外显子的终止坐标的最大值。
2，通过这个概念，可以把基因分成闭合基因和非闭合基因。 闭合基因:有一个最长转录本使得基因起始终止坐标等于该最长转录本的起始终止坐标。（这个是我乱说的，并没有这个定义）
3，如果基因只有一个转录本，那么基因的起始终止坐标，就是转录本的起始终止坐标！
4，一个基因的一个转录本的5'utr区域可以包括多个外显子区域，前者是翻译行为，后者是转录行为
5，起始密码子和终止密码子是CDS的起止处，是基于翻译的概念
6，一个基因的多个转录本的外显子坐标不一定会排列整齐，每个转录本的剪切位点并不一定要比其它转录本一致！
比如对这个ANXA1基因来说，非常多的转录本，但是基因的起始终止坐标，是所有转录本起始终止坐标的极大值和极小值！同时，它是一个闭合基因，因为它存在一个转录本的起始终止坐标等于该基因的起始终止坐标。可以看到它的外显子并不是非常整齐的，虽然多个转录本会共享某些外显子，但是也存在某些外显子比同区域其它外显子长的现象！
 anxa1-gene-structure
再比如下面这个例子：对DDX11L11这个基因来说，前两个外显子都不会翻译，直到第三个外显子才开始翻译，构成CDS。
有些转录本是没有utr的，所以该转录本的起始坐标，就是CDS的起始坐标
5UTR包括多个exon行为
首先把gtf文件格式化导入mysql数据库
我用的是broadinstitute提供的GENCODE的gtf文件，是hg19版本的
Modified GENCODE GTF file for human with contig names of form ("1","2", etc)
##description: evidence-based annotation of the human genome (GRCh37), version 7 (Ensembl 62)
##provider: GENCODE
##format: gtf
##date: 2011-03-23
很简单，下载数据，parse好之后，用R导入mysql即可！
下面是mysql代码：
--mysql code
use test;
drop table  if exists hg19_gtf;
create table hg19_gtf (
    gene_name VARCHAR(30),
    transcript_name VARCHAR(30) ,
    record  VARCHAR(15) NOT NULL ,
    chr VARCHAR(2) NOT NULL ,
    start INT NOT NULL ,
    end INT NOT NULL ,
    source VARCHAR(10) NOT NULL ,
    strand VARCHAR(1) NOT NULL ,
    gene_id VARCHAR(30) NOT NULL ,
    transcript_id VARCHAR(30) NOT NULL ,
    gene_status VARCHAR(30) ,
    gene_type VARCHAR(30)  ,
    transcript_type VARCHAR(30) ,
    transcript_status VARCHAR(30)
);
--我的网页好像不支持mysql代码高亮，大家凑合着看吧，反正就是一个简单的建表语句！
select * from hg19_gtf limit 100;
select * from hg19_gtf where gene_name='DMD';
select count(*) from hg19_gtf where gene_name='DMD' and record='start_codon';  --18 start condon
select count(distinct(transcript_name)) from hg19_gtf where gene_name='DMD' ;  --34 transcript
select count(distinct(transcript_name)) c ,gene_name from hg19_gtf where record='transcript' group by gene_name  order by c desc;
我是随意设计的一个表，主要是为了画图方便！
接下来是R code来具体的对基因进行画图！

[perl]
suppressMessages(library(ggplot2))
suppressMessages(library(RMySQL))
con <- dbConnect(MySQL(), host="127.0.0.1", port=3306, user="root", password="11111111")
dbSendQuery(con, "USE test")
gene='SOX10'
#gene='DDX11L11'
if (T){
query=paste("select * from hg19_gtf where gene_type='protein_coding' and gene_name=",shQuote(gene),sep="")
structure=dbGetQuery(con,query)
tmp_min=min(c(structure$start,structure$end))
structure$new_start=structure$start-tmp_min
structure$new_end=structure$end-tmp_min
tmp_max=max(c(structure$new_start,structure$new_end))
num_transcripts=nrow(structure[structure$record=='transcript',])
tmp_color=rainbow(num_transcripts)
x=1:tmp_max;y=rep(num_transcripts,length(x))
#x=10000:17000;y=rep(num_transcripts,length(x))
plot(x,y,type = 'n',xlab='',ylab = '',ylim = c(0,num_transcripts+1))
title(main = gene,sub = paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
j=0;
tmp_legend=c()
for (i in 1:nrow(structure)){
tmp=structure[i,]
if(tmp$record == 'transcript'){
j=j+1
tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
}
if(tmp$record == 'exon') lines(c(tmp$new_start,tmp$new_end),c(j,j),col=tmp_color[j],lwd=4)
}
# legend('topleft',legend=tmp_legend,lty=1,lwd = 4,col = tmp_color);

}
[/perl]

通过这个代码，，就能简单的得到我上面显示的基因结构图啦！

 

二 18

生物信息学分析过程中常见文件格式

Posted on by

刚开始接触生物信息学的时候我也很纠结什么fastq,fastq,sam,bam,vcf,maf,gtf,bed,psl等等，甚至还有过时了的NCBI,ENSEMBL格式，如果是我刚开始 学的时候，我倒是很愿意把他们全部搞透彻，写详细的说明书，但是现在成长了，这些东西感觉很low了，正好我看到了一篇帖子讲数据格式的收集大全，分享给大家，希望初学者能多花点时间好好钻研！

https://www.biostars.org/p/55351/

每种文件格式的定义，都是有它的道理的，大部分是因为一个比较流行的软件，少量的数据格式是因为国际组织广泛认可而流行的
二 15

用R获取芯片探针与基因的对应关系三部曲-bioconductor

Posted on by

现有的基因芯片种类不要太多了！

但是重要而且常用的芯片并不多！
一般分析芯片数据都需要把探针的ID切换成基因的ID，我一般喜欢用基因的entrez ID。
一般有三种方法可以得到芯片探针与gene的对应关系。
金标准当然是去基因芯片的厂商的官网直接去下载啦！！！
一种是直接用bioconductor的包

一种是从NCBI里面下载文件来解析好！
首先，我们说官网，肯定可以找到，不然这种芯片出来就没有意义了！
然后，我们看看NCBI下载的，会比较大
这两种方法都比较麻烦，需要一个个的来！
所以我接下来要讲的是用R的bioconductor包来批量得到芯片探针与gene的对应关系！
一般重要的芯片在R的bioconductor里面都是有包的，用一个R包可以批量获取有注释信息的芯片平台，我选取了常见的物种，如下：
        gpl           organism                  bioc_package
1     GPL32       Mus musculus                        mgu74a
2     GPL33       Mus musculus                        mgu74b
3     GPL34       Mus musculus                        mgu74c
6     GPL74       Homo sapiens                        hcg110
7     GPL75       Mus musculus                     mu11ksuba
8     GPL76       Mus musculus                     mu11ksubb
9     GPL77       Mus musculus                     mu19ksuba
10    GPL78       Mus musculus                     mu19ksubb
11    GPL79       Mus musculus                     mu19ksubc
12    GPL80       Homo sapiens                        hu6800
13    GPL81       Mus musculus                      mgu74av2
14    GPL82       Mus musculus                      mgu74bv2
15    GPL83       Mus musculus                      mgu74cv2
16    GPL85  Rattus norvegicus                        rgu34a
17    GPL86  Rattus norvegicus                        rgu34b
18    GPL87  Rattus norvegicus                        rgu34c
19    GPL88  Rattus norvegicus                         rnu34
20    GPL89  Rattus norvegicus                         rtu34
22    GPL91       Homo sapiens                      hgu95av2
23    GPL92       Homo sapiens                        hgu95b
24    GPL93       Homo sapiens                        hgu95c
25    GPL94       Homo sapiens                        hgu95d
26    GPL95       Homo sapiens                        hgu95e
27    GPL96       Homo sapiens                       hgu133a
28    GPL97       Homo sapiens                       hgu133b
29    GPL98       Homo sapiens                     hu35ksuba
30    GPL99       Homo sapiens                     hu35ksubb
31   GPL100       Homo sapiens                     hu35ksubc
32   GPL101       Homo sapiens                     hu35ksubd
36   GPL201       Homo sapiens                       hgfocus
37   GPL339       Mus musculus                       moe430a
38   GPL340       Mus musculus                     mouse4302
39   GPL341  Rattus norvegicus                       rae230a
40   GPL342  Rattus norvegicus                       rae230b
41   GPL570       Homo sapiens                   hgu133plus2
42   GPL571       Homo sapiens                      hgu133a2
43   GPL886       Homo sapiens                     hgug4111a
44   GPL887       Homo sapiens                     hgug4110b
45  GPL1261       Mus musculus                    mouse430a2
49  GPL1352       Homo sapiens                       u133x3p
50  GPL1355  Rattus norvegicus                       rat2302
51  GPL1708       Homo sapiens                     hgug4112a
54  GPL2891       Homo sapiens                       h20kcod
55  GPL2898  Rattus norvegicus                     adme16cod
60  GPL3921       Homo sapiens                     hthgu133a
63  GPL4191       Homo sapiens                       h10kcod
64  GPL5689       Homo sapiens                     hgug4100a
65  GPL6097       Homo sapiens               illuminaHumanv1
66  GPL6102       Homo sapiens               illuminaHumanv2
67  GPL6244       Homo sapiens   hugene10sttranscriptcluster
68  GPL6947       Homo sapiens               illuminaHumanv3
69  GPL8300       Homo sapiens                      hgu95av2
70  GPL8490       Homo sapiens   IlluminaHumanMethylation27k
71 GPL10558       Homo sapiens               illuminaHumanv4
72 GPL11532       Homo sapiens   hugene11sttranscriptcluster
73 GPL13497       Homo sapiens         HsAgilentDesign026652
74 GPL13534       Homo sapiens  IlluminaHumanMethylation450k
75 GPL13667       Homo sapiens                        hgu219
76 GPL15380       Homo sapiens      GGHumanMethCancerPanelv1
77 GPL15396       Homo sapiens                     hthgu133b
78 GPL17897       Homo sapiens                     hthgu133a
这些包首先需要都下载
gpl_info=read.csv("GPL_info.csv",stringsAsFactors = F)
### first download all of the annotation packages from bioconductor
for (i in 1:nrow(gpl_info)){
  print(i)
  platform=gpl_info[i,4]
  platform=gsub('^ ',"",platform) ##主要是因为我处理包的字符串前面有空格
  #platformDB='hgu95av2.db'
  platformDB=paste(platform,".db",sep="")
  if( platformDB  %in% rownames(installed.packages()) == FALSE) {
    BiocInstaller::biocLite(platformDB)
    #biocLite(platformDB )
  }
}
下载完了所有的包， 就可以进行批量导出芯片探针与gene的对应关系！
for (i in 1:nrow(gpl_info)){
  print(i)
  platform=gpl_info[i,4]
  platform=gsub('^ ',"",platform)
  #platformDB='hgu95av2.db'
  platformDB=paste(platform,".db",sep="")
  if( platformDB  %in% rownames(installed.packages()) != FALSE) {
    library(platformDB,character.only = T)
    #tmp=paste('head(mappedkeys(',platform,'ENTREZID))',sep='')
    #eval(parse(text = tmp))
###重点在这里，把字符串当做命令运行
    all_probe=eval(parse(text = paste('mappedkeys(',platform,'ENTREZID)',sep='')))
    EGID <- as.numeric(lookUp(all_probe, platformDB, "ENTREZID"))
##自己把内容写出来即可
  }
}

 

 

一 31

画基因的外显子覆盖度图

Posted on by
一般情况下，我们得到了测序reads在基因组的比对情况文件bam格式的，里面的信息非常多，如果我想特定的查看某个基因的情况，那么我们可以选择IGV等可视化工具，但它并不是万能的，因为即使是一个基因，它也会有多个转录本，多个外显子。
所以，我们可以画它的外显子覆盖图，如下：横坐标是外显子的长度，纵坐标是测序深度，每一个小图都是一个外显子
 DMD.NM_000109
根据这个图，我们就可以很明显的看出，DMD基因NM_000109转录本的1，10-17号外显子缺失，用IGV一个个的看这些外显子区域，是同样的结果！可能是芯片捕获不到，也可能是样本本身变异，造成的大片段缺失。但是这个图的信息就非常有用！
那么，我们该如何画这样的图呢？
首先，我们需要找到需要探究的基因的全部转录信息，及外显子信息！
在hg19_refGene.txt里面会有，在UCSC里面可以下载，新手可能会比较麻烦，实在不行你去annovar的目录也可以找到！
1
那么，我们根据这个信息，就可以判断该基因的起始终止位点啦
然后用samtools的depth命令去找这个基因的全部片段的测序深度信息
最后再格式化成下面的三列数据
第一列是该外显子的坐标，从1到该外显子的成都
第二列是该外显子在该坐标的测序深度，通过samtools的depth命令得到
最后一列是该外显子的标记，从exon:79一直倒推到exon:1，因为该基因在染色体的负链，所以外显子顺序是反着的！
1 84 exon:79
2 84 exon:79
3 84 exon:79
4 85 exon:79
5 85 exon:79
6 86 exon:79
7 85 exon:79
8 87 exon:79
9 89 exon:79
10 91 exon:79
11 92 exon:79
12 95 exon:79
13 96 exon:79
14 96 exon:79
15 99 exon:79
16 99 exon:79
17 97 exon:79
最后根据这个txt文档，用R语言，很容易就画出上面那样的图片了！
这里面的信息量还是蛮大的！
APOB.NM_000384