29

芯片探针注释基因ID或者symbol,并对每个基因挑选最大表达量探针

在R里面实现这个功能其实非常简单,难的是很多packages经常会出现安装问题,更有的人压根不看芯片平台是什么,芯片对应的package是什么,就开始到处发问,自学能力实在是堪忧!

我前面有写目前所有bioconductor支持的芯片平台对应关系:通过bioconductor包来获取所有的芯片探针与gene的对应关系

但那其实是一个很笨的办法,得到所有的各式各样的探针ID与基因的对应关系,以为它绕路了,正常情况只需要在GEO里面找到芯片对应基因关系即可,没必要下载那么多package的,但是这样做的好处也是很明显的, 对很多初学者来说,如果package能解决的话,就省心很多,比如下面这个转换关系:

suppressPackageStartupMessages(library(CLL))
## 这个package自带了一个数据,是我们需要用的
data(sCLLex)  ## 这个数据里面有24个样本,分成两组,可以直接拿来测试差异基因分析
library(hgu95av2.db)  ## 一定要搞清楚自己的芯片是什么数据包
exprSet=exprs(sCLLex)  ##得到表达数据矩阵,但是矩阵的行名,是探针ID,无法理解,需要转换
##首先你取出所有的探针ID,#这里可以用三种方法来得到symbol,或者得到entrezID也可以
probeset=rownames(exprSet)
Symbol=as.character(as.list(hgu95av2SYMBOL[probeset]))
#annotate包提供              getSYMBOL( probeset ,"hgu95av2" )
#还可以用lookUp函数     lookUp( probeset , "hgu95av2", "SYMBOL")
#这些只是技巧而已啦
a=cbind.data.frame(Symbol,exprSet)
## 下面这个函数是对每个基因挑选最大表达量探针
rmDupID <-function(a=matrix(c(1,1:5,2,2:6,2,3:7),ncol=6)){
  exprSet=a[,-1]
  rowMeans=apply(exprSet,1,function(x) mean(as.numeric(x),na.rm=T))
  a=a[order(rowMeans,decreasing=T),]
  exprSet=a[!duplicated(a[,1]),]
  #exprSet=apply(exprSet,2,as.numeric)
  exprSet=exprSet[!is.na(exprSet[,1]),]
  rownames(exprSet)=exprSet[,1]
  exprSet=exprSet[,-1]
  return(exprSet)
}
exprSet=rmDupID(a)
对每个基因挑选最大表达量探针,只是一种处理方法而已,只是我一般处理芯片是这样做的,并不一定就是最好的!
23

强烈推荐用shiny来制作交互式网页展现数据

shiny是Rstudio公司推出的一个网页服务器,可以直接用R语言制作客户端和服务端程序来交互式的展现数据,而且有越来越丰富的扩展包可以借鉴使用,我觉得它的前途很光明,值得大家入坑!

虽然Rstudio公司吹嘘初学者无需具备html和css,js基础,但是个人认为还是多了解一下比较好,可以在w3cschool里面在线学习!

新手安装好shiny包之后里面自带了12个app例子,一边调试一边学习,很快就可以入门了。推荐一个教程,用shiny构建网页中文教程,想查看更详细的介绍和实例,请访问Shiny的官方主页。当然,你肯定需要要会R,这里有个R的FAQ教程。

如果你完全看完了上面那些,说明你已经真正入门了!

你现在可以去搜索一个shiny cheetsheet,然后背诵下来!!!!

然后你可以去shiny的github里面找到108个示例程序,一个个运行了解它们!!!

这时候你已经是高手啦!!!

接下来你应该取了解一个叫做dashboard的东西,熟练UI界面设计。

最后你再了解一些辅助包,帮助你用shiny更好的展示数据,主要是JS绘图插件

里面最出名的就是DT包了,一定要学!!!

最后你可以了解一下在rstudio上面分享五个免费的shiny网页程序,需要你搜索一些。

如果你还感兴趣的话,就可以自己整一个虚拟机Ubuntu或者centos系统,试着安装shiny的server,这个比较考验技术。

总结:你可以需要200个小时左右来完全掌握shiny

23

我的生信资源网盘分享

我重新整理了一下我当初生信入门时候收集的资料,希望对后学者有帮助,如果你以前储存了我的分享,请删掉以前的, 重新转这个完整版,以后不会更新了,入门的资料其实很多都是一样的!

公开链接:http://pan.baidu.com/s/1jIvwRD8 ,没有密码,直接根据需求自行下载

Continue reading

23

我的优酷视频网盘分享

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

我的自频道-优酷视频 

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

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

Continue reading

15

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

这个程序是我在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模块终极解决方案-下

其实可以手动下载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模块终极解决方案-上

不管别人怎么说,反正我是非常喜欢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 cpanmlocal::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

基因组各种版本对应关系

我是受到了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

开发单位:华大,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数据分析病毒整合位点

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

Continue reading

15

新的比对工具MOSAIK

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

Continue reading

12

bioconductor中文社区招募站长

我已经构建好bioconductor中文社区的雏形,大家可以进去看看!

写在前面

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

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

 

Continue reading

10

THis is a test for markdown

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数据库里面下载芯片数据

这个包跟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")