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