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

T检验靠谱吗?

对于实验生物学研究者来说,经常遇到的一个场景是,对照组测量3个值,实验组测量3个值,然后用T检验得出均值是否有显著性差异,然后做出结论。如下图所示:
bio_t

 

 

 

 

 

 

 

 

 

 

 

 

那么T检验靠谱吗?
本文采用模拟抽样的方法来探讨这件事情,由于本人统计学知识有限,还请各位读者批评指正。

我们通过构建两个随机总体,总体1的均值固定为1(模拟control),总体2的均值我们从依次取0.5到2.0(模拟treatment相对与control的变化),标准差我们取0.01、0.05、0.10(模拟测量误差1%,5%,10%)。总体数量为1000,每次随机抽取3个样本,抽样1000次,考察通过T检验得出正确结论的频率。代码如下:

import numpy as np
import pandas as pd
import seaborn as sns
from scipy import stats
import matplotlib.pyplot as plt
 
plt.xkcd()
sns.set(style="white", context="talk")
 
def norm_t(loc1, loc2, scale1, scale2, n):
    normSpace1 = stats.norm.rvs(loc = loc1, scale = scale1, size = 1000)
    normSpace2 = stats.norm.rvs(loc = loc2, scale = scale2, size = 1000)
    count = 0
    for i in range(1000):
        sample1 = np.random.choice(normSpace1, n)
        sample2 = np.random.choice(normSpace2, n)
        t, p = stats.ttest_ind(sample1, sample2)
        if (p < 0.05 and (sample1.mean() > sample2.mean()) == (loc1 > loc2)):
            count = count + 1
    return count / 1000
 
def plot(stds):
    df = pd.DataFrame()
    for std in stds:
        x = np.arange(0.5, 2.1, 0.1)
        y = [norm_t(loc1 = 1, loc2 = i, scale1 = std, scale2 = std, n = 3) for i in x]
        z = [std for i in x]
        df_std = pd.DataFrame([x, y, z]).T
        df_std.columns = ["difference", "sensitivity", "std"]
        if df.shape == (0, 0):
            df = df_std
        else:
            df = df.append(df_std)
    sns.pointplot(x = "difference", y = "sensitivity", hue = "std", data = df)
 
plot([0.01, 0.05, 0.10])

std_sensitivity
通过上图,我们发现,当标准差我们取0.01(模拟测量误差1%)时,增大到1.1倍或者减小到0.9倍,正确率可以达到100%。而当标准差取0.05和0.10(模拟测量误差5%和10%),效果就不理想了,本来是有差异的,还是有很大几率通过T检验得出错误的结论。显而易见的是,随着标准差的增加,T检验的效果变差了。

那么当标准差取0.10时,随着抽样样本数的增加,T检验判断正确的频率会是怎么样的呢?请看如下代码:

def plot2(diffs):
    df = pd.DataFrame()
    for diff in diffs:
        x = list(range(3, 60, 3))
        y = [norm_t(loc1 = 1, loc2 =1.1, scale1 = 0.10, scale2 = 0.10, n = i) for i in x]
        z = [diff for i in x]
        df_diff = pd.DataFrame([x, y, z]).T
        df_diff.columns = ["samle number", "sensitivity", "difference"]
        if df.shape == (0, 0):
            df = df_diff
        else:
            df = df.append(df_diff)
    sns.pointplot(x = "samle number", y = "sensitivity", hue = "difference", data = df)
 
plot2([0.5, 0.8, 1.2, 2.0])

n_sensitivity
通过上图,我们发现,即使标准差达到0.10,当样本数量达到30以上时,T检验的效果就非常好了。

那么如果两个总体没有差异,T检验得出有差异的错误结论的情况是什么样子的呢?我们构建两个均值为1容量为1000的总体,每次随机抽取3个样本,在不同标准差及样本数目的情况下,考察T检验得出正确结论的频率。请看代码:

def norm_t2(loc, scale1, scale2, n):
    normSpace1 = stats.norm.rvs(loc = loc, scale = scale1, size = 1000)
    normSpace2 = stats.norm.rvs(loc = loc, scale = scale2, size = 1000)
    count = 0
    for i in range(1000):
        sample1 = np.random.choice(normSpace1, n)
        sample2 = np.random.choice(normSpace2, n)
        t, p = stats.ttest_ind(sample1, sample2)
        if p > 0.05:
            count = count + 1
    return count / 1000
 
def plot3(stds):
    stds = [0.1, 0.5, 1]
    df = pd.DataFrame()
    for std in stds:
        x = list(range(3, 60, 3))
        y = [norm_t2(loc = 1, scale1 = std , scale2 = std, n = i) for i in x]
        z = [std for i in x]
        df_std = pd.DataFrame([x, y, z]).T
        df_std.columns = ["samle number", "specificity", "std"]
        if df.shape == (0, 0):
            df = df_std
        else:
            df = df.append(df_std)
    sns.pointplot(x = "samle number", y = "specificity", hue = "std", data = df)
    plt.ylim(0, 1.2)
 
plot3([0.1, 0.5, 1])

specificity
通过上图,我们惊讶地发现,当两个总体没有差异时,不管标准差是多少和取样数目是多少,T检验都有非常好的表现,即得到差异不显著的的结论。

这里构建的总体都是正太总体,其实对于其他分布的总体,结论是一样的,这里就不展示了。

通过上面的探索,我们可以得到如下结论:
1. 对于两个未知总体,通过T检验考察均值是否有显著差异,如果得出结论是差异不显著,那么进一步分析这些数据的标准差是否太大了,考虑是否增加抽样样本数做进一步分析。
2. 对于两个未知总体,通过T检验考察均值是否有显著差异,如果得出结论是差异显著,那么请相信它吧!

如果差异不显著,我们通过增大样本数量,使得差异显著;如果差异显著,那就是一个好结果嘛!
哈哈,T检验是实验生物学家的利器啊!

使用JBrowse搭建基因组浏览器服务

JBrowse是一款快速的,可以嵌入的基因组浏览器,它采用JavaScript和HTML5开发,同时提供了perl脚本处理数据。JBrowse采用了Ajax技术,主要的计算在客户端上完成, 对服务器资源的占用非常低,速度非常快,可以轻松部署上G的数据。JBrowse支持GFF3, BED, FASTA, Wiggle, BigWig, BAM, VCF (with tabix), REST等众多数据格式,BAM, BigWig,VCF格式的数据可以不需要转换而直接使用。 JBrowse对系统依赖非常低,只需要一个web服务器,因此部署起来非常方便。

下载JBrowse最新版本,然后安装:

unzip JBrowse-1.12.1.zip
cd JBrowse-1.12.1
./setup.sh

./setup.sh这一步需要安装很多perl依赖库,需要的时间可能有点久。如果报错,请查看setup.log文件,一般是系统可能缺某个依赖库。

安装好之后,我们将JBrowse-1.12.1下所有文件拷贝到网站根目录,这里我们使用Nginx作为web服务器。有关Nginx的安装和配置,参加另一篇博客文章:http://blog.biochen.com/archives/797

sudo mv * /usr/local/nginx/html/

访问http://127.0.0.1可以看到欢迎页面:
JBrowse_0

访问http://127.0.0.1/index.html?data=sample_data/json/volvox可以看到示例页面:
JBrowse_1

BED格式转GTF格式

# -*- coding: utf-8 -*-
"""
Created on Fri Dec 25 19:43:00 2015
 
@author: biochen
"""
 
import sys, getopt, csv
 
# get parameter
opts, args = getopt.getopt(sys.argv[1:], "hi:o:s:")
input_file = ""
output_file = ""
source = ""
for op, value in opts:
    if op == "-i":
        input_file = value
    elif op == "-o":
        output_file = value
    elif op == "-s":
        source = value
    elif op == "-h":
        print("Usage: python3 BED2GTF.py -i input.bed -o output.gtf -s source")
        sys.exit()
 
 
BED_in_file = open(input_file, "r")
GTF_out_file = open(output_file, "w")
BED = csv.reader(BED_in_file, dialect = "excel-tab")
GTF = csv.writer(GTF_out_file, dialect = "excel-tab", quotechar = "'")
 
for item in BED:
    transcript = []
    transcript.append(item[0])
    transcript.append(source)
    transcript.append("transcript")
    transcript.append(int(item[1]) + 1)
    transcript.append(int(item[2]) + 1)
    transcript.append(item[4])
    transcript.append(item[5])
    transcript.append(".")
    transcript.append('transcript_id "' + item[3] + '";')
    GTF.writerow(transcript)
    exon_number = int(item[9])
    exon_size = list(map(int, (filter(lambda x:x and len(x.strip()) > 0, item[10].split(",")))))
    exon_start = list(map(int, (filter(lambda x:x and len(x.strip()) > 0, item[11].split(",")))))
    for i in range(exon_number):
        exon = []
        exon.append(item[0])
        exon.append(source)
        exon.append("exon")
        exon.append(int(item[1]) + exon_start[i] + 1)
        exon.append(exon[3] + exon_size[i] - 1)
        exon.append(item[4])
        exon.append(item[5])
        exon.append(".")
        exon.append('transcript_id "' + item[3] + '\";')
        GTF.writerow(exon)
 
BED_in_file.close()
GTF_out_file.close()

GTF格式转BED格式

 # -*- coding: utf-8 -*-
"""
Created on Wed Dec 16 11:03:05 2015
 
@author: biochen
"""
 
import sys, getopt, csv, re
 
opts, args = getopt.getopt(sys.argv[1:], "hi:o:")
for op, value in opts:
    if op == "-i":
        GTF_in_File_Name = value
    elif op == "-o":
        BED_in_File_Name = value
    elif op == "-h":
        print("Usage: python3 GTF2BED.py -i input.gtf -o output.bed")
        sys.exit()
 
GTF_in_file = open(GTF_in_File_Name, "r")
BED_out_file = open(BED_in_File_Name, "w")
GTFs = csv.reader(GTF_in_file, dialect = "excel-tab")
BEDs = csv.writer(BED_out_file, dialect = "excel-tab")
 
transcripts = {}
for GTF in GTFs:
    if GTF[2] == "exon":
        transcript_id = re.findall('transcript_id "([^;]*)"', GTF[8])[0]
        if transcript_id in transcripts:
            transcripts[transcript_id][6].append([int(GTF[3]), int(GTF[4])])
        else:
            transcripts[transcript_id] = []
            transcripts[transcript_id].append(GTF[0])
            transcripts[transcript_id].append(GTF[3])
            transcripts[transcript_id].append(GTF[4])
            transcripts[transcript_id].append(GTF[5])
            transcripts[transcript_id].append(GTF[6])
            transcripts[transcript_id].append(transcript_id)
            transcripts[transcript_id].append([[int(GTF[3]), int(GTF[4])]])
 
for transcript in transcripts:
    transcripts[transcript][6].sort()
    transcript_start = transcripts[transcript][6][0][0]
    transcript_end = transcripts[transcript][6][len(transcripts[transcript][6]) - 1][1]
    exon_sizes = ""
    exon_starts = ""
    for exon in transcripts[transcript][6]:
        exon_size = exon[1] - exon[0] + 1
        exon_start = exon[0] - transcript_start
        exon_sizes = exon_sizes + str(exon_size) + ","
        exon_starts = exon_starts + str(exon_start) + ","
    BED = []
    BED.append(transcripts[transcript][0])
    BED.append(transcript_start - 1)
    BED.append(transcript_end)
    BED.append(transcripts[transcript][5])
    BED.append(transcripts[transcript][3])
    BED.append(transcripts[transcript][4])
    BED.append(transcript_start - 1)
    BED.append(transcript_end)
    BED.append("0, 0, 0")
    BED.append(len(transcripts[transcript][6]))
    BED.append(exon_sizes)
    BED.append(exon_starts)
    BEDs.writerow(BED)
 
GTF_in_file.close()
BED_out_file.close()

HTSeq

安装C/C++编译环境
Ubuntu缺省情况下,并没有提供C/C++的编译环境,运行以下命令安装:
sudo apt-get install build-essential

安装setuptools
setuptools可以更方便的创建和发布Python包,HTSeq的安装依赖setuptools
下载setuptools-20.3.1.tar.gz
tar zxvf setuptools-20.3.1.tar.gz
cd setuptools-20.3.1/
sudo python setup.py install

安装Cython
Cython是Python的C扩增,可以让Python程序更加高效。HTSeq依赖cython。
sudo apt-get install cython

安装Numpy
Numpy是Python科学计算基本库,提供了许多高级的数值编程工具。HTSeq依赖Numpy。
sudo apt-get install python-numpy

安装Zlib
Zlib是zlib是提供数据压缩用的函数库。pysam依赖Zlib。
下载zlib128.zip
unzip zlib128.zip
cd zlib-1.2.8/
./configure
make
sudo make install

安装pysam
pysam可以操作很多基因组数据,包括SAM/BAM/VCF/BCF/BED/GFF/GTF/FASTA/FASTQ。HTSeq依赖pysam。
下载pysam-0.9.0.tar.gz
tar zxvf pysam-0.9.0.tar.gz
cd pysam-0.9.0/
sudo python setup.py install

安装HTSeq
下载HTSeq-0.6.1.tar.gz
tar zxvf HTSeq-0.6.1.tar.gz
cd HTSeq-0.6.1/
sudo python setup.py install

报错