首先，干数据分析，至少需要会三种语言。
1. 用于绘图、高级数学计算的高级语言 （R SAS Matlab 其一）
2. 用于文本处理、转换格式的脚本语言 （perl python ruby php + bash + awk）
3. 用于快速运算、复杂算法运算的低级语言 (C/C++ Java)
当然这一点并不绝对。比如本人认为目前的 python 在合理使用扩展包的情况下，完全可以胜任以上三种任务。
用本人的 mbp ，用 R Python C++ 各写了两个代码，计算 student- T test 的 p-value
首先是使用 R语言，apply p_value 在data_frame中：
Program1 R - apply
setwd("~/Documents/【】Computer-science/R/李程02_r_matrix")
mm_data <- read.delim("GSE6477_series_matrix_expr.txt.xls", row.names=1)
mm_info <- read.delim("GSE6477_sample _info.xls", row.names=1)
mm.group <- mm_info$Type == "New" # new MM patients
normal.group <- mm_info$Type == "Normal donor" # normal samples
mm_data.sub <- mm_data[1:10, ] # test new code using a subset
apply(mm_data.sub, 1, t.test) # run t-test for each row (1 for row)
# but it's one sample t-test
one.t.test <- function(vec) {
t.test(vec[normal.group], vec[mm.group])$p.value
}
system.time( # time the computation
p_vector1 <- apply(mm_data, 1, one.t.test)
)
如果不用 apply 逐行计算，
Program2 R - byline
setwd("~/Documents/【】Computer-science/R/李程02_r_matrix")
mm_data <- read.delim("GSE6477_series_matrix_expr.txt.xls", row.names=1)
mm_info <- read.delim("GSE6477_sample _info.xls", row.names=1)
mm.group <- mm_info$Type == "New" # new MM patients
normal.group <- mm_info$Type == "Normal donor" # normal samples
mm.vec <- as.numeric(mm_data[1, mm.group])
normal.vec <- as.numeric(mm_data[1, normal.group])
t_test_result <- t.test(normal.vec, mm.vec)
mm_data.sub <- mm_data[1:10, ] # test new code using a subset
apply(mm_data.sub, 1, t.test) # run t-test for each row (1 for row)
# but it's one sample t-test
one.t.test <- function(vec) {
t.test(vec[normal.group], vec[mm.group])$p.value
}
p_vector2<-NULL
system.time(
for (i in c(1:length(rownames(mm_data)))){
p_vector2[[i]]<-apply(mm_data[i,],1,one.t.test)
}
)
write.table(p_vector2,"test2.result.xls",sep="\t",quote=FALSE)
如果用 python，方法是逐行读取，或者将各行的 两组数值放入一个
（ [1,2,3] , [5,6,7,8]
.......
[1,3,4], [5,6,9,10] ）
这样的 2d-numpy.array 中，最后用2d-numpy.array 计算correlation
或者逐行计算correlation
2d-numpy array算法：
Program3 Python numpy-2D array:
from __future__ import division
import numpy as np
from scipy import stats
import sys
def show_help():
print >>sys.stderr, "\n\tpython ",sys.argv[0],"GSE6477_sample_info.xls GSE6477_series_matrix_expr.txt.xls >../test3.result.xls"
def load_samp_type(sam_info,samp_type):
f_saminfo = open(sam_info,"r")
f_saminfo.readline()
for line in f_saminfo:
line = line.strip('\n')
f = line.split('\t')
samp_type[f[0]] = f[2]
f_saminfo.close()
def main():
try:
sam_info = sys.argv[1]
sam_exp = sys.argv[2]
except IndexError:
show_help()
sys.exit(1)
samp_type = {}
load_samp_type(sam_info,samp_type)
f_sam_exp = open(sam_exp,"r")
in_h = f_sam_exp.readline()
in_h = in_h.strip('\n')
in_h = in_h.strip('\r')
f_h = in_h.split('\t')
samp = f_h[1:]
sam_norm = []
sam_new = []
# print samp
for i,sam in enumerate(samp):
if samp_type[sam] == 'New':
sam_new.append(i+1)
elif samp_type[sam] == 'Normal donor':
sam_norm.append(i+1)
l_samp_new = []
l_samp_norm = []
l_gene = []
for line in f_sam_exp:
line = line.strip('\n')
line = line.strip('\r')
f = line.split('\t')
l_gene.append(f[0])
ll_new = []
ll_norm = []
for i1 in sam_new:
ll_new.append(float(f[i1]))
for i2 in sam_norm:
ll_norm.append(float(f[i2]))
l_samp_new.append(ll_new)
l_samp_norm.append(ll_norm)
np_ll_new = np.array(ll_new)
np_ll_norm = np.array(ll_norm)
# t,pval = stats.ttest_ind(np.array(ll_new),np.array(ll_norm))
# print t,pval,line
f_sam_exp.close()
np_samp_new = np.array(l_samp_new)
np_samp_norm = np.array(l_samp_norm)
l_pval = stats.ttest_ind(np_samp_new.T,np_samp_norm.T,equal_var=False)[1]
for i in range(0,len(l_gene)):
print l_gene[i],l_pval[i]
if __name__ == '__main__':
main()
Program4 Python by-line :
from __future__ import division
import numpy as np
from scipy import stats
import sys
def show_help():
print >>sys.stderr, "\n\tpython ",sys.argv[0],"GSE6477_sample_info.xls GSE6477_series_matrix_expr.txt.xls >../test3.result.xls"
def load_samp_type(sam_info,samp_type):
f_saminfo = open(sam_info,"r")
f_saminfo.readline()
for line in f_saminfo:
line = line.strip('\n')
f = line.split('\t')
samp_type[f[0]] = f[2]
f_saminfo.close()
def main():
try:
sam_info = sys.argv[1]
sam_exp = sys.argv[2]
except IndexError:
show_help()
sys.exit(1)
samp_type = {}
load_samp_type(sam_info,samp_type)
f_sam_exp = open(sam_exp,"r")
in_h = f_sam_exp.readline()
in_h = in_h.strip('\n')
in_h = in_h.strip('\r')
f_h = in_h.split('\t')
samp = f_h[1:]
sam_norm = []
sam_new = []
# print samp
for i,sam in enumerate(samp):
if samp_type[sam] == 'New':
sam_new.append(i+1)
elif samp_type[sam] == 'Normal donor':
sam_norm.append(i+1)
for line in f_sam_exp:
line = line.strip('\n')
line = line.strip('\r')
f = line.split('\t')
ll_new = []
ll_norm = []
for i1 in sam_new:
ll_new.append(float(f[i1]))
for i2 in sam_norm:
ll_norm.append(float(f[i2]))
np_ll_new = np.array(ll_new)
np_ll_norm = np.array(ll_norm)
# print f[0],np.mean(np_ll_new),np.std(np_ll_new),np.mean(np_ll_norm),np.std(np_ll_norm)
t,pval = stats.ttest_ind(np.array(ll_new),np.array(ll_norm),equal_var=False)
print f[0],pval
f_sam_exp.close()
if __name__ == '__main__':
main()
使用 c++ boost 库计算：
Program4 C++ boost :
g++ 编译方法：
g++ -O -o study-licheng4 main.cpp
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <vector>
#include <algorithm>
#include <functional>
#include <numeric>
#include <map>
#include <set>
#include <cmath>
#include <inttypes.h>
#include <boost/unordered_map.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string.hpp>
#include <boost/math/distributions/students_t.hpp>
using namespace std;
double two_samples_t_test_unequal_sd(
double Sm1,
double Sd1,
unsigned Sn1,
double Sm2,
double Sd2,
unsigned Sn2)
{
//
// Sm1 = Sample Mean 1.
// Sd1 = Sample Standard Deviation 1.
// Sn1 = Sample Size 1.
// Sm2 = Sample Mean 2.
// Sd2 = Sample Standard Deviation 2.
// Sn2 = Sample Size 2.
//
// A Students t test applied to two sets of data.
// We are testing the null hypothesis that the two samples have the same mean and that any difference if due to chance.
using namespace boost::math;
// Degrees of freedom:
double v = Sd1 * Sd1 / Sn1 + Sd2 * Sd2 / Sn2;
v *= v;
double t1 = Sd1 * Sd1 / Sn1;
t1 *= t1;
t1 /= (Sn1 - 1);
double t2 = Sd2 * Sd2 / Sn2;
t2 *= t2;
t2 /= (Sn2 - 1);
v /= (t1 + t2);
// cout << setw(55) << left << "Degrees of Freedom" << "= " << v << "\n";
// t-statistic:
double t_stat = (Sm1 - Sm2) / sqrt(Sd1 * Sd1 / Sn1 + Sd2 * Sd2 / Sn2);
//
// Define our distribution, and get the probability:
//
students_t dist(v);
double q = cdf(complement(dist, fabs(t_stat)));
return 2*q;
}
void usage()
{
cout << "./study-licheng4 <samp_info> <samp_matrix> " << endl;
cout << "./study-licheng4 GSE6477_sample_info.xls GSE6477_series_matrix_expr.txt.xls " << endl;
cout << " -h get help information" << endl;
exit (0);
} // void usage
void load_samp_type(
string f_samp_info,
boost::unordered_map< string,string > & M_samp_type)
{
ifstream infile;
infile.open(f_samp_info.c_str());
if ( ! infile )
{
cerr << "fail to open input file" << f_samp_info << endl;
exit(0);
}
string lineStr;
getline(infile,lineStr,'\n');
while (getline(infile,lineStr,'\n'))
{
if (lineStr[0] == ' ' || lineStr[0] == '\n'){
continue;
}
vector<string> LineVec;
boost::split(LineVec,lineStr,boost::is_any_of("\t\r\n"),boost::token_compress_on);
M_samp_type[LineVec[0]] = LineVec[2];
// cout << LineVec[0] << "\t" << LineVec[2] << "\t" << M_samp_type[LineVec[0]] << endl;
}
} // void load_samp_type
int main(int argc,char *argv[])
{
int c;
while ( (c=getopt(argc,argv,"h")) !=-1 )
{
switch(c){
case 'h' : usage();break;
default : usage();
}
}
if (argc <2) usage();
string f_samp_info = argv[optind++];
string f_samp_matrix = argv[optind++];
boost::unordered_map< string,string > M_samp_type;
load_samp_type(f_samp_info,M_samp_type);
ifstream infile;
infile.open(f_samp_matrix.c_str());
if (! infile){
cerr << "fail to open input file " << f_samp_matrix << endl;
exit(0);
}
string lineStr_h;
getline(infile,lineStr_h,'\n');
vector<string> LineVec_h;
boost::split(LineVec_h,lineStr_h,boost::is_any_of("\t\r\n"),boost::token_compress_on);
vector<string> V_samp;
for ( int i=1; i < LineVec_h.size()-1; i++ ){
V_samp.push_back(LineVec_h[i]);
}
vector<double> V_sam_new;
vector<double> V_sam_norm;
for ( int i=0; i < V_samp.size(); i++ ){
if ( M_samp_type[V_samp[i]] == "New" ){
V_sam_new.push_back(i+1);
}
else{
if ( M_samp_type[V_samp[i]] == "Normal donor" ){
V_sam_norm.push_back(i+1);
}
}
}
string lineStr;
while (getline(infile,lineStr,'\n'))
{
if (lineStr[0] == ' ' || lineStr[0] == '\n'){
continue;
}
vector<string> LineVec;
boost::split(LineVec,lineStr,boost::is_any_of("\t\r\n"),boost::token_compress_on);
string gene = LineVec[0];
vector<double> V_new;
vector<double> V_norm;
for (int i=0; i<V_sam_new.size(); i++){
int i1 = V_sam_new[i];
V_new.push_back( boost::lexical_cast<double>(LineVec[i1]) );
}
for (int i=0; i<V_sam_norm.size(); i++){
int i2 = V_sam_norm[i];
V_norm.push_back( boost::lexical_cast<double>(LineVec[i2]) );
}
double sum_new = accumulate(V_new.begin(), V_new.end(), 0.0);
double mean_new = sum_new / V_new.size();
vector<double> diff1(V_new.size());
transform( V_new.begin(), V_new.end(), diff1.begin(), bind2nd(minus<double>(),mean_new ) );
double sq_sum_new = inner_product(diff1.begin(), diff1.end(), diff1.begin(), 0.0);
double stdev_new = sqrt(sq_sum_new / (V_new.size() - 1) );
unsigned cnt_new = boost::lexical_cast<unsigned>(V_new.size());
double sum_norm = accumulate(V_norm.begin(), V_norm.end(), 0.0);
double mean_norm = sum_norm / V_norm.size();
vector<double> diff2(V_norm.size());
transform( V_norm.begin(), V_norm.end(), diff2.begin(), bind2nd(minus<double>(),mean_norm) );
double sq_sum_norm = inner_product(diff2.begin(), diff2.end(), diff2.begin(), 0.0);
double stdev_norm = sqrt(sq_sum_norm / (V_norm.size() - 1) );
unsigned cnt_norm = boost::lexical_cast<unsigned>(V_norm.size());
cout << "\n" << gene << "\t" << mean_new << "\t" << stdev_new << "\t" << cnt_new << "\t" << mean_norm << "\t" << stdev_norm << "\t" << cnt_norm << endl;
double p_value = two_samples_t_test_unequal_sd( mean_new, stdev_new, cnt_new, mean_norm, stdev_norm, cnt_norm );
cout << gene << "\t" << p_value << endl;
}
infile.close();
}// int main
最后，可以用 Rcpp写 R 的 c++ 插件，加速 R 的运算。R 本身对C++高度兼容，这一点弥补了其计算速度慢的这一大缺点。
Program 6 Rcpp
setwd("~/Documents/【】Computer-science/R/李程02_r_matrix")
Rcpp::sourceCpp('~/Documents/【】Computer-science/R/李程02_r_matrix/test/rcpp_hello_world.cpp')
mm_data <- read.delim("GSE6477_series_matrix_expr.txt.xls", row.names=1)
mm_info <- read.delim("GSE6477_sample _info.xls", row.names=1)
mm.group <- mm_info$Type == "New" # new MM patients
normal.group <- mm_info$Type == "Normal donor" # normal samples
one.t.test <- function(vec) {
p_value(as.numeric(vec[normal.group]), as.numeric(vec[mm.group]))
}
p_vector3<-NULL
system.time( # time the computation
p_vector3<-apply(mm_data, 1, one.t.test)
)
write.table(p_vector3,"test5.result.xls",sep="\t",quote=FALSE)
该程序调用的 rcpp_hello_world.cpp 如下：
#include <Rcpp.h>
#include <boost/math/distributions/students_t.hpp>
using namespace Rcpp;
// [[Rcpp::export]]
double p_value(NumericVector V_new, NumericVector V_norm){
using namespace boost::math;
int Sn1 = V_new.size();
int Sn2 = V_norm.size();
double Sm1 = mean(V_new);
double Sm2 = mean(V_norm);
double Sd1 = sd(V_new);
double Sd2 = sd(V_norm);
double v = Sd1 * Sd1 / Sn1 + Sd2 * Sd2 / Sn2;
v *= v;
double t1 = Sd1 * Sd1 / Sn1;
t1 *= t1;
t1 /= (Sn1 - 1);
double t2 = Sd2 * Sd2 / Sn2;
t2 *= t2;
t2 /= (Sn2 - 1);
v /= (t1 + t2);
double t_stat = (Sm1 - Sm2) / sqrt(Sd1 * Sd1 / Sn1 + Sd2 * Sd2 / Sn2);
students_t dist(v);
double q = cdf(complement(dist, fabs(t_stat)));
return 2*q;
}
以上 六个 程序在本人 mbp 的运行时间：
real 0m10.780s
user 0m10.475s
sys 0m0.170s
real 1m16.729s
user 1m16.222s
sys 0m0.402s
real 0m1.737s
user 0m1.564s
sys 0m0.154s
real 0m6.881s
user 0m6.818s
sys 0m0.059s
real 0m1.067s
user 0m0.949s
sys 0m0.094s
real 0m12.077s
user 0m10.345s
sys 0m0.567s
可见 速度 C++ > Python > R 。
同时应注意，在尚未使用任何 ipython 魔法加速计算的情况下、在未使用cython的情况下，合理使用的python 的 numpy scipy 其计算速度已经十分接近 c++ 了，时间差距只有不到2倍，远高于 R 的10倍（apply）。同时 R 在不apply 的情况下，使用 for循环速度慢的简直坑爹，是c++ 用时的一百倍。
此外注意到最后一个程序（Rcpp）虽然占用了12s，比不用Rcpp还慢。但是仔细看下 R 的系统时间 log ，发现实际上最后一个程序的时间主要花在了 source(Rcpp)上。抛去其编译时间，只看三个R 程序的计算pvalue用时，如下：
用户 系统 流逝
5.596 0.045 5.658
用户 系统 流逝
71.373 0.277 71.667
用户 系统 流逝
0.632 0.064 0.696
可见用 Rcpp 的计算速度是 R-apply t-test 的 10 倍，同时是 R用 for 循环的 100 倍。在不考虑编译、IO的情况下，速度接近 c++的计算，而其代码却简单了不少，省去了复杂的文本输入、关联数组，省去了需要小心数据计算方差、平均值时的 过大、过小溢出方面的考虑。因此对于R的大数据计算，非常推荐用 Rcpp 写程序。