Author Archive for 写长城的诗

《R语言与Bioconductor生物信息学应用》已经出版

《R语言与Bioconductor生物信息学应用》一书主要是使用R语言来解决生物信息学问题。主要作者是两位在生物信息学界有极深造诣的科研人士,高山和欧剑虹,这二位的博客地址如下:
http://blog.sciencenet.cn/u/gaoshannankai
http://pgfe.umassmed.edu/ou/
我本人也完成了其中部分章节的写作。此书的特点在于:

  • 从实际课题出发,提出解决这个问题的思路,结合用到的原理或基础知识,但更偏重整个解决问题的框架和流程,选用R这种简单易学但功能强大的语言,把讲解延伸到具体程序代码,让读者100%经历整个课题研究过程。
  • 最大的创新点是:实际课题直接来作者发表的SCI文章,全部都是真枪实弹,不杜撰所谓“实际应用”。国内外尚未见到与SCI文章紧密结合的生物信息书籍。
  • 本书是多名R领域专家(全部都是一线科研工作者)通过互联网联手写作。在前期网上调研的基础上,尽量在本书内突出大多数人普遍关心而又难找相关资料的问题。
  • 所见即所得,学到的知识可以通过简单编程(仅仅代码拷贝粘贴)加以实现,印象深刻,学了不会忘。提出三板斧学习法,让无基础的人也能编程。
  • 本书作者通过QQ群直接面向读者答疑,并且共享了大量的参考资料和习题答案。使用正版书籍的读者都可以入群享受最好的服务。不在书中罗列基础知识凑字数,也不使用光盘,既能减轻读者负担,又能保护环境。



此书详细目录如下

第一章 R基础知识 4
1.1 什么是R 4
1.1.1 R语言的起源 4
1.1.2 R语言的特点 5
1.1.3 R语言的主要用途 7
1.1.4 R语言的应用现状和发展趋势 10
1.2 R的下载与安装 12
1.2.1 主程序的下载与安装 12
1.2.2 扩展包的下载与安装 14
1.2.3 R语言的集成开发环境 16
1.2.4 R主程序和扩展包的管理与升级 19
1.3 R语言快速入门 21
1.3.1 从哪里入手开始学习R 21
1.3.2 三板斧搞定R语言 21
1.3.3 一个例子来说明三板斧 22
1.4 一些简单的语法知识 23
1.4.1 什么是编程 23
1.4.2 变量 23
1.4.3 函数 24
1.4.4 综合案例 25
1.5 本章源代码详解及小结 25
1.5.1 例1-1 25
1.5.2 例1-2 26
1.5.3 例1-3 27
1.5.4 例1-4 28
1.5.5 小结 30
参考文献: 31
第二章 生物信息学基础知识 4
2.1中心法则-生物信息流 4
2.1.1 生物大分子 4
2.1.2 中心法则 7
2.1.3 基因组、转录组和蛋白质组 8
2.1.4 非编码RNA和microRNA 9
2.2测序与序列分析 10
2.2.1 DNA测序技术 10
2.2.2 第二代测序技术的应用领域 12
2.2.3 序列分析 13
2.2.4 序列比对和相似性搜索 14
2.2.5 分子进化和系统发生树 15
2.3基因表达分析 17
2.3.1基因表达的检测方法 17
2.3.2 基因表达数据分析 18
2.3.3基因表达差异的显著性分析 19
2.3.4基因本体论分析 20
2.3.5通路分析 22
2.4注释、统计与可视化 22
2.4.1 注释与ID映射 23
2.4.2 统计与可视化 23
参考文献: 24
第三章 R在生物信息学中的简单应用 2
3.1 一个序列分析课题 2
3.1.1 课题背景 2
3.1.2 研究目的与实验设计 2
3.1.3 数据获取与处理流程 3
3.2 用R包(非bioconductor)实现课题 4
3.2.1 定义全部函数(例3-1) 4
3.2.2 课题实现 13
3.2.3 源代码详解与小结 16
3.3 用R包(bioconductor)实现课题一 21
3.3.1 重新设计数据处理流程和全部函数(例3-2) 21
3.3.2 课题实现 25
3.3.3源代码详解与小结 26
3.4 用R包(bioconductor)实现课题二 26
3.41 重新设计数据处理流程和全部函数(例3-3) 26
3.4.2 课题实现 35
3.4.3源代码详解与小结 35
第四章 Bioconductor简介 1
4.1 什么是Bioconductor 2
4.1.1 Bioconductor的起源 2
4.1.2 Bioconductor主要特点 2
4.2 Bioconductor包的分类介绍 5
4.2.1 三大板块分类介绍 5
4.2.2 软件包的进一步介绍 6
4.2.3 按照应用领域分类 9
4.3 从R到Bioconductor的跨越:Biostrings, BiomaRt以及AnnotationDbi包 12
4.3.1 应用Biostrings处理生物序列 13
4.3.2 应用BiomaRt获取实验数据与注释信息 21
4.3.3 应用AnnotationDbi生成注释包 27
参考文献: 32
第五章Bioconductor分析基因芯片数据
5.1 快速入门 2
5.2 基因芯片基础知识 3
5.2.1 探针组 3
5.2.2 主要的芯片文件格式 4
5.3 基因芯片数据预处理 5
5.3.1 数据输入 6
5.3.2 质量控制 7
5.3.3 背景校正、标准化和汇总 17
5.3.4 预处理的一体化算法 20
5.4 基因芯片数据分析 24
5.4.1 选取差异表达基因 24
5.4.2 注释 27
5.4.3 统计分析及可视化 28
5.5 芯片处理实际课题一 39
5.5.1 课题背景 39
5.5.2 数据集与预处理 40
5.5.3 R程序与代码讲解 41
5.6 芯片处理实际课题二 42
5.6.1 课题背景 42
5.6.2 数据集与处理过程 43
5.6.3 R程序与代码讲解 43
5.7 芯片处理实际课题三 44
5.7.1 课题背景 44
5.7.2 数据集与处理过程 45
5.7.3 R程序与代码讲解 46
参考文献: 48
第六章 Bioconductor分析RNA-seq数据 4
6.1示例课题介绍 4
6.1.1课题背景 4
6.1.2数据集和处理过程 4
6.2高通量测序基础知识 5
6.2.1高通量测序原理 5
6.2.2测序的质量分数 9
6.2.3高通量测序文件格式 12
6.3 RNA-seq技术的特点 16
6.3.1 RNA-seq对芯片的优势 16
6.3.2 RNA-seq存在的问题 17
6.4 RNA-seq数据预处理 18
6.4.1 质量控制 18
6.4.2 读段清理 22
6.4.3 转录组组装 25
6.4.4 转录组定量和标准化 25
6.4.5 线性相关系数 27
6.5 RNA-seq数据分析 28
6.5.1基因表达差异的显著性分析 28
6.5.2 RNA-seq数据的其它分析 31
参考文献: 31
第七章 R的高级语法与如何创建R包
7.1 R的高级语法 2
7.1.1 数据类型及相互转换 2
7.1.2 向量运算 6
7.1.3 函数 8
7.1.4 循环与条件 9
7.1.5 输入输出 11
7.1.6 对象和类 12
7.2 创建及发布自己的R/Bioconductor包 14
7.2.1 在Windows下创建和发布R包 14
7.2.2 在Linux下创建和发布包 25
7.3 R包结构 26
7.3.1 R的源代码包 27
7.3.2 R的二进制包 28
参考文献: 29

附录A 进一步学习的资源
附录B R常用函数
附录C R的内存管理和帮助系统

购买方式:

支付宝账号: [email protected] 何静
留言:姓名+R语言书+册数
定书后,请将邮寄地址(一定要邮编)和电话发以上电子邮件
电子邮件题目是买R语言书,内容就是姓名+邮寄地址+电话
本书定价58元,群内读者预定50元(包快递费);团购5本以上45元(包快递费);50本40元(包快递费);
如果大学指定为教材,我可以提供更多材料,长期合作,打造精品课程

一份数据挖掘小书评

书看得多了,渐渐有了品味,所以看那些英文版的时候也不再有很膜拜的感觉。不过再差的书,总会有些许长处或教训,值得借鉴。下面列出一些在数据挖掘实现方面的书,主要是基于R和python的。

  • Data mining with R Learning with Case studies: 我的启蒙书,很好的案例教学,这种书应该出更多些。
  • Data Mining Applications with R: 风格类似上面,多人合作的案例讲解,非常不错,就是每个例子略为简短了一点。
  • Machine Learning for Hackers: 也是以案例为主,但是内容质量不能保持在同一水准上,代码讲的比较细。
  • An Introduction to Statistical Learning: 比较系统的讲解了统计学习的内容,以及相应的R包函数,优点在于有习题。
  • Machine Learning with R: 同样是系统的列出各机器学习算法对应的R包函数,直观理解很好,而且讲了一些caret包的使用。比较新的书。
  • Applied predictive modeling: 对预测算法讲的最细,涉及R包最多,由于作者是caret的开发者,所以要学习caret则必看此书。
  • Big Data Analytics with R and Hadoop: 刚出的大数据挖掘书籍,讲解了如何用RHadoop来整合R和hadoop,什么时候会有R和spark的书出来啊。

  • programming_collective_intelligence: 有些年头的书了,但确实是经典,原理讲的简洁,就讲python代码实现,喜欢这种风格。
  • Machine Learning in Action: 和上面风格类似,不过实现是主要基于python的numpy库,这样代码量要精简很多,也很好。
  • Python for data analysis: 主要讲pandas库,精彩例子不多,可以直接去看帮助文档。
  • Learning scikit-learn Machine Learning: 主要以scikit-learn扩展库为工具做机器学习,比较简单,可以直接去看帮助文档。
  • Building Machine Learning Systems with Python: 和上面那本类似,讲函数功能为主,原理直觉没涉及到。

用模拟来理解混合效应模型之一:random intercept model

混合效应模型(Mixed-effect Model)在之前文章提到过,但感觉仍是雾里看花。此番又研究了一些资料,准备来做一个系列讲讲清楚,也算是自己学习的一个总结。

普通的线性回归只包含两项影响因素,即固定效应(fixed-effect)和噪声(noise)。噪声是我们模型中没有考虑的随机因素。而固定效应是那些可预测因素,而且能完整的划分总体。例如模型中的性别变量,我们清楚只有两种性别,而且理解这种变量的变化对结果的影响。

那么为什么需要 Mixed-effect Model?因为有些现实的复杂数据是普通线性回归是处理不了的。例如我们对一些人群进行重复测量,此时存在两种随机因素会影响模型,一种是对某个人重复测试而形成的随机噪声,另一种是因为人和人不同而形成的随机效应(random effect)。如果将一个人的测量数据看作一个组,随机因素就包括了组内随机因素(noise)和组间随机因素(random effect)。这种嵌套的随机因素结构违反了普通线性回归的假设条件。

你可能会把人员ID(组间的随机效应)看作是一种分类变量放到普通线性回归模型中,但这样作是得不偿失的。有可能这个factor的level很多,可能会用去很多自由度。更重要的是,这样作没什么意义。因为人员ID和性别不一样,我们不清楚它的意义,而且它也不能完整的划分总体。也就是说样本数据中的路人甲,路人乙不能完全代表总体的人员ID。因为它是随机的,我们并不关心它的作用,只是因为它会影响到模型,所以不得不考虑它。因此对于随机效应我们只估计其方差,不估计其回归系数。

混合模型中包括了固定效应和随机效应,而随机效应有两种方式来影响模型,一种是对截距影响,一种是对某个固定效应的斜率影响。前者称为 Random intercept model,后者称为 Random Intercept and Slope Model。Random intercept model的函数结构如下:

Yij = a0 + a1*Xij +  bi + eij

a0: 固定截距
a1: 固定斜率
b: 随机效应(只影响截距)
X: 固定效应
e: 噪声

下在的代码就是人工生成这种数据,再用nlme包建模展示。

参考资料:
Regress by simulation
A very basic tutorial for performing linear mixed effects analyses


如何善用R的语法:三段示例

在改进R代码的运算效率时,我们常说的第一句话就是:要善用R的语法。那么问题是:HOW?下面用一个具体的问题来表达这句比较抽象的话。问题就是用R代码来模拟一个二维布朗运动,一个人从原点出发,随机的向任何一个方向前进或后退一步,不会停留在原地。这个问题是一位英国留学的朋友问我的,我觉得可能比较典型,就顺便写在博客上供大家参考比较吧。

下面就这个问题,给出了三个解答,slowwalk是运行最慢的,在我的电脑上大约要运算8秒,这段代码慢的原因在于用c来补进新的向量元素。而在R中,对一个向量进行新增的操作实际上是进行了复制重建,所以速度是很慢的。为了改进这一点,speedwalk1函数中,先用一个空向量占好了位置,再向这个向量中放东西,这样不需要反复新建,速度会增加10倍,只花了0.7秒左右。但R最重要的运算思路是向量化,也就是同时对向量各元素进行计算,而不是把向量元素一个个放到循环中。speedwalk2函数就是向量化计算的实现,这样下来代码速度又加快了十几倍,大约只需0.05秒,而且还没有用到多核心计算。

奇异值分解和潜在语义分析

奇异值分解是一种矩阵分解技术,如果将列看作特征,SVD也能看作是一种降维技术。如果将数据先进行标准化之后,再进行SVD操作,其左奇异向量U就是主成分得分,右奇异向量V就是主成分负载。但SVD的优势恰恰在于不进行标准化,当我们处理高维稀疏数据(例如文本数据)时,标准化会使矩阵不再稀疏,影响计算效率。所以一般是用SVD直接对原始数据进行操作。得到U的第一列是平均值的意义,U的后几列才有区别各样本的功能。

在文本挖掘中,潜在语义分析(LSA)假设在词项和文档之间有一个潜在语义层,词汇和语义产生关联,文档也和语义产生关联。这样通过语义概念,可以将词项和文档放在一个语义空间中观察分析。进一步我们可以将语义所为词项的“主成分”放到数据挖掘中进行分析。LSA的一般过程是先构造词项-文档矩阵(其转置为文档-词项矩阵,视需要而定),再计算TFIDF值,然后用SVD来分解这个矩阵,得到的奇异向量就是潜在语义。U的第一列可以认为是词项在文档中的平均出现频率,V的第一列可以认为是文档中包含的词项频率,我们更有兴趣的是后几列。关于SVD和LSA有一份精彩的入门文档,下面的R代码就是对该文档的模仿,但是SVD得到的结果却和文档不同,我很奇怪为什么会这样。(使用的数据

说英雄谁是英雄之光荣三国

本文的缘起是知乎的这一篇文章,老马这位三国历史专家让大家都知道了有这么一份三国志游戏数据的存在。这份数据很容易从网上找到,而本人只关心三国志四的各人物数据,为什么呢,因为只玩过三国志四嘛。游戏是多年前玩的,如今来玩玩数据吧。下面的问题主要围绕一个焦点,谁是三国志四中文武兼资第一人。

网上下载的数据是一个Excel表格,整理之后只取了三国志四的部分,去除了分值为0的人物,整理后数据你可以从这里下载。整个数据有454名人物,包括了六个变量,分别为人物’名字’,’统率’,’武力’,’智力’,’政治’,’魅力’。我们可以使用plyr包中的arrange函数找到各指标分别对应的前十名。不出意外,武神关二爷获得统率第一,人中吕布获得武力第一,诸葛亮获得智力第一,张昭获得政治第一,而会哭的刘备获得魅力第一,晕!

arrange(data,desc(tongshuai))[1:10,]
arrange(data,desc(wuli))[1:10,]
arrange(data,desc(zhili))[1:10,]
arrange(data,desc(zhengzhi))[1:10,]
arrange(data,desc(meili))[1:10,]
 
由于每个人物对应多个指标,我们需要一个综合指标来比较,最简单的,上PCA。
prinModel <- prcomp(x=data[,2:6],
retx=TRUE,
center=TRUE,
scale=TRUE)
 
prinModel$rotation[,1:2,drop=F]
summary(prinModel)
 
对数据进行主成份分析,发现前两个主成份的累积方差贡献率达到88%,还不错,第一PC以智力和政治为主,达到50%贡献率,第二PC以武力和统率为主,达到30%贡献率。那么第一主成份可以认为是文韬因素,第二主成份可以认为是武略因素。好的,下面可以根据两个主成份得分列出文韬前十和武略。

# 武略榜
data$name[order(prinModel$x[,2],decreasing=T)][1:10]
[1] "關羽" "趙雲" "曹操" "孫策" "孫堅" "黃忠"
[7] "周瑜" "陸遜" "夏侯惇" "姜維"
# 文韬榜
data$name[order(prinModel$x[,1],decreasing=T)][1:10]
[1] "郭嘉" "張紘" "張昭" "諸葛亮" "荀彧" "諸葛瑾"
[7] "荀攸" "曹植" "魯肅" "馬良"
 
下面要将这两个主成分合为一个综合指标,采用加权平均方法,权重是方差贡献率,大致上0.6*文韬+0.4*武略,下面就取出兼资文武之名将前十名,当当当当,他们是:
final <- prinModel$x[,1]*0.6 + prinModel$x[,2]*0.4
data$name[order(final,decreasing=T)][1:10]
 
[1] "諸葛亮" "曹操" "周瑜" "陸遜" "司馬懿" "趙雲"
[7] "魯肅" "姜維" "龐統" "徐庶"
 
看来在三国志的游戏数据中,诸葛亮是最强人物啊。最后来看看蜀国五虎排序,通常民间排序是关张赵马黄,这种序列显然没有考虑到智力因素问题啊,下面我们用前面的评价体系来为五虎重新排序。
shu5 <- c('趙雲','馬超','張飛','關羽','黃忠')
data.frame(round(prinModel$x[data$name %in% shu5,1:2],2),
final=round(final[data$name %in% shu5],2),
name=data$name[data$name%in% shu5])
 
PC1 PC2 final name
437 -1.51 2.22 -0.02 馬超
488 -2.65 1.31 -1.07 張飛
658 0.59 2.92 1.52 黃忠
741 1.68 3.31 2.34 趙雲
922 1.20 3.48 2.12 關羽
 
从上面这个表格可以看到,根据游戏数据正确的排序应该是赵关黄马张,赵云君,我看好你啊。我们还可以看看五虎将在群雄中的位置,将PC1和PC2做散点图如下,这张图中处于右上角的就是文韬武略兼备的人才。
library(ggplot2)
ggdata <- data.frame(prinModel$x[,1:2],name=data$name)
p <- ggplot(ggdata,aes(PC1,PC2))+
geom_point(size=3,alpha=0.5)+
geom_text(data=subset(ggdata,name %in% shu5),
aes(label=name))+
labs(x='文韬',y='武略')+
theme_bw()
print(p)


 
最后还可以用雷达图看这五虎的技能值


使用PCA来进行综合指标评价不一定是最好的方法,各位有什么好主意,请不吝赐教。

生存分析函数小结

生存分析(survival analysis)适合于处理时间-事件数据。例如中风病人从首次发病到两次复发,其中就涉及到时间和事件。此例中时间就是复发的时间间隔,事件就是是否复发。如果用普通的线性回归对复发时间进行分析,就需要去除那些没有复发的病人样本。如果用Logistic回归对是否复发进行分析,就没有用到时间这个因素。而生存分析同时考虑时间和事情这两个因素,效果会更好些。

在R语言中我们可以使用survival包进行生存分析,其中主要的函数功能罗列如下:

Surv:用于创建生存数据对象
survfit:创建KM生存曲线或是Cox调整生存曲线
survdiff:用于不同组的统计检验
coxph:构建COX回归模型
cox.zph:检验PH假设是否成立
survreg:构建参数模型


下面是使用一个实例来使用R中的生存分析函数,其中用到的数据集可以在这里下载