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

28

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

其实国际三大主流生物信息学数据库运营单位都出了自己的基因组坐标转换,它们分别是 (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

28

使用可视化工具MutationMapper来看看基因上面突变的分布

how to generate lollipop diagrams ?
这个工具是网页版的,不需要下载,只需上传两列数据即可,就是基因上面的第几个氨基酸是如何突变的这样的信息
网页就可以把这些信息画在基因上面,而且注释到domain和cosmic数据库,挺好用的!
Hugo_Symbol HUGO symbol for the gene TP53
Protein_Change Amino acid change V600E
结果也很容易理解!
我输入了4个基因,网页就会对每个基因画一个图,其实ABCA1这个基因有6个突变,就都画在该基因上面了,基因突变分布在基因的哪个domain也看得一清二楚!
1
我的输入文件是这样的!
2
非常好用:
  • Support mutation data with annotated protein effects
  • Mutation diagram/lollipop view
  • Mutation table view
  • 3D structure view if available
而且pfam数据库本身也提供了这个功能:

Pfam provides an online tool to not only generate the domain information in JSON format, but to draw the lollipop diagram using javascript as well. They have more information here: http://pfam.xfam.org/help#tabview=tab9

IMHO, not as pretty as cBioPortal's but it gets you close to a solution.

EDIT / SHAMELESS PLUG: After seeing the data available and how easy it'd be, I made my own quick tool to fetch the data and draw the diagram for me in a style similar to cBioPortal - feel free to fork it and add features: https://github.com/pbnjay/lollipops

Example output (w/ labels per the comments)

3

甚至还有高手用D3.JS写了一个能实现同样需求的模块

We found ourselves in the same need, we wanted such a plot (JavaScript). Thus, I add our solution, Mutations Needle Plot. The library creates an SVG image (with D3), which then may be downloaded.

You will npm in order to be able to install & run the library.

Examples may be found in the snippets folder or also the index.html - The one displayed here below

类似的高手很多:
My colleague @SolenaLS recently asked me to write something like this: (uniprot+SVG+javascript : )http://lindenb.github.io/pages/uniprot/paintsvg.html
而且ensembl数据库也提供类似的功能
Ensemble also gives similar plot.

 

22

用RankComp的思想来做差异基因分析

是福建医科大学的学者开发的,
文章里面详细讲解了他们的这个差异分析的统计学原理
大意就是找到同一组织的normal样本的表达量数据,几百个,这样就可以分析2万基因之间的互相配对,检测表达量是否在几百个样本里面稳定的不一样!

我现在还不是很确定这个方法,只是试一试,欢迎与我交流对该方法的讨论!

文章是:

Wang H, Sun Q, Zhao W, et al. Individual-level analysis of differential expression of genes and pathways for personalized medicine[J]. Bioinformatics, 2014: btu522.

他们把它写成了一个R包,可以下载使用,但是必须用R2.15.2版本,我用了一下,不好用!

We can download the R code for in http://bioinformatics.oxfordjournals.org/content/31/1/62/suppl/DC1

他们这个程序真心不好用,但是很容易看懂算法,可以自己用R语言写一个来实现同样的过程!

比如A基因在几百个样本里面表达量都是3左右,而B基因都是5左右,而且满足99%的A表达量高于B,那么这就是一个稳定的基因对!
一般2万基因之间可以配成2亿个基因对,其中稳定的大概有10%~40%
然后我们对每个疾病样本都可以进行检验,看看这样稳定的基因对是否被改变!
比如,我们拿到一个疾病样本的2万个基因的表达量,我们挑取一个,如果它有100个稳定的up的基因对,100个稳定的down的基因对
那么,我看看这些基因对是否被改变了,如果这样还有70基因对在该疾病样本里面仍然是up,60个是down,那么我用Fisher精确检验的结果是
4
这个基因在该疾病样本,相对于normal pool并没有差异表达!当然检验得到的P值最后可以做FDR校验。
依次这样,把所有的gene都分析完,就知道这个样本有哪些差异的gene了。
介绍完原理,我们拿一个具体的例子来看看吧:
首先我们下载一个2008年的一个人的肝脏表达数据(Gene Expression in Human Liver),都是正常组织,共427个样本。
不过这个芯片比较小众,是默克医药公司定制化的, 需要下载探针对应gene的文件!
我们读取GSE9588这个数据,得到表达矩阵,然后计算rank矩阵,然后计算得到comp矩阵
5
> table(rank_comp)
rank_comp
     down        no        up
    58479 465752098     58479
>
不知道为什么这个数据,stable的那么少,不知道是不是出了什么问题!
其实我的程序都是对的了,只是因为这个数据集已经不是纯粹的表达量数据了,而是这427个样本的数据都减去了某个样本的表达量。
这样每个个体的基因之间的表达量排序就会被干扰,导致得到的稳定基因对非常少!!!
但是,我后来下载了GTEx的表达数据,拿那里面的normal组织样本表达量来做,可以得到非常多的稳定基因对。
实际代码大概是:
得到正常组织的表达矩阵:
然后计算表达矩阵的rank,得到各个样本自己的基因排序情况,得到排序矩阵!
处理排序矩阵,每个基因对之间都算一下是否稳定,得到稳定性描述矩阵!
然后根据每个疾病个体的基因表达情况,来循环每个基因, 看看该基因是否差异!

 

14

用broad出品的软件来处理bam文件几次遇到文件头错误

报错如下:ERROR MESSAGE: SAM/BAM file input.marked.bam is malformed: SAM file doesn't have any read groups defined in the header.  The GATK no longer supports SAM files without read groups !

有些人遇到的是bam的染色体顺序不一样,还有可能是染色体的名字不一样,比如>1和>chr1的区别,虽然很傻,但是遇到这样问题的还不少!
还有一些人是遇到基因组没有dict文件,也是用picard处理一下就好。

大部分人是在GATK遇到的,我是在RNA-SeQC遇到的,不过原理都是一样的。
都是因为做alignment的时候并未添加头信息,比如:
bwa samse ref.fa my.sai my.fastq > my.sam
samtools view -bS my.sam > my.bam
samtools sort my.bam my_sorted
java -jar ReordereSam.jar I=/path/my_sorted.bam O=/path/my_reordered.bam R=/path/ref.fa
通过这个代码可以得到排序好的bam,但是接下来用GATK就会报错
java -jar GenomeAnalysisTK.jar -T DepthOfCoverage -R /paht/ref.fa -I /path/aln_reordered.bam
就是因为没有头信息,group相关信息,解决方法有两种:
bwa samse -r @RG\tID:IDa\tSM:SM\tPL:Illumina ref.fa my.sai my.fastq > my.sam
java -jar AddOrReplaceReadGroups I=my.bam O=myGr.bam LB=whatever PL=illumina PU=whatever SM=whatever
一种是比对的时候就加入头信息,这个需要比对工具的支持。
第二种是用picard工具来修改bam,推荐用这个!虽然我其实并不懂这些头文件信息是干嘛的, 但是broad开发的软件就是需要!希望将来去读PHD能系统性的学习一些基础知识!

 

14

用RNA-SeQC得到表达矩阵RPKM值

这个软件不仅仅能做QC,而且可以统计各个基因的RPKM值!尤其是TCGA计划里面的都是用它算的
一、程序安装
直接在官网下载java版本软件即可使用:http://www.broadinstitute.org/cancer/cga/tools/rnaseqc/RNA-SeQC_v1.1.8.jar
但是需要下载很多注释数据
clipboard

二、输入数据

clipboard

箭头所指的文件,一个都不少,只有那个rRNA.tar我没有用, 因为这个软件有两种使用方式,我用的是第一种
三、软件使用
软件的官网给力例子,很容易学习:
RNA-SeQC can be run with or without a BWA-based rRNA level estimation mode. To run without (less accurate, but faster) use the command:
java -jar RNASeQC.jar -n 1000 -s "TestId|ThousandReads.bam|TestDesc" -t gencode.v7.annotation_goodContig.gtf -r Homo_sapiens_assembly19.fasta -o ./testReport/ -strat gc -gc gencode.v7.gc.txt 
我用的就是这个例子,这个例子需要的所有文件里面,染色体都是没有chr的,这个非常重要!!!
代码如下:
 java -jar RNA-SeQC_v1.1.8.jar  \
-n 1000 \
-s "TestId|ThousandReads.bam|TestDesc" \
-t gencode.v7.annotation_goodContig.gtf \
-r ~/ref-database/human_g1k_v37/human_g1k_v37.fasta  \
-o ./testReport/ \
-strat gc \
-gc gencode.v7.gc.txt \
To run the more accurate but slower, BWA-based method :
java -jar RNASeQC.jar -n 1000 -s "TestId|ThousandReads.bam|TestDesc" -t gencode.v7.annotation_goodContig.gtf -r Homo_sapiens_assembly19.fasta -o ./testReport/ -strat gc -gc gencode.v7.gc.txt -BWArRNA human_all_rRNA.fasta
Note: this assumes BWA is in your PATH. If this is not the case, use the -bwa flag to specify the path to BWA
四、结果解读
运行要点时间,就那个一千条reads的测试数据都搞了10分钟!
出来一大堆突变,具体解释,官网上面很详细,不过,比较重要的当然是RPKM值咯,还有QC的信息
clipboard
TCGA数据里面都会提供由RNA-SeQC软件处理得到的表达矩阵!
Expression
  • RPKM data are used as produced by RNA-SeQC.
  • Filter on >=10 individuals with >0.1 RPKM and raw read counts greater than 6.
  • Quantile normalization was performed within each tissue to bring the expression profile of each sample onto the same scale.
  • To protect from outliers, inverse quantile normalization was performed for each gene, mapping each set of expression values to a standard normal.
软件的主页是:
 
 
 
 
11

批量运行GSEA,命令行版本

之前用过有界面的那种,那样非常方便,只需要做好数据即可,但是如果有非常多的数据,每次都要点击文件,点击下一步,也很烦,不过,,它既然是java软件,就可以以命令行的形式来玩转它!

能够命令行运行了,就很容易批量啦

一、程序安装

直接在官网下载java版本软件即可:http://software.broadinstitute.org/gsea/downloads.jsp

二、输入数据

需要下载gmt文件,自己制作gct和cls文件,或者直接下载测试文件p53

见:http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats

1
三、运行命令

hgu95av2的芯片数据,只有一万多探针,所以很快就可以出结果

 java -cp gsea2-2.2.1.jar  -Xmx1024m xtools.gsea.Gsea   -gmx c2.cp.kegg.v5.0.symbols.gmt  \
 -res P53_hgu95av2.gct  -cls P53.cls   -chip  chip/HG_U95Av2.chip  -out result -rpt_label p53_example
但是一般我们都用默认的即可

2
里面报错说有些探针找不到,不要管它

四、输出数据
3
自己看官网去理解这些结果咯!
需要下载的数据如下:

首先需要下载 Molecular Signatures Database (MSigDB),一般选择C2的kegg,BioCarta 还有Reactome

都是gmt格式的文件!
CP: Canonical pathways
(browse 1330 gene sets)
Gene sets from the pathway databases. Usually, these gene sets are canonical representations of a biological process compiled by domain experts. details Download GMT Files
original identifiers
gene symbols
entrez genes ids
CP:BIOCARTA: BioCarta gene sets
(browse 217 gene sets)
Gene sets derived from the BioCarta pathway database (http://www.biocarta.com/genes/index.asp). Download GMT Files
original identifiers
gene symbols
entrez genes ids
CP:KEGG: KEGG gene sets
(browse 186 gene sets)
Gene sets derived from the KEGG pathway database (http://www.genome.jp/kegg/pathway.html). Download GMT Files
original identifiers
gene symbols
entrez genes ids
CP:REACTOME: Reactome gene sets
(browse 674 gene sets)
Gene sets derived from the Reactome pathway database (http://www.reactome.org/). Download GMT Files
original identifiers
gene symbols
entrez genes ids
然后做出表达数据gct文件和cls表型文件~
然后就可以直接运行了
如果是芯片数据,第一列是芯片探针,那么还需要下载chip数据:ftp://ftp.broadinstitute.org/pub/gsea/annotations
06

R语言软件的各种旧版本下载

这种编程语言,我之前以为只有包,或者模块什么的才烦人,没想到最近还碰到了版本的问题。
一般来说, 去R的官网,可以下载到基于各种操作系统的最新版R软件,但是在某些特别的情况下,我们可能需要下载以前的旧版本。这时候,想用百度这个破东西搜索出来下载地址简直是痴心妄想。
比如,我想搜索2.15.2这个版本,在Google里面很容易找到了
如果进入bin目录,可以看到各种操作系统的软件https://cran.r-project.org/bin/ ,里面都是可执行的文件,有下面几种操作系统
[DIR] linux/ 23-Jan-2008 19:47 -
[DIR] macos/ 19-Apr-2005 09:45 -
[DIR] macosx/ 12-Dec-2015 09:04 -
[DIR] windows/ 24-Feb-2012 18:41 -
大多数情况下面我们会用window这个平台的,毕竟,可视化嘛!
很多时候,我们也会需要linux平台的https://cran.r-project.org/bin/linux/ 但是linux平台本身就比较复杂
[DIR] debian/ 15-Dec-2015 02:06 -
[DIR] redhat/ 27-Jul-2014 21:12 -
[DIR] suse/ 16-Feb-2012 15:09 -
[DIR] ubuntu/ 06-Jan-2016 04:05 -
我用的Ubuntu比较多,即使选择了ubuntu,里面还有好几个选项,因为Ubuntu也有各种版本!
[DIR] precise/ 06-Jan-2016 04:03 -
[DIR] trusty/ 06-Jan-2016 04:04 -
[DIR] vivid/ 06-Jan-2016 04:04 -
[DIR] wily/ 06-Jan-2016 04:05 -

所以如果,大家是想在linux系统里面安装旧版本的R,建议大家直接下载c源码,直接三部曲就可以安装啦!

[DIR] R-0/ 04-Oct-2004 10:20 -
[DIR] R-1/ 04-Oct-2004 19:02 -
[DIR] R-2/ 01-Mar-2013 09:11 -
[DIR] R-3/ 10-Dec-2015 09:13 -
windows系统里面一般是是exe的可执行文件,直接安装,不需要源码变异,只需要不停的下一步即可
R 3.2.2 (August, 2015)
R 3.2.1 (June, 2015)
R 3.2.0 (April, 2015)
R 3.1.3 (March, 2015)
R 3.1.2 (October, 2014)
R 3.1.1 (July, 2014)
R 3.1.0 (April, 2014)
R 3.0.3 (March, 2014)
R 3.0.2 (September, 2013)
R 3.0.1 (May, 2013)
R 3.0.0 (April, 2013)
R 2.15.3 (March, 2013)
R 2.15.2 (October, 2012)
R 2.15.1 (June, 2012)
R 2.15.0 (March, 2012)
R 2.14.2 (February, 2012)
R 2.14.1 (December, 2011)
R 2.14.0 (November, 2011)
R 2.13.2 (September, 2011)
R 2.13.1 (July, 2011)
R 2.13.0 (April, 2011)
R 2.12.2 (February, 2011)
R 2.12.1 (December, 2010)
R 2.12.0 (October, 2010)
R 2.11.1 (May, 2010)
R 2.11.0 (April, 2010)
R 2.10.1 (December, 2009)
R 2.10.0 (October, 2009)
R 2.9.2 (August, 2009)
R 2.9.1 (June, 2009)
R 2.9.0 (April, 2009)
R 2.8.1 (December, 2008)
R 2.8.0 (October, 2008)
R 2.7.2 (August, 2008)
R 2.7.1 (June, 2008)
R 2.7.0 (April, 2008)
R 2.6.2 (February, 2008)
R 2.6.1 (November, 2007)
R 2.6.0 (October, 2007)
R 2.5.1 (July, 2007)
R 2.5.0 (April, 2007)
R 2.4.1 (December, 2006)
R 2.4.0 (October, 2006)
R 2.3.1 (June, 2006)
R 2.3.0 (April, 2006)
R 2.2.1 (December, 2005)
R 2.2.0 (October, 2005)
R 2.1.1 (June, 2005)
R 2.1.0 (April, 2005)
R 2.0.1 (November, 2004)
R 2.0.0 (October, 2004)
R 1.9.1 (June, 2004)
R 1.8.1 (November, 2003)
R 1.7.1 (June, 2003)
R 1.6.2 (January, 2003)
Installer for R 1.5.1 (June, 2002)
Installer for R 1.4.1 (January, 2002)
Installer for R 1.3.1 (September, 2001)
Binary files for R 1.2.2 (March, 2001)
Binary files for R 1.0.0 (February, 2000)

 

06

picnic对拷贝数变异检测芯片数据进行分析

这里的拷贝数变异检测芯片指的是Affymetrix Genome-Wide Human SNP Array 6.0

cel数据,处理成segment及genotype数据

一、程序安装
这本来是一个matlab程序,但是有linux版本,需要安装matlab编译环境
   clipboard
下载解压之后首先安装matlab环境:
./MCRInstaller.bin -console
因为我的服务器没有界面,所以我用了-console执行安装程序,很简单就安装好了
clipboard2
其中,cdf目录下面是我们从http://www.affymetrix.com 里面下载的Library FilesGenome-Wide Human SNP Array 6.0 (zip, 246 MB)
celConverter目录下面是一个java程序,可以把芯片出来的cel文件转为flat file
其中Matlab_running 是我执行安装的matlab环境
其余两个脚本run_HMM.sh  run_preprocessing.sh是picnic程序的第二步和第三步
二、输入数据准备
我随便找了两个snp6.0芯片的raw data
-rw-rw-r-- 1 jmzeng jmzeng 66M Dec 30 06:30 GSM1949207.CEL
-rw-rw-r-- 1 jmzeng jmzeng 66M Sep  9 11:08 GSM887898.CEL
三、程序使用

程序主要分为三个步骤来实现
Step 1: Convert the binary *.cel file to a flat file
Step 2: Normalise and estimate the ploidy of the genome and the level of normal contamination
Step 3: Segment the data and produce the genotype
第一步:
暂时还不需要matlab工作环境
示例的code是:java -Xmx2G -jar CelFileConverter.jar -m Snp6FeatureMappings.csv -c 'cdf_file_including_path' -s 'directory name
of cel files' -t rootDir/outdir/raw
我的code是:

jmzeng@ubuntu:/home/jmzeng/bio-soft/picnic/c_code/celConverter$ java -jar CelFileConverter.jar -m Snp6FeatureMappings.csv -c ../cdf/GenomeWideSNP_6.Full.cdf -s ../celFiles/ -t ../result/ 
Loading CDF file: ../cdf/GenomeWideSNP_6.Full.cdf
Loading feature mapping file: Snp6FeatureMappings.csv
Processing CEL file: /home/jmzeng/bio-soft/picnic/c_code/celConverter/../celFiles/GSM1949207.CEL
Created feature intensity file: /home/jmzeng/bio-soft/picnic/c_code/celConverter/../result/GSM1949207.feature_intensity
Processing CEL file: /home/jmzeng/bio-soft/picnic/c_code/celConverter/../celFiles/GSM887898.CEL
Created feature intensity file: /home/jmzeng/bio-soft/picnic/c_code/celConverter/../result/GSM887898.feature_intensity

然后在输出目录就有了一个feature_intensity后缀的文本文件,约110M
-rw-rw-r-- 1 jmzeng jmzeng 105M Dec 30 06:35 GSM1949207.feature_intensity
-rw-rw-r-- 1 jmzeng jmzeng 109M Dec 30 06:35 GSM887898.feature_intensity
第二步:
这里需要用matlab啦,但是这里经常出现库的问题!!!
示例的code是
sh run_preprocessing.sh mcr_dir cell_name info_dir feature_int_dir normalised_outdir outdir sample_type in_pi
我的code是:
jmzeng@ubuntu:/home/jmzeng/bio-soft/picnic/c_code$ sh run_preprocessing.sh Matlab_running/v710/ GSM1949207 info/ result/ result/output result/
------------------------------------------
Setting up environment variables
---
LD_LIBRARY_PATH is .:Matlab_running/v710//runtime/glnxa64:Matlab_running/v710//bin/glnxa64:Matlab_running/v710//sys/os/glnxa64:Matlab_running/v710//sys/java/jre/glnxa64/jre/lib/amd64/native_threads:Matlab_running/v710//sys/java/jre/glnxa64/jre/lib/amd64/server:Matlab_running/v710//sys/java/jre/glnxa64/jre/lib/amd64/client:Matlab_running/v710//sys/java/jre/glnxa64/jre/lib/amd64
My Own Exception: Fatal error loading library /home/jmzeng/bio-soft/picnic/c_code/Matlab_running/v710/bin/glnxa64/libmwmclmcr.so Error: libXp.so.6: cannot open shared object file: No such file or directory
第三步:
也需要用matlab,库的问题必须解决!!!
我另外一台服务器的结局办法是用v714的matlab
示例的code是
sh run_HMM.sh /nfs/team78pc2/kwl_temp/segments/PICNIC/C/release/Matlab_Compiler_Runtime/v710 A01_CGP_PD3945a.feature_intensity '/nfs/team78pc3/KWL/segments/PICNIC/matlab/C/release/info/' '/nfs/team78pc2/kwl_temp/segments/PICNIC/data/normalized/' '/nfs/team78pc2/kwl_temp/segments/PICNIC/data/' '10' '0.33598' '1.9915' '0.40997'
我的code是:无所谓了,这个服务器不知道怎么回事,总是出现库文件的问题,而这个问题需要root权限,我懒得弄了,我在其它的服务器上面都木有问题的!!!
所以我就换了MATLAB,毕竟,这个软件本来就是matlab版本的,第一步照旧,第二步在matlab里面运行!
我这里只用GSM1949207做测试吧!!!
Genomic DNA was extracted from saliva, peripheral blood, or fibroblast cell lines using the QIAamp DNA Blood Mini Kit or QIAamp DNA Mini Kit. DNA quality and quantity was assessed using a Nanodrop Spectrophotometer and agarose gel electrophoresis.
打开matlab,进入picnic目录
重复第二步,输入:
preprocessing('GSM1949207.feature_intensity','info\','result\raw\','result\output\','result\')
需要7个参数,我这个是cell lines数据,所以后面两个参数省略不写!info文件夹自己下载放在picnic目录,其中result\raw 存放你第一步的结果文件!
这一步运行,会比较久!
3
同时可以看到程序在我的result目录里面新增加了两个目录用来存放结果,不过这一步的结果还是中间文件,就不解释了!
然后再重复第三步,输入:

 HMM('GSM1949207.feature_intensity','info\','result\output\','result\',10,0,2.0221,0.40997)

这个软件的参数很没规律,最好不要像说明书那些用output做文件夹名字,不然,很容易出错。
这一步好像也很耗时间!
4

重要的结果有两个:
56
十二 29

用firehose_get 来下载所有TCGA寄存在broad的数据

该软件是broad institute写的一个数据接口,主要是供他人下载TCGA的所有寄存在broad institute的免费数据,主要是level3,level4的数据。(说错了,好像只有level4的数据,就是可以发文章的分析结果及图片)
软件下载地址:https://confluence.broadinstitute.org/display/GDAC/Download

懂它的使用规则,编码规则即可:
就是一个很简单的shell脚本而已,根据几个用户自定义参数来选择性的下载数据。
clipboard
我们可以用-t这个参数来指定下载的数据类型,可以是mut/rna/mutsig/gistic等各种数据,至于这些单词代表什么意义,需要自己去看说明书啦
还可以指定时间,截止到什么时间的数据!
它支持的癌症种类:

ACC  BLCA  BRCA  CESC  COAD  COADREAD  DLBC  ESCA  
	GBM  HNSC  KICH  KIRC  KIRP  LAML  LGG  LIHC  
	LUAD  LUSC  OV  PAAD  PANCANCER  PANCAN8  PANCAN12  PRAD  
	READ  SARC  SKCM  STAD  THCA  UCEC  UCS
这些癌症种类的简称,也是可以去官网里面看到的!官网:http://gdac.broadinstitute.org

 

十二 29

自动化出网页报告的R语言包-Nozzle

根据测试代码输出的报告如下:
clipboard
这些报告里面的元素,都是测试代码添加的,很容易就能理解
十二 29

用Mutation-Assessor软件来看突变位点对基因或者蛋白功能的影响

这是一个在线工具,非常好用,对snp位点进行注释,看看该突变是否影响蛋白功能,一定要收藏!!!

官网:http://mutationassessor.org/

也应该是有standalone版本,我没有去找,不过网页版就很好用了,只需要复制粘贴进去自己想分析的数据,按照一定的格式即可,比如:
clipboard
该软件就能分析出该突变位点发生在BRCA2这个基因上面,对氨基酸的改变也能写出来,对蛋白功能的改变等选项都是可以自由定制化的。
输入数据非常广泛:
The server accepts list of variants, one variant per line, plus optional text describing your variants,
in genomic coordinates, "+" strand assumed :
<genome build>,<chromosome>,<position>,<reference allele>,<substituted allele>
Genome build is optional (build 18 assumed), accepted values: 'hg18' and 'hg19'
Examples:

hg19,13,32912555,G,T   BRCA2
hg18,7,55178574,G,A   GBM
7,55178574,G,A   GBM

or in protein space: <protein ID> <variant> <text>, where <protein ID> can be :

1. Uniprot protein accession (i.e. EGFR_HUMAN)
2. NCBI Refseq protein ID (i.e. NP_005219)

examples:

EGFR_HUMAN R521K
EGFR_HUMAN R98Q Polymorphism
EGFR_HUMAN G719D disease
NP_000537 G356A
NP_000537 G360A dbSNP:rs35993958
NP_000537 S46A Abolishes phosphorylation

ID types can be mixed in one list in any way.

 

十二 11

用excel表格做差异分析

其实主要要讲的不是用excel来做差异分析,只是想讲清楚差异分析的原理,用excel可视化的操作可能会更方便理解,而且想告诉大家,其实生物信息学分析,本来就很简单的,那么多软件,只有你理解了原理,你自己就能写出来的!

首先,还是得到表达矩阵,下面绿色的样本是NASH组,蓝色的样本是normal组

image001

我们进行差异分析,很简单,就是看两组的表达值,是否差异,而检验的方法就是T检验。

=AVERAGE(D2:L2)    ##求NASH组的平均表达量

=AVERAGE(M2:S2)    ###求normal的平均表达量

=T2-U2             ##计算得到logFOLDchange值

=AVERAGE(D2:S2)    ###得到所有样本的平均表达量

=T.TEST(D2:L2,M2:T2,2,3)  ###用T检验得到两个组的表达量的差异显著程度。

简单检查几个值就可以看到跟limma包得到的结果差不多。

image002

 

十二 11

用limma包对芯片数据做差异分析

下载该R语言包,然后看说明书,需要自己做好三个数据(表达矩阵,分组矩阵,差异比较矩阵),总共三个步骤(lmFit,eBayes,topTable)就可以啦

image001

首先做第一个数据,基因表达矩阵!

自己在NCBI里面可以查到下载地址,然后用R语言读取即可

exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)

rownames(exprSet)=exprSet[,1]

exprSet=exprSet[,-1]

image002

然后做好分组矩阵,如下

image003

然后做好,差异比较矩阵,就是说明你想把那些组拿起来做差异分析,如下

image004

最后输出结果:

我进行了6次比较,所以会输出6次比较结果

image005

最后打开差异结果,解读,说明书如下!

 

 

image006

在我的github有完整代码

 

十一 05

在linux系统里面安装matlab运行环境mcr

matlab毕竟是收费软件,而且是有界面的。所以搞生物信息的都用R和linux替代了,但是很多高大上的单位,比如大名鼎鼎的broadinstitute,是用matlab的,所以他们开发的程序也会以matlab代码的形式发布。但是考虑到大多研究者用不起matlab,或者不会用,所以就用linux系统里面安装matlab运行环境来解决这个问题,我们仍然可以把人家写的matlab程序,在linux命令行下面,当做一个脚本来运行!

比如,这次我就需要用broadinstitute的一个软件:Mutsig,找cancer driver gene的,http://www.broadinstitute.org/cancer/cga/mutsig_run,但是我看了说明才发现,它是用matlab写的,所以我要想在我的服务器用,就必须按照安装matlab运行环境,在官网可以下载:http://www.mathworks.com/products/compiler/mcr/

Capture3

我这里选择的是R2013a (8.1),下载之后解压是这样的,压缩包约四百多M

Capture1

然后直接在解压后的目录里面运行那个install即可,然后如果你的linux可以传送图像,那么就会想安装windows软件一样方便!如果你的linux是纯粹的命令行,那么,就需要一步步的命令行交互,选择安装地址,等等来安装了。

记住你安装之后,会显示一些环境变量给你,请千万要记住,然后自己去修改自己的环境变量,如果你忘记了,就需要搜索来解决环境变量的问题啦!安装之后是这样的:

Capture2

请记住你的安装目录,以后你运行其它matlab相关的程序,都需要把这个安装目录,当做一个参数传给你的其它程序的!!!

如果你没有设置环境变量,就会出各种各样的错误,用下面这个脚本可以设置

其中MCRROOT一般是$path/biosoft/matlab_running/v81/ 这样的东西,请务必注意,LD_LIBRARY_PATH非常重要,非常重要,非常重要!!!!

MCRROOT=$1
echo ---
LD_LIBRARY_PATH=.:${MCRROOT}/runtime/glnxa64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/bin/glnxa64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRROOT}/sys/os/glnxa64;
MCRJRE=${MCRROOT}/sys/java/jre/glnxa64/jre/lib/amd64 ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/native_threads ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/server ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE}/client ;
LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${MCRJRE} ;
XAPPLRESDIR=${MCRROOT}/X11/app-defaults ;
export LD_LIBRARY_PATH;
export XAPPLRESDIR;
echo LD_LIBRARY_PATH is ${LD_LIBRARY_PATH};

如果你没有设置正确,那么会报一下的错误!

error while loading shared libraries: libmwmclmcrrt.so.8.1: cannot o pen shared object file: No such file or directory

error while loading shared libraries: libmwlaunchermain.so: cannot o pen shared object file: No such file or directory

等等!!!

 

24

用freebayes来call snps

step1:,软件安装
wget http://clavius.bc.edu/~erik/freebayes/freebayes-5d5b8ac0.tar.gz
tar xzvf freebayes-5d5b8ac0.tar.gz
cd freebayes
make
一个小插曲,安装的过程报错:/bin/sh: 1: cmake: not found
所以我需要自己下载安装cmake,然后把cmake添加到环境变量

首先下载源码包http://www.cmake.org/cmake/resources/software.html

wget http://cmake.org/files/v3.3/cmake-3.3.2.tar.gz

 解压进去,然后源码安装三部曲,首先 ./configu  然后make 最后make install  

cmake 会默认安装在 /usr/local/bin 下面,这里需要修改,因为你可能没有 /usr/local/bin 权限,安装到自己的目录,然后把它添加到环境变量!

step2:准备数据

an alignment file (in BAM format)
a reference genome in (uncompressed) FASTA format.
正好我的服务器里面有很多
不过,该软件也可以出了一个测试数据集
wget http://bioinformatics.bc.edu/marthlab/download/gkno-cshl-2013/chr20.fa
wget ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/data/NA12878/alignment/NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam
wget ftp://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/data/NA12878/alignment/NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam.bai

用这个代码就可以下载千人基因组计划的NA12878样本的第20号染色体相关数据啦

step3:运行命令
网站给出的实例是:
freebayes -f chr20.fa \
    NA12878.chrom20.ILLUMINA.bwa.CEU.low_coverage.20121211.bam >NA12878.chr20.freebayes.vcf

一般就用默认参数即可

 
step4:输出结果解读
没什么好解读的了,反正是vcf文件,都看烂了,就那些东西
不过该软件的作者倒是拿该软件与broad用GATK做出的NA12878样本的突变数据做了比较

 

13

使用Seq2HLA进行HLA分型

基于高通量测序数据进行HLA分型的软件挺多的,比较老的有三个,作者分别是Boegel et al.Kim et al.Major et al.,然后他们都被OptiType这个软件的作者被批评了,我这里先介绍Kim et al的seq2HLA使用方法,以下是它的一些链接。

功能概述

seq2HLA is a computational tool to determine Human Leukocyte Antigen (HLA) directly from existing and future short RNA-Seq reads. It takes standard RNA-Seq sequence reads in fastq format as input, uses a bowtie index comprising known HLA alleles and outputs the most likely HLA class I and class II types, a p-value for each call, and the expression of each class.

软件简介

Type of tool     Program

Nature of tool          Standalone

Operating system   Unix/Linux, Mac OS X

Language        Python, R

Article     (Boegel et al., 2013) HLA typing from RNA-Seq sequence reads. Genome medicine.

PubMed http://www.ncbi.nlm.nih.gov/pubmed/23259685

URL          https://bitbucket.org/sebastian_boegel/seq2hla

源代码,下载并安装

https://bitbucket.org/sebastian_boegel/seq2hla/src

http://tron-mainz.de/tron-facilities/computational-medicine/seq2hla/

第一版是这样的

image001

第二版是这样的

image002

只有第二版才支持gz压缩包格式的fastq,而且不需要指定length了

其中reference文件夹下面的是发布这个软件的团体已经制备好来的HLA库文件

image003

下载即可使用,前提是你的系统其它环境都OK

用法:

python seq2HLA.py -1 <readfile1> -2 <readfile2> -r "<runname>" [-p <int>]* [-3 <int>]**

image004

很简单,-1和-2指定我们的双端测序数据即可,可以是压缩包格式的(自动调用gzip),-r的输出目录,会输出7个文件,需要一个个解读,-p指定线程数给bowtie用的,-3是指定需要trim几个低质量碱基。

但是运行这个软件的要求非常多,需要安装好python和R,而且还有版本限制,需要安装好biopython而且还必须是双端测序,而且当前文件夹下面的reference文件夹下面必须有参考基因组的bowtie索引,而且系统必须安装好了bowtie,还需要在快捷方式里面!

我这里用的是第二版的

image006

所以,我用的也是第二版改进的命令。非常好用,我这里用的是一个外显子测序数据,是hiseq2500测的PE100

python seq2HLA.py -1 ../../6-exon/PC3-1.read1_Clean.fastq.gz -2 ../../6-exon/PC3-1.read2_Clean.fastq.gz -r PC3

貌似输出文件太多了一点

#Output:#The results are output to stdout and to textfiles. Most important are:

#i) <prefix>-ClassI.HLAgenotype2digits => 2 digit result of Class I

#ii) <prefix>-ClassII.HLAgenotype2digits => 2 digit result of Class II

#iii) <prefix>-ClassI.HLAgenotype4digits => 4 digit result of Class I

#iv) <prefix>-ClassII.HLAgenotype4digits => 4 digit result of Class II

#v) <prefix>.ambiguity => reports typing ambuigities (more than one solution for an allele possible)

#vi) <prefix>-ClassI.expression => expression of Class I alleles

#vii) <prefix>-ClassII.expression => expression of Class II alleles

根据文献,我简单看了一下,文件的确好复杂,不过我们只需要看输出日志即可

-----------2 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

A       A*68   7.287148e-05   A*24   0.03680272

B       B*52   0.1717737       B*53   0.3952319

C       C*12   0.03009331     hoz("C*14")     0.6783964

Calculation of locus-specific expression ...

BC1-1/BC1-1-ClassI.bowtielog

A: 7.93 RPKM

C: 9.75 RPKM

B: 8.35 RPKM

The digital haplotype is written into BC1-1/BC1-1-ClassI.digitalhaplotype3

-----------4 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

!A     A*68:01 7.287148e-05   A*24:02 0.03680272

!B     B*52:01 0.1717737       B*53:01'       0.6542288

!C     C*12:02 0.03371717     C*12:02 0.6783964

上面的HLA的class I的数据结果

接下来是class II的数据结果,是不是很简单呀!

-----------2 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

DQA     DQA1*01 0.1511134       DQA1*02 0

DQB     DQB1*02 0.02321615     DQB1*05 0.42202

DRB     DRB1*15 2.595144e-05   DRB1*07 0.321219

Calculation of locus-specific expression ...

BC1-1/BC1-1-ClassII.bowtielog

DQB1: 4.47 RPKM

DRB1: 5.59 RPKM

DQA1: 0.44 RPKM

-----------4 digit typing results-------------

#Locus Allele 1       Confidence     Allele 2       Confidence

!DQA   DQA1*01:02'     0.1511134       DQA1*02:01     0.0

!DQB   DQB1*02:01'     0.02321615     DQB1*05:01     0.42202

!DRB   DRB1*15:02'     2.595144e-05   DRB1*07:01     0.321219