Illumina公司有一张特别牛x的poster，叫做 For all you seq，全面记录了目前问世的所有二代测序技术的建库protocol，但是这张poster是加密的，无法分割、提取进行打印。然后，他们又出了一个手册，来介绍For all you seq里的建库方法，叫做Sequencing Methods Review，全手册加上illumina的广告将近200页。我打算从头到尾都看一遍，有兴趣的就写出来分享一下。

对于做生物信息的人，为啥要知道实验上的建库原理？

因为测序数据的第一道处理步骤就是排除并分析建库和测序时产生的错误和误差。为什么会产生这些问题？如何在mapping时排除相应问题，做到准确alignment？

第一部分就从RNA转录开始说起，RNA转录的度量方法分为很多种：

• post-translational modifications 后转录修饰
• RNA splicing 剪接
• RNA bound to RNA binding proteins (RBP) 与RNA相关蛋白结合
• RNA expressed at various stages 在不同阶段的表达量
• unique RNA isoforms 转录本异构体
• RNA degradation 降解
• regulation of other RNA species 调控其他RNA “物种” （我不清楚这句话的意思）

Chromatin isolation by RNA purification

Chromatin isolation by RNA purification (ChIRP-Seq) RNA纯化染色质分离技术是用于定位非编码RNA，例如长非编码RNA及其蛋白同染色质相互作用的位置。

首先进行交联，超声打断，用Biotinylated寡核苷酸对感兴趣的RNA进行杂交，并用链霉抗生物素蛋白磁珠捕获复合物。用RNase H（核糖核酸酶H）处理后，对DNA进行提取和测序。

技术优点：

1. 可在基因组全部区域找寻结合位点
2. 可发现新结合位点
3. 可只对感兴趣的RNA进行相应研究

技术缺点：

1. 非特异性的寡核苷酸相互作用会导致对结合位点的错误解读
2. 染色质会在预处理阶段降解
3. 感兴趣的RNA序列信息必须已知 （这不是废话吗，感兴趣的RNA肯定是已知的，已知的能不知道序列？感兴趣的未知RNA的作用位点，知道了也没法研究啊）

Accumulation Distribution 乍一看可以翻译成累积分布，但累积分布的正确英文写法是cumulative distribution。 今天要说的这个Accumulation Distribution是证券交易中的技术指标——累积/派发线。

难道我这个浓眉大眼的家伙也开始炒股了？ -__,- 当然不是，我只是想用空闲时间了解一些经济学基本知识，从网上找到了Online Trading Concepts ，这个和计算有点关系的技术指标介绍网站，所以打算好好看一看，为了记录自己的学习成果，有机会就编译点网站上内容。

累积/派发线通过成交量来确认价格趋势，对于可能导致价格波动反转的弱趋势给出预警。

累积：当天收盘高于前一天收盘价时，将成交量计入前一日的成交量累积值中，作为今日的累积/派发线值。 派发：当天收盘低于前一天收盘价时，将成交量从前一日的成交量累积值中减去，作为今日的累积/派发线值。

所以累积/派发线是用来检测价格矩估计和成交量矩估计之间的差异。

两个时间点之间的成交量曲线是上涨的，那么说明价格会上涨，泛指则应认定为价格会有波动（下降）。

• High #1 到 High #2

图中的High #1 到 High #2 的价格是一样的，但是累积/派发线呈下降趋势，说明价格会有波动（下降）。

• High #3 到 High #4

图中的High #3 到 High #4 的价格是上涨的，但是累积/派发线呈下降趋势，说明价格会有波动（下降）。

• Low #1 to Low #2

a sign of strength， 翻译不好这个strength的意思。但是看图，应该还是波动（下降）趋势。

累积/派发线是一个实用的检测价格趋势反转的指标。其他类似的分析成交量和价格的指标是Chaikin Oscillator (蔡金摆动指标), Money Flow Index (资金流量指标)，以及 Price Volume Trend indicator (价格交易量趋势指标)。

R 代码 Williams Accumulation / Distribution

library(TTR)
library(magrittr)
library(plotly)
data(ttrc)
ttrc$ad <- williamsAD(ttrc[,c("High","Low","Close")])
df=tail(ttrc,30)
i <- list(line = list(color = '#FFD700'))
d <- list(line = list(color = '#0000ff'))

p <- df %>%
  plot_ly(x = ~Date, type="candlestick",
          open = ~Open, close = ~Close,
          high = ~High, low = ~Low,
          increasing = i, decreasing = d)%>%
  add_lines(y= ~ad, name = "williamsAD")  %>%
  layout(yaxis = list(title = "Price"))

水一篇。我喜欢这种类似回文串的日期。 有三个多月没更新博客了，马上满血复活，要填的坑还是比较多，希望能早日填完。

坑之列表

近期

• animation package bug
• Freshman21 theme bug
• 欠的blog文章 (10+), see todo list in my driver
• 某bot （五月初完成一个bot，但还差一个-_-）

远期

• 某R包编写

久远

• 某专题

其余部分包括Git的使用，可重复性研究（Rmarkdown），以及一些简单的生物信息学分析（Blast，DNA mutation，RNA expression），内容不连贯，我就不一一罗列了。

我当时有个笔记草稿，现已放到了gist上，欢迎参观，但应该没人能完全看懂。。。

2017年更新

最新知识，比较RNA表达量的分布，需要知道最少用几个样本，有了样本才能知道表达量的分布，那么究竟表达量分布长什么样子，需要多少个样本？

在R语言部分，老师推荐了这本书An Introduction to Statistical Learning with Applications in R。

还稍微介绍了一点编译相关的内容。

接下来就是简单的实战，我发现没有编成经验的学生最容易弄错的一段代码是随机数，他们一般都会问，为什么自己的结果同老师演示的结果不同。

x=100
myfunction1<-function(x){
  v=rnorm(x)
  mean(v)
}
myfunction1(x)

全局变量

（20170410更新，我觉得文章这么长，你们肯定不会看到最后的留言）

在课堂中老师忘了如何将函数中的变量定义为全局变量，其实这个功能我在R中基本就没有用到过。不过我还是查到了方法，用双箭头<<-来定义全局变量。

谢老大留言：双箭头 «- 表示全局变量这个说法不严谨，但很多人都被误导了。其实它的精确含义是向上层环境中的变量赋值，问题只是如果上层变量中不存在这个变量的话，就继续再向上，如果一直都找不到这个变量，最上层的环境是全局环境 .GlobalEnv，所以有时候就变成了创建全局变量。我自己经常用双箭头魔法向上层环境赋值，比如 knitr 中的 opts_chunk$set()就是基于双箭头的：https://github.com/yihui/knitr/blob/master/R/defaults.R

具体问题参考Define global variable using function argument in R

x=100

myfunction=function(x){
  v<<-rnorm(x)
  mean(v)
}
myfunction(x)
v

杀掉Rstudio进程

有时由于程序写错了，进入了死循环，并且内存占用越来越大，最简单的处理方式是在终端中杀掉Rstudio进程（虽然不推荐）。

killall rstudio

R的数据结构

• vector
• list
• data.frame （还有一种data.table）
• matrix
• character

课上练习

words=c("Hello","and hi","World")
l=list()
for(i in 1:length(words)){
  if(i==3 && (words[i]=="World" || words[i]=="world")){
    l[[i]]="America!"
  }else{
    l[[i]]=words[i]
  }
}
for(i in l){
  cat(i," ")
}

实用网站和关键函数介绍

• rseek 查询
• R操作算子顺序
• R逻辑操作符
• R句法

正则表达式

没教\w，老师的观念是用\S,\s,\d,\d+来推导一切表达式。

fh=file("/media/sf_CB2/dm/dm.faa",open="r")
while (length (line=readLines(fh,warn=F,n=1))>0){
  pattern="^\\>gi\\|(\\d+)\\|ref\\|(\\S+)\\|"
  #m=regexec(pattern,line,perl=T)
  #m=mymatch(pattern,line)
  #if(m[[1]][1]==-1){
  #  next
  #}
  #if(!(ifmatched(m))){
  #  next
  #}
  m=mymatch(commapattern,line)
  if(!ifmatched(m)){
    my=mymatch(pattern,line)
  }
  if(!(ifmatched(m))){
    next
  }
  cat(getstring(line,m,2))
  cat("\n")
  #v=regmatches(line,m)
  #vec=v[[1]]
  #cat(vec[2],vec[3],sep="\t")
  #cat("\n")
}
close(fh)

pattern="\\|\\s+(.\*)\\s+\\["

commapattern="\\|\\s+(.\*)\\,\\s+\\["

getstring=function(text,match,index){
  v=regmatches(text,match)
  vec=v[[1]]
  return(vex[[index]])
}

mymatch=function(pattern,text){
  m=regexec(patter,text,perl=T)
  return(m)
}

ifmatched=function(m){
  if(m[[1]][1]==-1){
    return(FALSE)
  }else{
    return(TRUE)
  }
}

v=rnorm(100)
median(v)
mean(v)
hist(v)

R画图

用最简单的图形命令来画图，在画图之前还穿插了Git的介绍（我会在后续文章中介绍）。

data("airquality")
mean(airquality$Ozone)
mean(airquality$Ozone[is.na(airquality$Ozone)])

d=airquality$Ozone
index=is.na(airquality$Ozone)
clean_idx=!index
clean_value=d[clean_index]
mean(clean_value)
mean(airquality$Ozone,na.rm=T)
summary(airquality)


plot(airquality$Temp,airquality$Ozone, xlab="Temp", ylab="", main="Temp vs Ozone")

plot(airquality$Temp,airquality$Ozone,type="n",xlab="",ylab="",axes=F)
points(airquality$Temp,airquality$Ozone,pch=16)
axis(1)
axis(2)
box()
title(main="Temp vs Ozone",xlab="Temp",ylab="Ozone")
legend("topleft",legend="Some points",pch=16)

par(mfrow=c(1,2))
plot(airquality$Temp,airquality$Ozone, xlab="Temp", ylab="", main="Temp vs Ozone")
plot(airquality$Temp,airquality$Ozone, xlab="Temp", ylab="", main="Temp vs Ozone")