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

Python计算Rank(秩)

在统计学中,有时候我们需要求一组数据的Rank,如非参检验、Spearman相关性等。Rank中文意思叫做顺序,更专业点的叫法是”秩”。这里介绍一下用Python求Rank。
先看一个简单的例子:

>>> a = [1.0, 5.6, 3.5, 9.7, 7.8]
>>> # rank_a = [1, 3, 2, 5, 4]
...
>>> rank_a = sorted(range(len(a)), key = a.__getitem__)
>>> rank_a
[0, 2, 1, 4, 3]
>>> rank_a = [i + 1 for i in rank_a]
>>> rank_a
[1, 3, 2, 5, 4]

上面的例子中,列表a中的元素是没有重复的,直接排序即可。如果有重复元素怎么办呢?最简单的办法是使用SciPy提供的rankdata()函数。

>>> from scipy.stats import rankdata
>>> rankdata([0, 2, 3, 2])
array([ 1. ,  2.5,  4. ,  2.5])
>>> rankdata([0, 2, 3, 2], method='min')
array([ 1,  2,  4,  2])
>>> rankdata([0, 2, 3, 2], method='max')
array([ 1,  3,  4,  3])
>>> rankdata([0, 2, 3, 2], method='dense')
array([ 1,  2,  3,  2])
>>> rankdata([0, 2, 3, 2], method='ordinal')
array([ 1,  2,  4,  3])

rankdata()函数返回一个array对象,默认的方法是average,即相同元素的Rank值取算术平均数,同时rankdata()函数还提供了其他方法。
如果觉得安装SciPy库太麻烦,我们也可以自己实现。下面的代码计算重复元素的Rank使用的方法是average,返回一个列表。

def rank_simple(vector):
    return sorted(range(len(vector)), key=vector.__getitem__)
 
def rankdata(vector):
    n = len(vector)
    ivec = rank_simple(vector)
    svec = [vector[rank] for rank in ivec]
    sumranks = 0
    dupcount = 0
    newvector = [0] * n
    for i in range(n):
        sumranks = sumranks + i
        dupcount = dupcount + 1
        if i == n - 1 or svec[i] != svec[i + 1]:
            averank = sumranks / float(dupcount) + 1
            for j in range(i - dupcount + 1, i + 1):
                newvector[ivec[j]] = averank
            sumranks = 0
            dupcount = 0