佛祖的灯芯

《大话西游》里有一段:
紫霞来到了驿站的外面,在院落里趴着一只黑狗,正在盯着紫霞,一个老者正在来回的奔波着,不断的在院子里取火,然后又跑进屋里去点灯。
“各位客官,请稍等一下,马上就好了。”
灯里面只有灯油,没有灯芯,当然没办法点着了,紫霞走了过去:“老伯啊,你这样是点不着的,你只有油,没有灯芯哎。”
“我知道,可是灯芯不见了,有什么办法呢?”
“再买一根喽。”
紫霞和二郎神化作的老头说了几句,老头忽然直起身来,看着紫霞,声音骤然一沉。
“你说得对,可是如来佛祖的那盏日月神灯的灯芯不见了,到哪里去买?”
没有灯芯点不着灯,那人心呢?

PhastCons、PhyloP和保守性

全基因组多序列比对的主要目的是揭示基因组中进化上保守的和非保守的区段以及它们的分布。但是,在碱基层次观察多个全基因组的多序列比对结果和序列保守性非常不方便;其次,许多功能区段的保守性介乎高度保守和完全不保守之间;再者,保守性要能予以方便的显示。上述三点使UCSC基因组浏览器使用两个软件PhastCons和PhyloP把由MULTIZ产生的多序列比对结果转换成两种单一的保守性计分和显示,这两个方法都基于已知的种系树结构和利用一个称作phylo-HMM的隐马尔可夫模型种系分析方法。

PhyloP只考虑比对的当前列而PhastCons同时也考虑比对列的相邻列,这使得PhastCons对于保守区段的出现更敏感而PhyloP对于保守区段的界定更有效。PhyloP和PhastCons之间的另一个区别是,PhyloP能够识别加速进化和保守进化的位点,它们分别产生正的和负的计分,而PhastCons的计分总是一个0至1之间的正值。

PhyloP和PhastCons多物种比对文件可以从UCSC基因组浏览器的网站上下载,他们还有提供了bigWigAverageOverBed这个工具。

用法非常简单,只是输入文件in.bw、in.bed和输出out.tab必须是这个名字。

使用MCScanX分析基因组共线性区块

MCScanX采用改进了的MCScan算法,分析基因组内或者基因组间的共线性区块。它利用两个物种蛋白质blastp比对结果,再结合这些蛋白质基因在基因组中的位置,得到两个物种基因组的共线性区块。如果是分析基因组内的共线性区块,物种内蛋白质自己比对自己就好了。
一、MCScanX安装
下载MCScanX的安装包MCScanX.zip,可以在linux系统和Mac OS进行编译:

$ unzip MCscanX.zip
$ cd MCScanX
$ make

如果是在64位系统进行编译,msa.h、dissect_multiple_alignment.h、detect_collinear_tandem_arrays.h这三个源文件的开头需要添加#nclude 。
MCScanX包括MCScanX、MCScanX_h、duplicate_gene_classifier三个主程序,位于主文件夹中;还有12个下游分析程序位于downstream_analyses文件夹中。MCScanX 检测共线性区域,并比对到参考染色体上。MCScanX_h 和MCScanX类似,只不过输入文件是成对的用tab分隔的同源基因。下游分析程序主要用于可视化输出。
二、输入文件的准备
我们需要两个物种的蛋白信息,还有对应蛋白在基因组中的位置。推荐从Ensembl下载蛋白质信息,它的注释信息非常全。如下是Homo_sapiens.GRCh38.pep.all.fa(下载于Ensembl,人的全部蛋白质信息)的一条蛋白质的信息:
>ENSP00000452494.1 pep chromosome:GRCh38:14:22449113:22449125:1 gene:ENSG00000228985.1 transcript:ENST00000448914.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:TRDD3 description:T cell receptor delta diversity 3 [Source:HGNC Symbol;Acc:HGNC:12256]
该蛋白质名ENSP00000452494.1,位于14号染色体上,起始位置22449113,终止位置22449125。
blast建库:

$ makeblastdb -in a.fa -dbtype prot -parse_seqids -out adb

blast比对:

$ blastp -query b.fa -db adb -out xyz.blast -evalue 1e-10 -num_threads 16 -outfmt 6 -num_alignments 5

按照软件说明书的要求,-evalue 1e-10是将e value卡在1e-10,-num_alignments 5是取最好的5个比对结果,-outfmt 6是输出格式为tab分隔的比对结果。
蛋白质基因在基因组上的位置文件为gff文件。这里的gff文件不是平常用到的那种,需要修改一下。文件是一个tab分隔的,只有4列:
sp# gene starting_position ending_position
sp是两个字母,后面的#号表示染色体编号,这个文件自己生成一下好了。
三、使用MCScanX分析基因组共线性区块
假设我们的blast文件为xyz.blast,gff文件为xyz.gff,且都放在一个文件夹下,我们可以简单运行一下命令:

$ MCScanX xyz

程序输出的结果有两个xyz.collinearity文件和xyz.html目录。
xyz.collinearity,成对的共线性区域,如下图所示:

xyz.html目录里面有很多小文件,文件名称是参考基因组染色体编号。第一列是每个基因位点的复制深度,第二列是基因参考染色体,红色部分是串联基因,后面黄色部分的内容是比对上的共线性区域,只有哪些比对上的基因会被展现出来。html文件可以用网页浏览器打开。如下图所示:

下面是一些更详细的参数:

四、下游分析程序
下游分析程序很多,主要是用来出图的,就不一一介绍了。
dot_plotter.java画散点图

$ java dot_plotter -g gff_file -s collinearity_file -c control_file -o output_PNG_file

circle_plotter.java画圆圈图

$ java circle_plotter -g gff_file -s collinearity_file -c control_file -o output_PNG_file

输入除了MCScanX的结果和gff文件,还需要一个ctl文件,这个文件自己按照MCScanX的说明去写即可,告诉程序图要画成什么样子。

更多用法参见官方的说明书。
http://chibba.pgml.uga.edu/mcscan2/documentation/manual.pdf