One Tip Per Day

fgrep MAN Page

2012-01-05T08:12:00.000-08:00

fgrep MAN Page: "fgrep

Search file(s) for lines that match a fixed string

Syntax
fgrep ...

fgrep is the same as `grep -F'

all other options are the same as grep
It's a popular fallacy that fgrep stands for fast-grep, in fact fgrep means fixed-string grep.

fgrep implements the Aho–Corasick string matching algorithm which is very fast at matching multiple strings in the input stream/file. So if fgrep is used like this:

fgrep -f patternlist.txt largetextfile.txt

it is much faster than

grep -f patternlist.txt largetextfile.txt

Matching a single pattern with fgrep is not significantly faster than grep.

"The knowledge imposes a pattern, and falsifies, For the pattern is new in every moment..." - T. S. Eliot"

'via Blog this'

an easy way to writing data.frame to Excel

2011-11-16T08:13:00.001-08:00

you can write it as

write.table(r.data.frame, "excel.file.xls", sep="\t", na="", row.names=F)

which I can usually open in Excel just by clicking on it.

Credit: http://tolstoy.newcastle.edu.au/R/help/05/04/3388.html

match vs. %in%

2011-11-04T12:59:00.000-07:00

match and %in% are two very commonly-used function in R. So, what's the difference of them?

First, how to use them -- (copy from R manual)

match returns a vector of the positions of (first) matches of its first argument in its second.
%in% is a more intuitive interface as a binary operator, which returns a logical vector indicating if there is a match or not for its left operand.

match(x, table, nomatch = NA_integer_, incomparables = NULL)x %in% table

Examples:

> a
[1] 1 1 0 1 5 1 2 4
> b
 [1] 10  9  8  7  6  5  4  3  2  1
> match(a,b)
[1] 10 10 NA 10  6 10  9  7
> a %in% b
[1]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE

So, if two vectors are overlapped like

a ---------------

b -----------------------------

To get the overlapped part in order of a, use a[a %in% b], even though there are duplicates in the overlapped part. However, this does not work for match, since match() only returns the first match of a in b. For example,

> match(b,a)
 [1] NA NA NA NA NA  5  8 NA  7  1
> match(b,a, nomatch=0)
 [1] 0 0 0 0 0 5 8 0 7 1
> a[match(b,a, nomatch=0)]
[1] 5 4 2 1

even using 'nomatch=0', the final command still returns 4 elements, not the overlapped ones.

用粤语数数

2011-10-25T13:27:00.001-07:00

小曦最近在中文学校教kindergarten的小朋友们中文。根据教材的设置，我们是从中文的数字开始学起的。学习一种新的语言，不是从吃喝拉撒衣住行这种最生活的东西开始，而是从数数开始，这确实让人说不出道理却又总习惯去做。我很小的时候就开始学习英语，老人见了我，问我英语学得怎么样，让我表演表演，我开口就是“汪，土，水，火……”。而到了挪威，想跟朋友学学挪威语，也是劈头就问人家怎么用挪威语数一到十。粤语也是一样，我看到很多人都说自己会一些粤语，让他们说说看的时候，他们就开始“呀，咦，三，噻……”。我自己也一样，认为把粤语里的一到十都说得很溜了，后来认真学习粤语拼音才发现，粤语的发音可是有玄机的。粤语不但声调比普通话多，而且还有很多带韵尾的字，而韵尾的概念是普通话里根本没有的。就拿数字一到十举例，除了“二”、“四”、“九”这三个字的发音跟普通话类似——也是“声母+韵母”的模式之外，其余几个都是很特别的。以下是我的一些心得，不一定完全准确。

“一”，拼作“jat1”，格式可以解读为声母“j”+韵母“a”+韵尾“t”，调号为“1”；因为有了韵尾t，所以收尾短促，并且以舌顶上齿、欲发“t”音但尚未发出的状态作结。因此，读“一”的时候，要短促而有力，而不是很多人读的“呀aaaaaa”；

“二”，拼作“ji6”，可以读得拖一些；

“三”，拼作“sam1”，收尾的时候要闭嘴，做出发“m”音的状态；

“四”，拼作“sei3”；

“五”，拼作“ng5”，这个音最有意思，因为发这个音的时候，嘴从始至终都没有张开过，而是从鼻子里哼出一个调子，类似于普通话里的“嗯”。粤语里常用“唔”取代“不”，这个“唔”也是和“五”类似的发音（拼作“ng4”）。比如“不用”在粤语里是“唔使”，我最喜欢和大黄说这个词，说的时候，就是鼻子里憋一口气先哼一声读出“唔”，再顺势张开嘴读出“使”；但大黄说的时候，就会说成“毛塞”，就是读“唔”的时候没有把嘴巴闭起来的原因。

“六”，拼作“luk6”，这个音带着韵尾“k”，和“一”一样，也是一个入声字。虽然“k”不真的发出来，但是已经有这个趋势了。如果这个字后面再跟上一个字一起读的话，有些时候就能听出韵尾“k”了。比如，英语“taxi”，粤语翻译成“的士”，这个音译就比普通话的音译“出租车”要更接近原来语言的发音，因为“的”（dik1），当它和“士”（si3）拼在一起的时候，就读作“diksi”，是可以听出“k”音的！

“七”，拼作“cat1”，参看“一”；

“八”，拼作“baat”，这个韵母“aa”比“a”的嘴型要大些，韵尾的读法同上；

“九”，拼作“gau”；

“十”，拼作“sap”，韵尾是“p”，收尾的时候又是要闭嘴的。我想起一部港片里，一个内地人甲跟香港人乙赔礼道歉后说了句，“这样的话，你接不接受啊？”结果那个香港人一听就火了，为什么呢？因为甲不懂韵尾，他把“接”按照普通话的习惯读成了“zi”，而不是“zip”！“zip sau”和“zi sau”，仅仅一“p”之差，前者是“接受”，后者听起来就像是“折寿”了。谁听了不生气！

Sites about info visualization

2011-10-24T12:46:00.000-07:00

It's not all, definitely. But these are sites I found very interesting to read and helpful to learn:

http://visualeconomics.creditloan.com/ (what I newly found today)
http://www.datavizchallenge.org/ (a contest site for Vis)
http://vis.stanford.edu/ (A pioneer group in this study)
http://mbostock.github.com/d3/ (one of library from above group)
http://flowingdata.com/ (Nathan's website)
http://www.google.com/publicdata/home (Google Public Data)
http://www.gapminder.org/ (where Prof. Rosling's original TED talk mentioned)
http://periscopic.com (A company in this field)
adding...

let me know if you know any else.

如何把文字排版成异形

2011-10-10T16:35:00.000-07:00

如何把一篇文字排版成喜欢的形状呢？PS的路径工具可以帮你！

是不是有点儿意思？其实做法很简单。以下是步骤：1、下载一个心形的图片；2、在ps里生成心形的路径（可以用路径铅笔手动勾勒出路径，也可以用魔术棒创建选区之后生成路径）；3、在生成的路径下使用文字工具T，写出预设的文字段落；4、调整字体类型、大小和颜色；搞定！

再做一个，这次用sin wave旋转之后来做，文字是最近喜欢的一首悲情林教头的歌——《自从你走后》

data frame extraction

2011-10-10T12:03:00.000-07:00

This is a trivial but very useful tip:

> x=data.frame(a=1:4, c=5)
> x
  a c
1 1 5
2 2 5
3 3 5
4 4 5
> x[1,] 
  a c
1 1 5
> x[,1]
[1] 1 2 3 4
> x[,1, drop=F]
  a
1 1
2 2
3 3
4 4

where you can see that:

to avoid a[, i] become a vector, rather than a one-column data frame, we can use drop=F option.

Here is full explanation for that:

http://stat.ethz.ch/R-manual/R-patched/library/base/html/Extract.data.frame.html

k-mean clustering + heatmap

2011-10-10T11:07:00.000-07:00

Several R functions in this topic:

1. dist(X) -- calculate the distance of rows of data matrix X. The default distance method is euclidean. It can be maximal, manhattan, binary etc.

> a=matrix(sample(9),nrow=3)
> a
     [,1] [,2] [,3]
[1,]    5    2    9
[2,]    8    7    1
[3,]    6    4    3
> dist(a, diag=T, method='max')
  1 2 3
1 0    
2 8 0  
3 6 3 0
 
> dist(a, diag=T, method='euc')
         1        2        3
1 0.000000                  
2 9.899495 0.000000         
3 6.403124 4.123106 0.000000

2. hclust(D) -- hierarchical clustering of a distance/dissimilarity matrix (e.g output of dist function): join two most similar objects (based on similarity method) each time until there is one single cluster.

hclust(D) can be displayed in a tree format, using plot(hclust(D)), or plclust(hclust(D))

3. heatmap(X, distfun = dist, hclustfun = hclust, ...) -- display matrix of X and cluster rows/columns by distance and clustering method.

One enhanced version is heatmap.2, which has more functions. For example, you can use

key, symkey etc. for legend,
"col=heat.colors(16)" or "col='greenred', breaks=16" to specify colors of image
cellnote (text matrix with same dim), notecex, notecol for text in grid
colsep/rowsep to define blocks of separation, e.g. colsep=c(1,3,6,8) will display a white separator at columns of 1, 3, 6, 8 etc.

Both have 'ColSideColors/RowSideColors', a color vector with length of cols/rows. Here is an example(http://chromium.liacs.nl/R_users/20060207/Renee_graphs_and_others.pdf).

Another enhanced version is pheatmap, which produced pretty heatmap with additional options:

cellwidth/cellheight to set the size of cell
treeheight_row/treeheight_col: height of tree
annotation: a data.frame, each column is an annotation of columns of X. So, nrow(annotation)==ncol(X)
legend/annotation_legend: whether to show legend
filename: save to file

4. kmeans(X, centers=k) -- partition points (actually rows of X matrix) into k clusters . For example:

# a 2-dimensional example
x <- rbind(matrix(rnorm(100, sd = 0.3), ncol = 2),
matrix(rnorm(100, mean = 1, sd = 0.3), ncol = 2))
colnames(x) <- c("x", "y")
(cl <- kmeans(x, 2))
plot(x, col = cl$cluster)
points(cl$centers, col = 1:2, pch = 8, cex=2)

The number of cluster can be determined by plot of sum of squares, eg.

# Determine number of clusters
wss <- (nrow(x)-1)*sum(apply(x,2,var))
for (i in 2:20) wss[i] <- sum(kmeans(x,centers=i)$withinss)
plot(1:20, wss, type="b", xlab="Number of Clusters",ylab="Within groups sum of squares")

Using hclust and cutree can also set the number of clusters:

hc <- hclust(dist(x), "ward")
plot(hc) # the plot can also help to decide the # of clusters
memb <- cutree(hc, k = 2)

Note: kmean is using partition method to cluster, while hclust is to use hierarchical clustering method. Here is a series of nice lectures for this. A more detail for cluster can be found here: CRAN Task View: Cluster Analysis

about methylation （甲基化）

2011-09-27T13:43:00.000-07:00

'Methylation' has been more and more widely used in the field of gene regulation. Here is what I cracked:

Methylation means addition of a methy- group to a subtract (see the right figure).
Methylation can mean DNA methylation or protein methylation.
DNA methylation occurs in a CpG (C followed by G) site in DNA sequences, where cytosine (C) becomes 5-methylcytosine (see right figure, counting counterclockwise from the NH).
About 80-90% of CpG sites in human DNA are methylated.
Most CpG islands (i.e. CG% >65%) are un-methylated.
Around 56% of mammalian genes has (un-methylated) CpG island in their promoter, like most housekeeping genes.
In normal tissue, methylation of a gene is mainly localized to the coding region, which is CpG-poor. In contrast, the promoter region of the gene is unmethylated, despite a high density of CpG islands in the region.
Methyl-seq is a technique combining DNA digestion by a methyl-sensitive enzyme HpaII and its methyl-insensitive isoschizomer MspI with the Illumina DNA sequencing platform. Sequence tags present in MspI libraries but not in HpaII libraries are derived from methylated regions. Conversely, sequence tags that occur in HpaII libraries come from at least partially unmethylated regions.
HpaII only digests 5'-CCGG-3' sites when the cytosine in the central CG dinucleotide is unmethylated, so the HpaII representation is enriched for the hypomethylated fraction of the genome. The MspI representation is a control for copy number changes and PCR amplification difficulties. Refer to HELP assay wiki page.
Reduced Representation Bisulfite Sequencing (RRBS) is another method for DNA methylation. It has base-pair resolution. More detail can be refer her: http://genome-test.cse.ucsc.edu/cgi-bin/hgc?hgsid=3051149&o=181430182&t=181430183&g=wgEncodeHaibMethylRrbsK562HaibSitesRep1&i=K562_Rep3_RRBS

DTW: dynamic time warping 动态时间规整

2011-09-18T17:34:00.000-07:00

Basically, DTW (dynamic time warping) is an algorithm to output cumulative distance of two time sequences, which is widely used e.g. for classification and clustering.

For example, when using k-mean for clustering, we can use DTW as distance function. Here is one of such nice instances (using R: http://www.rdatamining.com/examples/ts-mining)

Relevant information from Anshul's email.

A review of DTW http://csdl.ics.hawaii.edu/techreports/08-04/08-04.pdf

Code:

Python code: https://mlpy.fbk.eu/
R code: http://dtw.r-forge.r-project.org/
JAVA code: http://java-ml.sourceforge.net/ or http://code.google.com/p/fastdtw/
MATLAB: the samplealign() function in the bioinformatics toolbox does DTW

There is also the global alignment kernel that is faster and more accurate that DTW. It can be used to compute distance between time series to be used in clustering allowing for all possible global alignments. Here is code http://www.iip.ist.i.kyoto-u.ac.jp/member/cuturi/GA.html

Here is an example (in Chinese) from Ckary's blog:

在日常的生活中我们最经常使用的距离毫无疑问应该是欧式距离，但是对于一些特殊情况，欧氏距离存在着其很明显的缺陷，比如说时间序列，举个比较简单的例子，序列A：1,1,1,10,2,3，序列B：1,1,1,2,10,3，如果用欧氏距离，也就是distance[i][j]=(b[j]-a[i])*(b[j]-a[i])来计算的话，总的距离和应该是128，应该说这个距离是非常大的，而实际上这个序列的图像是十分相似的，这种情况下就有人开始考虑寻找新的时间序列距离的计算方法，然后提出了DTW算法，这种方法在语音识别，机器学习方便有着很重要的作用。这个算法是基于动态规划（DP）的思想，解决了发音长短不一的模板匹配问题，简单来说，就是通过构建一个邻接矩阵，寻找最短路径和。还以上面的2个序列作为例子，A中的10和B中的2对应以及A中的2和B中的10对应的时候，distance[3]以及distance[4]肯定是非常大的，这就直接导致了最后距离和的膨胀，这种时候，我们需要来调整下时间序列，如果我们让A中的10和B中的10 对应，A中的1和B中的2对应，那么最后的距离和就将大大缩短，这种方式可以看做是一种时间扭曲，看到这里的时候，我相信应该会有人提出来，为什么不能使用A中的2与B中的2对应的问题，那样的话距离和肯定是0了啊，距离应该是最小的吧，但这种情况是不允许的，因为A中的10是发生在2的前面，而B中的2则发生在10的前面，如果对应方式交叉的话会导致时间上的混乱，不符合因果关系。接下来，以output[6][6](所有的记录下标从1开始，开始的时候全部置0)记录A，B之间的DTW距离，简单的介绍一下具体的算法，这个算法其实就是一个简单的DP，状态转移公式是output[i][j]=Min(Min(output[i-1][j],output[i][j-1]),output[i-1][j-1])+distance[i][j];最后得到的output[5][5]就是我们所需要的DTW距离.

The C code there is also much helpful in understanding the algorithm.

Mozaika: a nice tool to make mosaic

2011-09-17T19:42:00.000-07:00

If you want to make a photographic mosaic, Mazaika would be one of the solutions (at least from what I knew). Here is the website for Mac version: http://www.mazaika.com/mac/

Here is what I tested on a profile photo I took when travelling in Daocheng.

(click to see the bigger photo)

ifelse function in R only returns the first element

2011-09-16T14:50:00.000-07:00

If you also favor to use the function, be aware of the returned value. For example:

> ifelse(1>0, 3, 4) 
[1] 3
> ifelse(1>0, c(2, 3), c(4, 5))   # only the first element returned.
[1] 2
 
> ifelse(c(1:10)>5, 'on', 'off')
 [1] "off" "off" "off" "off" "off" "on"  "on"  "on"  "on"  "on"

Here is another nice example (from http://rwiki.sciviews.org/doku.php?id=tips:programming:ifelse)

> varD <- log(0:9)  # A vector like this one that contains -Inf could be troublesome in subsequent manipulations
> varD
 [1]      -Inf 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379 1.7917595
 [8] 1.9459101 2.0794415 2.1972246
> varD <- ifelse(is.finite(varD), varD, NA)   # Remove the -Inf, replace with NA
> varD
 [1]        NA 0.0000000 0.6931472 1.0986123 1.3862944 1.6094379 1.7917595
 [8] 1.9459101 2.0794415 2.1972246

normalization, normalization, normalization!

2011-09-16T13:17:00.000-07:00

Basic Microarray Analysis:

Experiment (Raw Data)
Image Analysis (Pick Spots)
Pre-normalization analysis
1. Background subtraction
2. Spatial intensity
3. Chip to Chip variance
Normalization and Filtering
1. Scale normalization
2. Intensity dependent location normalization
Post-normalization analysis

Clustering
Hypothesis testing
Classification
More complicated statistical analysis

I found this list when I was searching the interquartile range (IQR) normalization. Actually Microarray folks might have been done the whole bunch of analysis decades ago. But the question for me is still there: how to compare two samples?

Let's say, I have two vectors: X=[x1, x2, x3, ..., xn] and Y=[y1, y2, y3, ... yn]. Before I could say X is correlated with Y simply because cor(X, Y) is large, I should make sure that each Xi and Yi are sampled from the same distribution. Right? That's the point. If they are not, I should do normalization ahead to assure that.

Maybe I should go thru the above list (esp. the normalization part) in some weekend.

Human Body Map, seq/chip

2011-09-16T09:39:00.000-07:00

Sequencing assays provide digital measures of sequence
abundance, i.e., read counts. In contrast, microarrays provide
analog measures of sequence abundance, i.e., intensities.

I read this in Sandrine Dudoit's presentation. Quite agree!

-------

Another piece of note for today might be old new, but still useful. -- Ensembl v62 has embeded all "Human Body Map 2.0" RNAseq data. Here is a nice blog for that.

http://www.ensembl.info/blog/2011/05/24/human-bodymap-2-0-data-from-illumina/

about high-throughput sequencing

2011-09-11T21:41:00.000-07:00

multiplexing: This term refers to the ability of sequencing multiple samples at the same time. For example, for small genome like yeast, C.elegant and drosophila, the number of reads generated in a sequencing unit (e.g. one of the 8 lanes in Illumina Genome Analyzer) may be several times of reads needed to provide a sufficient coverage of the genome.
barcode/index: To distinguish the different samples in same run/lane, a short sequence, call barcode/index, is added to the front of adaptor ligated together to different samples during sample preparation. "This is generally added to the 3' end of the upstream adapter so that this is the first 4-6 bases read during the sequencing run, then you sort the data by these first bases into your groups. "(http://www.umassmed.edu/DeepSeq_FAQ.aspx).

Here is the barcoding system for SOLiD:

For Illumina, see the slides in presentation:

https://www.illumina.com/documents/seminars/presentations/2010-06_sq_21_morris_multiplex_target_enrichment.pdf

Two related papers:

http://www.nature.com/nrg/journal/v10/n10/abs/nrg2641.html

http://www.nature.com/nmeth/journal/v5/n10/full/nmeth.1251.html

Several basic statistics terms

2011-09-07T20:40:00.000-07:00

Deviation: difference between observation and mean of the population, e.g Xi-E(X)

Residual: the difference between observed value and predicted/estimated value.

Error: the difference btw the observed value and the real value (which is usually unobservable)

MSE (Mean square error): the E[(X-X')^2], where X' is the real value of X.

Standard deviation: square root of the mean squared deviation. eg. sqrt(E[(Xi-u)^2]). It's square root of variance (see below)

Variance: a measurement of the dispersion for a set of numbers/observations (e.g. how far they are spread out from each other). If they have a mean, then variance is defined as E[(Xi-u)^2], where u is E(X). It's square of standard deviation.

Density curve of histogram plot in R

2011-09-01T20:23:00.000-07:00

Ref: http://casoilresource.lawr.ucdavis.edu/drupal/book/export/html/23

To add density curve on a histogram, like the green curve above, use code below:

#plot the distribution
hist(slope, breaks=1000, freq=F, main=main, xlab="Slope Value (percent)", xlim=c(0,150), ylim=c(0,.05) )
lines(density(slope, bw=1), col="green")

Add text aligned to legend in R plot

2011-09-01T12:45:00.000-07:00

What I meant is to add text on a R plot when there is already legend on it. Like the left plot in above figure, another piece of text was put exactly below the legend "Pearson'r ...RMSE = 1.9". Here is the code for that:

l=legend("topleft", paste("Pearson's r =", R, "(p-value:", pvalue , ")\nRMSE =", RMSE), bty="n", cex=2)

text(x=l$text$x, y=l$text$y-l$rect$h/2, "text you want to put", adj=c(0,1))

Ensemble learning, bagging, boosting

2011-08-29T12:35:00.001-07:00

Ensemble learning - Wikipedia, the free encyclopedia:

"an ensemble is a technique for combining many weak learners in an attempt to produce a strong learner"

Bagging, abbreviation of 'bootstrap aggregating', trains each model in the ensemble using a randomly-drawn subset of the training set. So, each model is independent. Random forest algorithm, for instance, combines random decision tree with bagging to get a high classification accuracy.

Boosting, is to build the ensembl by incrementally including new model instances to emphasize the training instances that were mis-classified by previous models. Common implementation of Boosting is Adaboost.

Other understanding from the randomForest documentation:

Recently there has been a lot of interest in “ensemble learning” — methods that generate many classiﬁers and aggregate their results. Two well-known methods are boosting (see, e.g., Shapire et al., 1998) and bagging (Breiman, 1996) of classiﬁcation trees. In boosting, successive trees give extra weight to points incorrectly predicted by earlier predictors. In the end, a weighted vote is taken for prediction. In bagging, successive trees do not depend on earlier trees— each is independently constructed using a bootstrap sample of the data set. In the end, a simple majority vote is taken for prediction.

Google Map API Exercise 1: display a static google map via URL

2011-08-25T20:20:00.000-07:00

这是一张以俺家为中心的static map，center的经纬度是通过javascript:void(prompt('',gApplication.getMap().getCenter()));获取的。其代码如下：

<img src=http://maps.googleapis.com/maps/api/staticmap?center=42.286959,-71.765538&zoom=17&size=512x512&maptype=satellite&markers=icon:http://bit.ly/pu5LFb|42.286959,-71.765538&sensor=true>

其中被shorten的地址(http://bit.ly/pu5LFb)的原地址为：http://chart.apis.google.com/chart?chst=d_map_spin&chld=0.84|30|FFFF88|10|_|%E4%BF%BA%E5%AE%B6

Estrogen receptor(ER) and hormone

2011-08-25T19:24:00.000-07:00

Sex hormone is one class of hormones, which includes estrogen, androgen, and progestogen. Hormone are important for inter-cell signaling.
Estradiol (E2) is one of estrogen, which can enter cell freely and bind to its receptor.
Estrogen receptor (ER) is a group of receptor that can be activated by estrogen. There are two types of ER, ERα and ERβ, both of which can form a dimer after activated by hormone. So, there are ERα (αα) or ERβ (ββ) homodimers or ERαβ (αβ) heterodimer. ER-alpha mostly occurs in endometrium, breast cancer cells, ovarian stroma cells, and the hypothalamus, where ER-beta are found in kidney, brain, bone, heart, lungs, intestinal mucosa, prostate, and endothelial cells.
Without sex hormone, ER are largely located in cytosol. When it's bound by hormone, a numer of events "starting with migration of the receptor from the cytosol into the nucleus, dimerization of the receptor, and subsequent binding of the receptor dimer to specific sequences of DNA (called hormone response elements)" are triggered. The receptor/DNA complex can then recruit other factor to initiate/activate the transcription of downstream genes.

display time series data in R

2011-08-19T08:04:00.000-07:00

Thanks to the Revolutions blog, several things learned here:
1. R code for heat calendar
2. generate SVG from R
3. pretty-R tool

OK. Let's explain it by plotting the fund WASCX (IVY ASSET STRATEGY FUND CLASS C) from 2009-03-14:

# source code of calendarHeat (by Paul Bleicher)
source("http://blog.revolution-computing.com/downloads/calendarHeat.R")
 
# fetch stock time series data from Yahoo finance
stock <- "WASCX"
start.date <- "2009-03-14"
end.date <- Sys.Date()
quote <- paste("http://ichart.finance.yahoo.com/table.csv?s=",
                stock,
                "&a=", substr(start.date,6,7),
                "&b=", substr(start.date, 9, 10),
                "&c=", substr(start.date, 1,4), 
                "&d=", substr(end.date,6,7),
                "&e=", substr(end.date, 9, 10),
                "&f=", substr(end.date, 1,4),
                "&g=d&ignore=.csv", sep="")             
stock.data <- read.csv(quote, as.is=TRUE)
 
# display in SVG format
svg("WASCX.svg",width=14,height=7)
calendarHeat(stock.data$Date, stock.data$Adj.Close, 
  varname="WASCX Adjusted Close")
dev.off()

Created by Pretty R at inside-R.org

One R Tip A Day: How to draw a plot with two Y axises and one X axis

2011-08-17T11:22:00.000-07:00

One R Tip A Day: How to draw a plot with two Y axises and one X axis

plot(1:10)
par("usr")
# [1] 0.64 10.36 0.64 10.36
# Now resetting y axis' usr coordinates:
par(usr=c(par("usr")[1:2], 101, 105))
points(1:5, 105:101, col="red")
axis(4)

par(mar=c(4, 5, 4, 5) + 0.1)
plot(y1, type='b', col='gray', ylim=range(c(y1, y2)), ylab="y")
lines(y2, type='b', pch=19)
# draw the other Y axis
par(usr=c(par("usr")[1:2], range(Y)))
lines(Y, type='b', pch=20, col='blue')
axis(4, col.axis='blue')
# manually label the other Y axis
mtext(side=4, line=1, 'Y', col='blue')
legend('topleft', c("y1", "y2", "Y"), pch=c(19, 1, 20), col=c('black', 'gray', 'blue'), lty=1)

A very smart way, rather than transforming one Y axis range to the other Y axis range by math formula.

#Y = (Y-min(Y))*(max(y1, y2) - min(y1, y2))/(max(Y)-min(Y)) + min(y1, y2)

How to run Meme on cluster

2011-08-17T09:23:00.000-07:00

The latest version of Meme, v4.6.1, was built with mpich, while the earlier versions were built with openmpi. You need to set the path correctly as in the examples below.

Your input database should consist of a file containing sequences in fasta format. In the example below, the file is 'mini-drosoph.s'.
Maxsize parameter: The maximum dataset size in characters. Determine the number of characters in your dataset by typing 'wc -c filename'. e.g.

[user@biowulf mydir]$ wc -c mini-drosoph.s 

506016 mini-drosoph.s

For this dataset, the maxsize parameter has to be set to greater than 506,016, so we will use 600000.Set up a batch script along the lines of the ones below:

Batch script for Meme 4.6.1

----  this file is called meme.batch ---------

#!/bin/bash

#This batch script can be used with Meme 4.6.1

#PBS -N Meme

#PBS -m be

#PBS -j oe



export PATH=/usr/local/mpich-1.2.7p1-gcc4_64/bin:$PATH

cd /data/user/meme/

time mpirun -machinefile $PBS_NODEFILE -np $np /usr/local/meme_4.6.1/bin/meme_p \

     /data/user/meme/test.fa -oc /data/user/meme/meme_out \

     -maxsize 10000000 -p $np

// The above tips are from http://biowulf.nih.gov/apps/meme.html

btw, the above script use PBS (Portable Batch System) submission system. Alternatively, there are other options, like LoadLeveler (Champion), LSF (Lonestar). There are syntax comparison here. Also, more detail for SGE and PBS here.

Difference between FIMO and MAST

2011-08-15T09:41:00.000-07:00

FIMO and MAST are two applications in the MEME suites.

FIMO is looking for the best individual matches to motifs. MAST is looking for sequences the have the best overall match to a collection of motifs.

FIMO's task is simple: given a set of motifs and a database of sequences, compute the match score to each motif at each position in each sequence, and report all the motif matches that pass the p-value/q-value threshold.

MAST's algorithm is more complex. For each sequence it carries out an initial scoring that is quite similar to FIMO's. MAST then picks the best match for each motif in the sequence. The p-values for these top matches are multiplied together to create an overall score for the full sequence. MAST reports the sequences that have the most significant overall scores. Typically MAST would be used to look for regulatory regions in DNA, or structures in proteins, where several motifs might occur near each other.

-- Thanks to Charles for the explanation

Two ways of javascript function declaration

2011-08-09T19:40:00.000-07:00

1. The first way uses a function name, like:

2. The second way assigns a function without name to a variable, for example:

Here are more details:
http://stackoverflow.com/questions/336859/javascript-var-functionname-function-vs-function-functionname

An efficient way to do dataset intersection

2011-07-27T09:33:00.000-07:00

The main message is to use "match" to get index of needed rows and then get the rows by the index, instead of using the row names to select, which is much slower. Here is example:

In example above, we know that the same values of column 2nd have same values of columns from 4th to the end. So, instead of doing unique on whole matrix, getting the unique of column 2nd and then getting the index of unique ones by match. Match(a,b) only return the index of first occurrence of a in b. For example

This tips also help in intersecting two big dataframes. For example,

Linux Tutorial - Apache Web Login Authentication:

2011-07-24T20:30:00.000-07:00

Linux Tutorial - Apache Web Login Authentication::

Password protection by a single login:

Password files:

Create the directory you want to password protect (example: membersonly)

Create a file /home/domain/public_html/membersonly/.htaccess in that director that looks something like this:

AuthName "Add your login message here."

AuthType Basic

AuthUserFile /home/domain/public_html/membersonly/.htpasswd

AuthGroupFile /dev/null

require user name-of-user

In this case the "name-of-user" is the login name you wish to use for accessing the web site.

Create (or clobber if it already exists) the password file /home/domain/public_html/membersonly/.htpasswd using the program htpasswd:

htpasswd -c .htpasswd name-of-user

How to choose a statistical test?

2011-07-24T20:26:00.000-07:00

Nice table from http://www.graphpad.com/www/book/choose.htm

Table 37.1. Selecting a statistical test

	Type of Data
Goal	Measurement (from Gaussian Population)	Rank, Score, or Measurement (from Non- Gaussian Population)	Binomial (Two Possible Outcomes)	Survival Time
Describe one group	Mean, SD	Median, interquartile range	Proportion	Kaplan Meier survival curve
Compare one group to a hypothetical value	One-sample ttest	Wilcoxon test	Chi-square or Binomial test **
Compare two unpaired groups	Unpaired t test	Mann-Whitney test	Fisher's test (chi-square for large samples)	Log-rank test or Mantel-Haenszel*
Compare two paired groups	Paired t test	Wilcoxon test	McNemar's test	Conditional proportional hazards regression*
Compare three or more unmatched groups	One-way ANOVA	Kruskal-Wallis test	Chi-square test	Cox proportional hazard regression**
Compare three or more matched groups	Repeated-measures ANOVA	Friedman test	Cochrane Q**	Conditional proportional hazards regression**
Quantify association between two variables	Pearson correlation	Spearman correlation	Contingency coefficients**
Predict value from another measured variable	Simple linear regression or Nonlinear regression	Nonparametric regression**	Simple logistic regression*	Cox proportional hazard regression*
Predict value from several measured or binomial variables	Multiple linear regression* or Multiple nonlinear regression**		Multiple logistic regression*	Cox proportional hazard regression*

Remove zero bytes files in folder

2010-11-20T19:38:00.000-08:00

Remove zero bytes files in folder - The UNIX and Linux Forums: "

find . -type f -size 0 | xargs rm

or

find . -type f -size 0 -exec rm {} \;

awk.gsub, awk.system

2010-11-20T15:47:00.000-08:00

Tips in Computer: "awk '{system('grep ' $4 ' file > result.txt'); getline result < 'result.txt'; close('result.txt'); print result;}' infile

e.g File
/scratch/dongx/jobid_2622600/18_2_62LPFAAXX_092810_MyersLab_4242.fastq.tophat.mm_1.unique.hg19.out/accepted_hits.bam
/scratch/dongx/jobid_2622600/18_2_62LPFAAXX_092810_MyersLab_4242.fastq.tophat.mm_1.unique.hg19.out/junctions.bed
/scratch/dongx/jobid_2622604/20_2_62MA2AAXX_101210_MyersLab_4430_4.fastq.tophat.mm_1.unique.hg19.out/accepted_hits.bam
/scratch/dongx/jobid_2622604/20_2_62MA2AAXX_101210_MyersLab_4430_4.fastq.tophat.mm_1.unique.hg19.out/junctions.bed
/scratch/dongx/jobid_2622606/20_5_62MA2AAXX_101210_PhiX.fastq.tophat.mm_1.unique.hg19.out/accepted_hits.bam

You want to copy the .bam and .bed file to another folder...

$ cat > tophat.log
$ grep "s\.b" tophat.log | awk '{a=$1; gsub(".*[0-9]/","", a); sub("/",".", a); system("cp " $1 " " a);}'

- Sent using Google Toolbar"

asparapee (芦笋尿)

2010-11-08T18:28:00.000-08:00

asparapee (芦笋尿)_我和你_百度空间: "今日看到一个很有趣的词——asparapee。这个俚语是asparagus+pee构成的一个复合词，直译就是“芦笋尿”。国外不少人喜欢吃芦笋，好吃而且对健康有好处。但是吃了芦笋之后，尿会有恶臭，臭到自己都受不了:) 科学上还没有一个确切的解释，大概是因为asparagus中含有硫化物或者某些氨基酸(amino acids)，可能和农民用硫来肥沃土壤有关。专家相信带有某种特殊基因的人能制造某种消化酶，这种酶可以把asparagus分解成某种有异味的氨基酸；而且带有某种基因的人，才能闻到这种气味。

让我们看看Urban Dictionary的解释：

http://www.urbandictionary.com/define.php?term=Asparapee

The noxious smell in urine experienced by most people immediately after consuming asparagus. Studies have shown that while many can’t smell the pungent stench, most do extrude it. This phenomenon is most likely caused by farmers adding sulfur to fertilize their soil. The odor is sometimes accompanied by a neon yellow or green hue.

Person A: Damn, your pee is freakin' reekin'! Did you have asparagus for dinner last night?

Person B: Yes, I have asparapee.

- Sent using Google Toolbar"

Parallel file downloading with wget | commandlinefu.com

2010-11-05T13:10:00.000-07:00

Parallel file downloading with wget | commandlinefu.com:

"Parallel file downloading with wget

$ wget -nv http://en.wikipedia.org/wiki/Linux -O- | egrep -o "http://[^[:space:]]*.jpg" | xargs -P 10 -r -n 1 wget -nv

xargs can be used in this manner to download multiple files at a time, and xargs will in this case run 10 processes at a time and initiate a new one when the number running falls below 10.

- Sent using Google Toolbar"

Make grep stop after the first match - The macosxhints Forums

2010-11-05T12:23:00.000-07:00

Make grep stop after the first match - The macosxhints Forums:

Make grep stop after the first match

Question: I'm trying to use grep to find all the href links in my asp files (except in-page links that start with #).
I want it to show me the href="blahblah" part only - nothing before href or after href's closing quote.
Here's what I have so far. This:

grep -i -r -o 'href="[^#].*"' rentals.asp

gives me this (green=good, blue=bad):

href="<%= webaddr %>/rentalStyle.css" type="text/css"
href="<%= webaddr %>/rentals_s.asp"
href="<%= webaddr %>/rentals_o.asp">rentals
href="<%= webaddr %>/search.asp">Search

Solution:
"I think the question has to do with the 'greediness' of grep when multiple matches occur in the same line (in this case, the double-quote character). To match everything up to the first closing quotes, try:
Code:

grep -i -r -o 'href='[^#][^']*'' rentals.asp

----------------
Similar question to get URL links in html source code like:

< href="ftp://encodeftp.cse.ucsc.edu:21/users/akundaje/rawdata/peaks/jul2010/idr0_02/narrowPeak_blacklistFiltered/wgEncodeBroadHistoneGm12878CtcfStdAlnRep0.bam_VS_wgEncodeBroadHistoneGm12878ControlStdAlnRep0.bam.regionPeak.gz">wgEncodeBroadHistoneGm12878CtcfStdAlnRep0.bam_VS_wgEncodeBroadHistoneGm12878ControlStdAlnRep0.bam.regionPeak.gz< / a >

hint code from commandlinefu :
egrep -o "ftp://[^[:space:]]*.gz"
will get
ftp://encodeftp.cse.ucsc.edu:21/users/akundaje/rawdata/peaks/jul2010/idr0_02/narrowPeak_blacklistFiltered/wgEncodeSydhTfbsK562bTr4UcdAlnRep0.bam_VS_wgEncodeSydhTfbsK562bInputUcdAlnRep1.bam.regionPeak.gz">wgEncodeSydhTfbsK562bTr4UcdAlnRep0.bam_VS_wgEncodeSydhTfbsK562bInputUcdAlnRep1.bam.regionPeak.gz

So, it should be corrected into :
egrep -o "ftp://[^[:space:]^\"]*.gz"

Smart, ha!
- Sent using Google Toolbar"

One R Tip A Day: Sorting/ordering a data.frame according specific columns

2010-11-01T14:36:00.000-07:00

One R Tip A Day: Sorting/ordering a data.frame according specific columns:

x = rnorm(20)
y = sample(rep(1:2, each = 10))
z = sample(rep(1:4, 5))

data.df <- data.frame(values = x, labels.1 = y, labels.2 = z)
print(data.df)

# data ordered according to 'labels.1' column
# and then 'labels.2' column
nams <- c('labels.1', 'labels.2')
data.df.sorted = data.df[do.call(order, data.df[nams]), ]
print(data.df.sorted)

===================== OR ====================

dd <- data.frame(b = factor(c("Hi", "Med", "Hi", "Low"),
      levels = c("Low", "Med", "Hi"), ordered = TRUE),
      x = c("A", "D", "A", "C"), y = c(8, 3, 9, 9),
      z = c(1, 1, 1, 2))
dd
    b x y z
1  Hi A 8 1
2 Med D 3 1
3  Hi A 9 1
4 Low C 9 2

R> dd[with(dd, order(-z, b)), ]
    b x y z
4 Low C 9 2
2 Med D 3 1
1  Hi A 8 1
3  Hi A 9 1

- Sent using Google Toolbar"

Interesting packages of new BioC2.7

2010-11-01T10:10:00.000-07:00

NuPoP

Nucleosome positioning prediction

Mulcom

Differential expression and false discovery rate calculation through multiple comparison

ontoCAT

Ontology parsing

BHC

Bayesian Hierarchical Clustering

iSeq

Bayesian Hierarchical Modeling of ChIP-seq Data Through Hidden Ising Models

URL: http://bioconductor.org/news/bioc_2_7_release/

btw, other interesting ones in BioC 2.6:

New sequence analysis tools address infrastructure (GenomicRanges, Rsamtools, girafe); ChIP-seq (BayesPeak, CSAR, PICS); digital gene expression and RNA-seq (DESeq, goseq, segmentSeq); and motif discovery (MotIV, rGADEM).

Running R on cluster in a job_array way

2010-10-28T16:09:00.000-07:00

Simple-Job-Array-Howto - GridWiki: "Example: R Scripts with Grid Engine Job Arrays

All of the above applies to well-behaved, interactive programs. However, sometimes you need to use R to analyze your data. In order to do this, you have to hardcode file names into the R script, because these scripts are not interactive. This is a royal pain. However, there is a solution that makes use of HERE documents in bash. HERE documents also exist in perl, and an online tutorial for them in bash is at http://www.tldp.org/LDP/abs/html/here-docs.html. The short of it is that a HERE document can represent a skeleton document at the end of a shell script. Let’s concoct an example. You have 100 data files, labeled data.1 to data.10. Each file contains a single column of numbers, and you want to do some calculation for each of them, using R. Let’s use a HERE document:

#!/bin/sh
#$ -t 1-10
WORKDIR=/Users/jl566/testing
INFILE=$WORKDIR/data.$SGE_TASK_ID
OUTFILE=$WORKDIR/data.$SGE_TASK_ID.out
# See comment below about paths to R
PATHTOR=/common/bin
if [ -e $OUTFILE ]
then
rm -f $OUTFILE
fi
# Below, the phrase "EOF" marks the beginning and end of the HERE document.
# Basically, what’s going on is that we’re running R, and suppressing all of
# it’s output to STDOUT, and then redirecting whatever’s between the EOF words
# as an R script, and using variable substitution to act on the desired files.
$PATHTOR/R --quiet --no-save > /dev/null << EOF
x<-read.table("$INFILE")
write(mean(x\$V1),"$OUTFILE")
EOF

So now you can use the cluster to analyze your data – just write the R script within the HERE document, and go from there. As I’ve only just figured this out, some caveats are necessary. If anyone experiments and figures out something neat, let me know. Be aware of the following:

1. In my limited experience, indenting is important for HERE documents. In particular, it seems that the beginning and end (i.e. both lines containing the term EOF in the above example), must be aligned with the left-hand edge of the buffer (i.e. not indented at all). So, if you use a HERE document in a conditional or control statement, be mindful of this.
2. In the mean command, I escaped the dollar sign with a backslash. In my limited experiments, both mean(x\$V1) and mean(x$V1) seem to work. However, escaping the dollar sign for the read.table command prevents the variable substitution from occurring in the shell, causing R to fail, because the input file named $INFILE cannot be found. In other words, escaping in that context causes the HERE doc to pass $INFILE as a string literal to R, rather than the value stored in the shell variable.
3. This is more useful than just array jobs on an SGE system. If you know bash well enough, you can write a shell script that takes a load of arguments, and processes them with a HERE document. This solves a major limitation with R scripts themselves. You can do the same in perl, too, on your workstation, but you must use a shell language on the cluster.

- Sent using Google Toolbar"

Linux tip: Bash parameters and parameter expansions

2010-10-11T12:35:00.000-07:00

Linux tip: Bash parameters and parameter expansions: "Table 2. Shell expansion pattern matchingExpansion Purpose
${PARAMETER#WORD} The shell expands WORD as in filename expansion and removes the shortest matching pattern, if any, from the beginning of the expanded value of PARAMETER. Using '@' or '$' results in the pattern removal for each parameter in the list.
${PARAMETER##WORD} Results in removal of the longest matching pattern from the beginning rather than the shortest.
${PARAMETER%WORD} The shell expands WORD as in filename expansion and removes the shortest matching pattern, if any, from the end of the expanded value of PARAMETER. Using '@' or '$' results in the pattern removal for each parameter in the list.
${PARAMETER%%WORD} Results in removal of the longest matching pattern from the end rather than the shortest.
${PARAMETER/PATTERN/STRING} The shell expands PATTERN as in filename expansion and replaces the longest matching pattern, if any, in the expanded value of PARAMETER. To match patterns at the beginning of the expanded value of PARAMETER, prefix PATTERN with # or prefix it with % if the match should be done at the end. If STRING is empty, the trailing / may be omitted and the matches are deleted. Using '@' or '$' results in the pattern substitution for each parameter in the list.
${PARAMETER//PATTERN/STRING} Performs the substitution for all matches instead of just the first.

Listing 11 shows some basic usage of the pattern matching expansions.

Listing 11. Pattern matching examples

[ian@pinguino ~]$ x='a1 b1 c2 d2'
[ian@pinguino ~]$ echo ${x#*1}
b1 c2 d2
[ian@pinguino ~]$ echo ${x##*1}
c2 d2
[ian@pinguino ~]$ echo ${x%1*}
a1 b
[ian@pinguino ~]$ echo ${x%%1*}
a
[ian@pinguino ~]$ echo ${x/1/3}
a3 b1 c2 d2
[ian@pinguino ~]$ echo ${x//1/3}
a3 b3 c2 d2
[ian@pinguino ~]$ echo ${x//?1/z3}
z3 z3 c2 d2

- Sent using Google Toolbar"

parameters and options

2010-10-11T11:16:00.000-07:00

In BASH shell script:
$0 : the program itself
$1 : the first parameter
$2 : the second parameter
...
A string enclosed in single or double quotes will be passed as a single parameter, and the quotes will be stripped.
$# : the number of total parameters, not including the program itself
$@ : the array of parameters. Double-quoted "$@" is equal to "$1", "$2", "$3".... This is the most common usage. But there are other variations, like $*, "$*", $@, or "$@". They can be explained distinctly. For more detail, see examples here.

For options, using
getopt optstring optname

where optstring is a double-quote string of option letters, like ":p:q". A colon (:) after an option letter indicates that the option
requires a value; The leading colon tells getopts to be silent and suppress the
normal error messages.

The second parameter, optname, is the name of a variable which will receive the name of the option found. If an option is expected to have a value, the value, if present, will be placed in the variable OPTARG. In silent mode, either of the following two error conditions may occur.

If an unrecognized option is found, then optname will contain a ? and OPTARG will contain the unknown option.
If an options that requires a value is found but the value is not, then optname will contain a : and OPTARG will contain the name of the option whose argument is missing.

For example, getopt ":p:q" optnames
will require an option (-p) with value, and -q without having value.

Example: 

#!/bin/bash
echo "OPTIND starts at $OPTIND"
while getopts ":pq:" optname
do
  case "$optname" in
    "p")
      echo "Option $optname is specified, with value"
      ;;
    "q")
      echo "Option $optname has value $OPTARG"
      ;;
    "?")
      echo "Unknown option $OPTARG"
      ;;
    ":")
      echo "No argument value for option $OPTARG"
      ;;
    *)
    # Should not occur
      echo "Unknown error while processing options"
      ;;
  esac
  echo "OPTIND is now $OPTIND"
done

file/folder existence

2010-10-11T10:52:00.000-07:00

General syntax in BASH shell is:

[ parameter FILE ]
OR
test parameter FILE

Where parameter can be any one of the following:

-e: Returns true value if file exists
-f: Return true value if file exists and regular file
-r: Return true value if file exists and is readable
-w: Return true value if file exists and is writable
-x: Return true value if file exists and is executable
-d: Return true value if exists and is a directory

For example,

Find out if file /etc/passwd file exists or not:

$ [ -f /etc/passwd ] && echo "File exists" || echo "File  does not exists"

# Make the directory for the job ID you are running if it does not exist
[ -d $HOME/scratch/jobid_$JOB_ID ] || mkdir -p $HOME/scratch/jobid_$JOB_ID

You can use conditional expressions in a shell script:

#!/bin/bash
FILE=$1

if [ -f $FILE ];
then
  echo "File $FILE exists"
else
  echo "File $FILE does not exists"
fi

See if a directory exists or not with NOT operator:
[ ! -d directory ]
OR
! test directory

Reorder one data.frame using two columns from another data.frame in R - Stack Overflow

2010-09-28T20:26:00.000-07:00

Reorder one data.frame using two columns from another data.frame in R - Stack Overflow: "b[match(a$id,b$id),]

Unfortunately, I'm not sure how to accomplish the same thing when I need to order by two columns.

a:

id lob
1 1+ X
2 3 X
3 2 X
4 1 X
5 1 Y
6 1+ Y
7 1+ X
8 3 X
9 3 X

b:

id lob val
1 1+ X 1
2 1+ Y 9
3 1+ X 1
4 3 X 5
5 3 X 5
6 3 X 5
7 2 X 4
8 1 X 3
9 1 Y 2

I want to get this:

id lob val
1 1+ X 1
2 3 X 5
3 2 X 4
4 1 X 3
5 1 Y 2
6 1+ Y 9
7 1+ X 1
8 3 X 5
9 3 X 5

- Sent using Google Toolbar"

How to replace space using sed in shell script

2010-09-28T13:11:00.000-07:00

How to replace space using sed in shell script: "Use the + modifer to match one or more occurences of the previous regular expression. It needs to be escaped:
Code:

$ echo '12:Mit 23627 83763 tty message' | sed 's/ \+/:/g'
12:Mit:23627:83763:tty:message

$echo '12:Mit 23627 83763 tty message' | sed 's/ /:/g'| cat -A
12:Mit:23627:::83763:tty::::message$

$ echo '12:Mit 23627 83763 tty message' | sed 's/ +/:/g'
12:Mit 23627 83763 tty message

- Sent using Google Toolbar"

PASS parameter to AWK - The UNIX and Linux Forums

2010-09-25T12:42:00.000-07:00

PASS parameter to AWK - The UNIX and Linux Forums: "You can also assign variables on the command line, e.g....

awk -v a=1 -v b=2 -v c=3 -f test.awk

The -v option means that the variables are assigned before the BEGIN section.

- Sent using Google Toolbar"

if you..., I would like to pay cold hard cash?

2010-09-25T10:15:00.000-07:00

Answers.com - What does the phrase 'cold hard cash' mean: "Back when coinage used to be made from precious metals like gold and silver, the value of the metal itself was used to guarantee the worth of the coin; some people didn't trust paper money, as inflation or war could quickly make it worthless.

These days, the idea has evolved into the idea of using money, including paper money, rather than a check or a credit card -- both of which have a higher risk of being fraudulent, or which can be more easily canceled after the sale. Money, once handed over, can't be taken back so easily.

- Sent using Google Toolbar"

说白了，就是中国话的“真金白银”

firefox add on blogspot - Google Search

2009-06-10T06:43:00.000-07:00

firefox add on blogspot - Google Search: "#
Url Shortening Firefox Addon with Google Analytics Support ~ Web Upd8
25 Apr 2009 ... Snip-n-Tag is a Firefox addon to easily shorten URLs, plus add Google ... on Post Author for a Wordpress / Blogger (blogspot.com) Blog ...
webupd8.blogspot.com/2009/04/url-shortening-firefox-addon-with.html - Cached - Similar pages -
#
JDownloader 0.5.859 - New Interface, Firefox Addon and Many other ...
JDownloader 0.5.859 - New Interface, Firefox Addon and Many other Changes. ... The new version also comes with a Firefox plugin which can automatically handle .... Depening on Post Author for a Wordpress / Blogger (blogspot.com) Blog ...
webupd8.blogspot.com/2009/05/jdownloader-05859-new-interface-firefox.html -"

Illustration of linux command

2008-09-23T04:46:00.000-07:00

AWK learning note

2008-09-10T02:17:00.000-07:00

1. first update the internal variables when reading one line:

當 AWK 從資料檔中讀取一筆資料列時, AWK 會使用內建變數$0 予以記錄.
AWK 會立刻重新分析 $0 的欄位情況, 並將 $0 上各欄位的資料用 $1, $2, ..予以記錄.


例如 : AWK 從資料檔 emp.dat 中讀入第一筆資料列
"A125 Jenny 100 210" 之後, 程式中:
$0 之值將是 "A125 Jenny 100 210"
$1 之值為 "A125" $2 之值為 "Jenny"
$3 之值為 100 $4 之值為 210
NF 之值為 4 $NF 之值為 210
NR 之值為 1 FILENAME 之值為 ``emp.dat''
where NF: Number of Fields in current $0
NR: Number of Records of currently having been read.
FILENAMEAWK: filename of current proceeding

2. 'PATTERN{ACTION}' or -f script.awk
the following two ways are same:

$awk -f pay1.awk emp.dat
$awk ' { print $2, $3 * $4 } ' emp.dat

if you save the script into a file named pay1.awk.
讀者可使用``-f''參數,讓AWK主程式使用其它僅含 AWK函數 的
檔案中的函數
其語法如下:
awk -f AWK主程式檔名 -f AWK函數檔名 資料檔檔名

3. BEGIN/END and array in AWK
for example, we have a data file like:

Mary O.S. Arch. Discrete
Steve D.S. Algorithm Arch.
Wang Discrete Graphics O.S.
Lisa Graphics A.I. Lily Discrete Algorithm

---------------------------------------

{for( i=2; i<>
END{
for(coursein Number)
printf("\%-10s %d\n", course, Number[course] )
}

---------------------------------------

comment:
a. NF=4 in this case, line number
b. END is a AWK之保留字, 為{ Pattern}之一種, like BEGIN. The only difference is END only run after all lines are proceeded, while BEGIN works initially before the script, and only one time (both BEGIN and END).
c. $i represents the ith elements in the line array, which is different from Perl program (in which, the $i is a variable name, in AWK, variable name cannot begin with $.)

4. Shell command and awk command
for example:

---------------------------------------

BEGIN {
while ( "who" | getline ) n++
print n
}

---------------------------------------
where the who is a system command used in shell, and the getline is an awk command for input;

5. Filename in the script should be quoted by "",
for example,

---------------------------------------

BEGIN {
print `` ID Number Arrival Time'' > ``today_rpt1''
print ``==========================='' > ``today_rpt1''
}

{ printf(" %s %s\n", $1,$2 ) > "today_rpt1" }

---------------------------------------

$awk -f reformat1.awk arr.dat

Note:
a. if today_rpt1 is not quoted by "", then it will be taken as a variable (which default value is 0, or Null String in AWK.)
b. the redirection mark is '>', not '>>‘, even you want to append to the end of the file. The only difference between them is, for '>>', it will append to the end of the file if it's open first time and the file exists. For '>', AWK will create a new file when it occurs first time, then append to the end (like '>>'). This is little bit different from Unix.

6. Input and output command in Awk
AWK input command: getline
AWK output command: print, printf

7. three ways to run awk
a. $awk '{print}' file1.txt file2.txt
b. $awk -f myscript.awk file1.txt file2.txt
save {print} into a file(myscript.awk) first
c. $myshell file1.txt file2.txt
save awk '{print}' $* into a shell file(named myshell. Here $* means all parameters after the shell command. You also can use $1 represents the first parameter, and $2 the second one.

8. FS(Field Separator) and RS(Record Separator)
By default, the FS is any empty character (space, \t, ), RS is newline '\n'. But they can be changed, like

--------------------------------------- make_report.awk -------------------------

BEGIN {
FS = "\n"
RS = ""
split( "一. 二. 三. 四. 五. 六. 七. 八. 九.", C_Number, " " )
}
{
printf("\n%s 報告人 : %s \n",C_Number[NR],$1)
for( i=2; i<= NF; i++)
printf(" %d. %s\n", i-1, $i)
}

--------------------------------------- week.rpt ------------------------------

張長弓

GNUPLOT 入門



吳國強

Latex 簡介

VAST-2 使用手冊

mathematica 入門









李小華

AWK Tutorial Guide Regular Expression
--------------------------------------- Output ------------------------
[xianjund@douglasgran data]$ awk -f make_report week.rpt

一. 報告人 : 張長弓
1. GNUPLOT    入門

二. 報告人 : 吳國強
1. Latex 簡介
2. VAST-2 使用手冊
3. mathematica 入門

三. 報告人 : 李小華
1. AWK Tutorial Guide Regular Expression
---------------------------------------
9. ARGC and ARGV[]

like C, but
a. ARGC does not include the -v, -f and their options. for example, in 
$awk -vx=36 -f program1 data1 data2
or
$awk '{ print $1 ,$2 }' data1 data2

ARGC=3
ARGV[0]= "awk"
ARGV[1]="data1"
ARGV[2]="data2"

png to ico

2008-09-09T16:10:00.000-07:00

How to convert PNG to ICO icon file?

1. download the png2ico source code

http://www.winterdrache.de/freeware/png2ico/

2. install

tar, cd, then make

3. prepare resized png, for example

convert image.png -resize 32x32 image.icon.png

4. png2ico

png2ico favicon.ico image.icon.png

You also can make an icon with multiple solution images, use

png2ico favicon.ico image.16x16.icon.png image.32x32.icon.png

reverse lines of file

2008-09-09T16:06:00.000-07:00

use rev filename. For example,

[xianjund@douglasgran hypotest]$ head test11.data

ENSG00000007372 ENSG00000109911 ENSG00000121690

ENSG00000007372 ENSG00000149100 ENSG00000170959

ENSG00000043355 ENSG00000175198 ENSG00000102452

ENSG00000172845 ENSG00000115840 ENSG00000091428

[xianjund@douglasgran hypotest]$ rev test11.data

09612100000GSNE 11990100000GSNE 27370000000GSNE

95907100000GSNE 00194100000GSNE 27370000000GSNE

25420100000GSNE 89157100000GSNE 55334000000GSNE

82419000000GSNE 04851100000GSNE 54827100000GSNE

It's cool, hmm?

to show multi line around grep result

2008-09-09T15:37:00.000-07:00

You can grep multiple lines before or after matching the keywords. Here is a simple tips, that what I discover grep capable of. A is after, B is before.

For example, test.data is like this:

ENSG00000007372 ENSG00000109911 ENSG00000121690

ENSG00000007372 ENSG00000149100 ENSG00000170959

ENSG00000043355 ENSG00000175198 ENSG00000102452

ENSG00000172845 ENSG00000115840 ENSG00000091428

ENSG00000172845 ENSG00000138430 ENSG00000128708

ENSG00000103449 ENSG00000103494 ENSG00000121274

ENSG00000104313 ENSG00000182674 ENSG00000140396

ENSG00000117707 ENSG00000136643 ENSG00000143499

ENSG00000121297 ENSG00000105176 ENSG00000178904

>grep ENSG00000138430 test.data

ENSG00000172845 ENSG00000115840 ENSG00000091428

while,

>grep ENSG00000138430 test.data -B1 -A3

ENSG00000172845 ENSG00000115840 ENSG00000091428

ENSG00000172845 ENSG00000138430 ENSG00000128708

ENSG00000103449 ENSG00000103494 ENSG00000121274

ENSG00000104313 ENSG00000182674 ENSG00000140396

ENSG00000117707 ENSG00000136643 ENSG00000143499

grep based on multiple words

2008-09-09T15:18:00.000-07:00

For example, grep "you" and "me" in the file,

grep "you\|me" filename

You need to put ESCAPE STRING ( \ ) for OR ( | ), else it will treat it as a simbol you want to search instead of regular expression symbol. or

grep -E "you|me" filename

egrep "you|me" filename

But, to get line with both "you" and "me", you can use

grep "you" filename | grep "me"

or,

egrep "you.*me" filename

but this will include those lines like "your lovely meebo", which is not what we want sometime. So, to get the exact words matched, use

egrep "\.*\" filename

For more info about egrep, use "man egrep"

The caret ^ and the dollar sign $ are metacharacters that respectively match the empty string at the beginning and end of a line. The symbols \<> respectively match the empty string at the beginning and end of a word. The symbol \b matches the empty string at the edge of a word, and \B matches the empty string provided it not at the edge of a word.

replace word in a file

2008-09-09T15:11:00.000-07:00

To replace a word in a file, use

perl -pi -e 's/abc/def/;' xyz

sed -e 's/abc/def/;' xyz > xyz_new

To make the replacement in place, use sed -i filename. then the file will be replaced after the command.

sed -e 's/abc/def/;' -i xyz

HOW TO REDIRECT OUTPUT TO A FILE AS WELL AS DISPLAY IT OUT

2008-09-09T14:45:00.000-07:00

To redirect standard output to a file is easy, you just need to use the redirection symbol, for example:

echo "hello world" > test.txt

But what if I want to display it out as well as store into a file?
Answer: tee

echo "hello world" | tee test.txt

Okay it seems very easy, how about append?

To append the standard output to a file, you do this:

echo"hello world" >> test.txt

Append to file and display it out as well?

echo"hello world" | tee -a test.txt

Okay, how about dealing with standard output(stdout) and standard error(stderr)?
There are two different output stream, one is stdout and another one is stderr. Normal print usually goes to stdout and error related message will goes to stderr. Lets make a simple python script to print 1 line to stdout and 1 line to stderr.

#!/usr/bin/env python  import sys  sys.stdout.write("I am stdout\n") sys.stderr.write("I am stderr\n")

Ok, lets save the python script as sout.py and try to redirect the output to a file.

$ ./sout.py > test.txt I am stderr

Standard output is redirect to test.txt but stderr is print out.

What if I want stderr to be redirect and display the stdout?

 ./sout.py 2> test.txt

I want both stored into the file.

 ./sout.py 2&> test.txt

At last, I want both display and redirect to a file:

 ./sout.py 2>&1 | tee test.txt

Interesting isn’t it?

Openning test

2008-09-09T14:34:00.000-07:00

Hope I can insist on sharing my learning points with all, in HTML, Perl, Unix/Linux and Bioinformatics :)