使用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

配置Visual Studio Code的C++开发环境

Windows下最好的IDE是Visual Studio,可惜它太臃肿,在我的笔记本电脑上不能流畅运行。和Atom、Sublime类似,微软推出了轻量级的编辑器Visual Studio Code,免费且跨平台。

1. 安装TDM-gcc
本来想装MinGW的,可是在线安装总是掉线,就换成了了。TDM-gcc可以下载离线安装包,而且自动添加环境变量(path),非常省心!

2. 安装Visual Studio Code并安装cpptools扩展
在官网下载Visual Studio Code,安装好之后,Ctrl+Shift+X打开扩展管理器,搜索到cpptools并安装。
cpptools

3. 配置c_cpp_properties.json,添加include路径
随便打开一个C++源文件,点击一下#include旁边的黄灯,然后点击Add include path to settings,VS code会弹出c_cpp_properties.json文件,我们在这个文件中添加include路径。
include1
include2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4. 配置tasks.json
按下Ctrl+shift+B,按照VS Code提示点击“配置任务运行程序”,选择“Others 运行任意外部命令的示例”,这时会弹出tasks.json文件。类似如下,修改tasks.json文件。
{
“version”: “0.1.0”,
“command”: “g++”,
“isShellCommand”: true,
“showOutput”: “always”,
“args”: [“-g”, “main.cpp”]
}
task_json1
task_json2
 
 
 
 
 
 
 
 
 
 
 
然后就可以通过按下Ctrl+shift+B编译C++程序了。

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

Ubuntu14.04 安装、配置nginx

Nginx是一款轻量级、高性能web服务器,它占用内存少,并发能力强。越来越多的网站使用Nginx作为web服务器。
Nginx依赖以下模块:
1. gzip模块需要 zlib 库
2. rewrite模块需要 pcre 库
3. ssl 功能需要openssl库
Nginx可以使用系统中已经安装的库,也可以自行指定。

下载Nginx,pcre,openssl,zlib源文件,然后通过如下命令安装:

tar zxvf nginx-1.11.0.tar.gz
tar zxvf pcre-8.38.tar.gz
tar zxvf openssl-1.0.1j.tar.gz
tar zxvf zlib-1.2.8.tar.gz
cd nginx-1.11.0
./configure --with-pcre=../pcre-8.38 --with-openssl=../openssl-1.0.1j --with-zlib=../zlib-1.2.8
sudo make & make install

Nginx默认安装在/usr/local/下,也可以在configure时使用–prefix参数指定安装路径。

Nginx安装后,Nginx文件夹的目录结构如下:

.
|-- conf
|   |-- fastcgi.conf
|   |-- fastcgi.conf.default
|   |-- fastcgi_params
|   |-- fastcgi_params.default
|   |-- koi-utf
|   |-- koi-win
|   |-- mime.types
|   |-- mime.types.default
|   |-- nginx.conf
|   |-- nginx.conf.default
|   |-- scgi_params
|   |-- scgi_params.default
|   |-- uwsgi_params
|   |-- uwsgi_params.default
|   `-- win-utf
|-- html
|   |-- 50x.html
|   `-- index.html
|-- logs
|   |-- access.log
|   `-- error.log
`-- sbin
    `-- nginx

conf文件夹下是配置文件,html文件夹是网站默认根目录,logs文件夹下是日志文件,sbin文件夹下是Nginx可执行程序。

启动Nginx:

sudo /usr/local/nginx/sbin/nginx             #启动Nginx
sudo /usr/local/nginx/sbin/nginx -s stop     #停止Nginx
sudo /usr/local/nginx/sbin/nginx -s reload   #重启Nginx

启动Nginx后,访问127.0.0.1,可以看到欢迎页面。
nginx

nginx.conf是主配置文件
worker_process表示工作进程的数量,一般设置为cpu的核数
worker_connections表示每个工作进程的最大连接数
server{}块定义了虚拟主机
listen监听端口
server_name监听域名
location{}指定网站文件的位置
index指定首页index文件的名称,可以配置多个,以空格分开。如有多个,按配置顺序查找。

Python运算符和优先级

运算符 说明
yield X 生成器函数结果(返回send()值)
lambda args1:X 匿名函数生成器(调用后返回X)
X if Y else Z 三元选择(当Y为真时计算X,否则计算Z)
X or Y 逻辑或
X and Y 逻辑与
not X 逻辑非
X in Y, X not in Y 成员是否在可迭代对象中
X is Y, X is not Y 对象等价测试
X Y, X >= Y 大小比较
X == Y, X != Y 相等运算符
X | Y 集合并
X ^ Y 集合对称差
X ^ Y 集合交
X << Y, X >> Y 将X左/右移Y位
X + Y, X – Y 加,减
X * Y, X / Y, X // Y, X % Y 乘,除,整除,取余
-X, +X 负,正
~X 按位非补(倒置)
X ** Y 乘方
X[i] 索引(序列,映射,其他)
X[i:j:k] 切片
X(args) 调用(函数,方法,类和其他可调用对象)
X.attr 属性引用
(…) 元组,表达式,生成器表达式
[…] 列表,列表内涵
{…} 字典,集合,字典与集合内涵

参考:
《Python袖珍指南》(第五版)

Python常用转义字符

转义字符 含义 转义字符 含义
\newline 忽略换行 \t 水平制表符
\\ 反斜杠(\) \v 竖直制表符
\’ 单引号(‘) \N{id} Unicode dbase id
\” 双引号(“) \uhhhh Unicode 16位十六进制
\a 蜂鸣(Bell) \Uhhhhhhhh Unicode 32位十六进制
\b 空格 \xhh Hex(最多两位数字)
\f 换页 \ooo Octal(达到三位数字)
\n 换行 \0 Null(非字符串结尾)
\r 回车 \other 非转义字符

这里的大部分转义字符跟linux系统保持一致!
\Uhhhhhhhh严格采用8个十六进制数字(h)。\u和\U只可以用在Unicode字符串文字中。

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()