evolving all we are

2013年7月12日星期五

choose colors for your plots

1. The HTML 4.01 specification^[9] defines sixteen named colors, as follows (names are defined in this context to be case-insensitive):
http://en.wikipedia.org/wiki/Web_colors

CSS 1–2.0 / HTML 3.2–4 / VGA color names
Name	Hex (RGB)	Red (RGB)	Green (RGB)	Blue (RGB)	Hue (HSL/HSV)	Satur (HSL)	Light (HSL)	Satur (HSV)	Value (HSV)	CGA number (name); alias
White	#FFFFFF	100%	100%	100%	0°	0%	100%	0%	100%	15 (white)
Silver	#C0C0C0	75%	75%	75%	0°	0%	75%	0%	75%	7 (light gray)
Gray	#808080	50%	50%	50%	0°	0%	50%	0%	50%	8 (dark gray)
Black	#000000	0%	0%	0%	0°	0%	0%	0%	0%	0 (black)
Red	#FF0000	100%	0%	0%	0°	100%	50%	100%	100%	12 (high red)
Maroon	#800000	50%	0%	0%	0°	100%	25%	100%	50%	4 (low red)
Yellow	#FFFF00	100%	100%	0%	60°	100%	50%	100%	100%	14 (yellow)
Olive	#808000	50%	50%	0%	60°	100%	25%	100%	50%	6 (brown)
Lime	#00FF00	0%	100%	0%	120°	100%	50%	100%	100%	10 (high green); green
Green	#008000	0%	50%	0%	120°	100%	25%	100%	50%	2 (low green)
Aqua	#00FFFF	0%	100%	100%	180°	100%	50%	100%	100%	11 (high cyan); cyan
Teal	#008080	0%	50%	50%	180°	100%	25%	100%	50%	3 (low cyan)
Blue	#0000FF	0%	0%	100%	240°	100%	50%	100%	100%	9 (high blue)
Navy	#000080	0%	0%	50%	240°	100%	25%	100%	50%	1 (low blue)
Fuchsia	#FF00FF	100%	0%	100%	300°	100%	50%	100%	100%	13 (high magenta); magenta
Purple	#800080	50%	0%	50%	300°	100%	25%	100%	50%	5 (low magenta)

2. using colors in R

http://research.stowers-institute.org/efg/Report/UsingColorInR.pdf

Generate Distinct Colors for Your MATLAB Plots

http://blogs.mathworks.com/pick/2010/12/23/generate-distinct-colors-for-your-matlab-plots/

4. Color Palette and the 56 Excel ColorIndex Colors

interior	font	*HTML*	bgcolor=	*Red*<	Green	Blue	Color
Black	[Color 1]	#000000	#000000	0	0	0	[Black]
White	[Color 2]	#FFFFFF	#FFFFFF	255	255	255	[White]
Red	[Color 3]	#FF0000	#FF0000	255	0	0	[Red]
Green	[Color 4]	#00FF00	#00FF00	0	255	0	[Green]
Blue	[Color 5]	#0000FF	#0000FF	0	0	255	[Blue]
Yellow	[Color 6]	#FFFF00	#FFFF00	255	255	0	[Yellow]
Magenta	[Color 7]	#FF00FF	#FF00FF	255	0	255	[Magenta]
Cyan	[Color 8]	#00FFFF	#00FFFF	0	255	255	[Cyan]
[Color 9]	[Color 9]	#800000	#800000	128	0	0	[Color 9]
[Color 10]	[Color 10]	#008000	#008000	0	128	0	[Color 10]
[Color 11]	[Color 11]	#000080	#000080	0	0	128	[Color 11]
[Color 12]	[Color 12]	#808000	#808000	128	128	0	[Color 12]
[Color 13]	[Color 13]	#800080	#800080	128	0	128	[Color 13]
[Color 14]	[Color 14]	#008080	#008080	0	128	128	[Color 14]
[Color 15]	[Color 15]	#C0C0C0	#C0C0C0	192	192	192	[Color 15]
[Color 16]	[Color 16]	#808080	#808080	128	128	128	[Color 16]
[Color 17]	[Color 17]	#9999FF	#9999FF	153	153	255	[Color 17]
[Color 18]	[Color 18]	#993366	#993366	153	51	102	[Color 18]
[Color 19]	[Color 19]	#FFFFCC	#FFFFCC	255	255	204	[Color 19]
[Color 20]	[Color 20]	#CCFFFF	#CCFFFF	204	255	255	[Color 20]
[Color 21]	[Color 21]	#660066	#660066	102	0	102	[Color 21]
[Color 22]	[Color 22]	#FF8080	#FF8080	255	128	128	[Color 22]
[Color 23]	[Color 23]	#0066CC	#0066CC	0	102	204	[Color 23]
[Color 24]	[Color 24]	#CCCCFF	#CCCCFF	204	204	255	[Color 24]
[Color 25]	[Color 25]	#000080	#000080	0	0	128	[Color 25]
[Color 26]	[Color 26]	#FF00FF	#FF00FF	255	0	255	[Color 26]
[Color 27]	[Color 27]	#FFFF00	#FFFF00	255	255	0	[Color 27]
[Color 28]	[Color 28]	#00FFFF	#00FFFF	0	255	255	[Color 28]
[Color 29]	[Color 29]	#800080	#800080	128	0	128	[Color 29]
[Color 30]	[Color 30]	#800000	#800000	128	0	0	[Color 30]
[Color 31]	[Color 31]	#008080	#008080	0	128	128	[Color 31]
[Color 32]	[Color 32]	#0000FF	#0000FF	0	0	255	[Color 32]
[Color 33]	[Color 33]	#00CCFF	#00CCFF	0	204	255	[Color 33]
[Color 34]	[Color 34]	#CCFFFF	#CCFFFF	204	255	255	[Color 34]
[Color 35]	[Color 35]	#CCFFCC	#CCFFCC	204	255	204	[Color 35]
[Color 36]	[Color 36]	#FFFF99	#FFFF99	255	255	153	[Color 36]
[Color 37]	[Color 37]	#99CCFF	#99CCFF	153	204	255	[Color 37]
[Color 38]	[Color 38]	#FF99CC	#FF99CC	255	153	204	[Color 38]
[Color 39]	[Color 39]	#CC99FF	#CC99FF	204	153	255	[Color 39]
[Color 40]	[Color 40]	#FFCC99	#FFCC99	255	204	153	[Color 40]
[Color 41]	[Color 41]	#3366FF	#3366FF	51	102	255	[Color 41]
[Color 42]	[Color 42]	#33CCCC	#33CCCC	51	204	204	[Color 42]
[Color 43]	[Color 43]	#99CC00	#99CC00	153	204	0	[Color 43]
[Color 44]	[Color 44]	#FFCC00	#FFCC00	255	204	0	[Color 44]
[Color 45]	[Color 45]	#FF9900	#FF9900	255	153	0	[Color 45]
[Color 46]	[Color 46]	#FF6600	#FF6600	255	102	0	[Color 46]
[Color 47]	[Color 47]	#666699	#666699	102	102	153	[Color 47]
[Color 48]	[Color 48]	#969696	#969696	150	150	150	[Color 48]
[Color 49]	[Color 49]	#003366	#003366	0	51	102	[Color 49]
[Color 50]	[Color 50]	#339966	#339966	51	153	102	[Color 50]
[Color 51]	[Color 51]	#003300	#003300	0	51	0	[Color 51]
[Color 52]	[Color 52]	#333300	#333300	51	51	0	[Color 52]
[Color 53]	[Color 53]	#993300	#993300	153	51	0	[Color 53]
[Color 54]	[Color 54]	#993366	#993366	153	51	102	[Color 54]
[Color 55]	[Color 55]	#333399	#333399	51	51	153	[Color 55]
[Color 56]	[Color 56]	#333333	#333333	51	51	51	[Color 56]
Excel only recognizes names for Color 1 through 8 (Black, White, Red, Green, Blue, Yellow, Magenta, and Cyan). The colors 1-16 are widely understood color names from the VGA color palette. Of the 56 colors only 40 colors appear on the palette. The 40 colors names indicated on the Excel color palette (see below) are for descriptive purposes only.

2013年7月2日星期二

parsing arguements in a R script

http://ggorjan.blogspot.ca/2013/07/parse-arguments-of-r-script.html?utm_source=feedburner&utm_medium=feed&utm_campaign=Feed:+GregorGorjancsBlog+(Gregor+Gorjanc's+blog)

R can be used also as a scripting tool. We just need to add shebang in the first line of a file (script):

#!/usr/bin/Rscript

and then the R code should follow.

Often we want to pass arguments to such a script, which can be collected in the script by the commandArgs() function. Then we need to parse the arguments and conditional on them do something. I came with a rather general way of parsing these arguments using simply these few lines:

## Collect arguments
args <- a="" href="http://inside-r.org/r-doc/base/commandArgs" style="color: #b81909; text-decoration: none;">commandArgs

(TRUE) ## Default setting when no arguments passed if(length(args) < 1) { args <- a="" href="http://inside-r.org/r-doc/base/c" style="color: #b81909; text-decoration: none;">c("--help") } ## Help section if("--help" %in% args) { cat(" The R Script Arguments: --arg1=someValue - numeric, blah blah --arg2=someValue - character, blah blah --arg3=someValue - logical, blah blah --help - print this text Example: ./test.R --arg1=1 --arg2="output.txt" --arg3=TRUE \n\n") q(save="no") } ## Parse arguments (we expect the form --arg=value) parseArgs <- a="" href="http://inside-r.org/r-doc/base/function" style="color: #b81909; text-decoration: none;">function(x) strsplit(sub("^--", "", x), "=") argsDF <- a="" href="http://inside-r.org/r-doc/base/as.data.frame" style="color: #b81909; text-decoration: none;">as.data.frame(do.call("rbind", parseArgs(args))) argsL <- a="" href="http://inside-r.org/r-doc/base/as.list" style="color: #b81909; text-decoration: none;">as.list(as.character(argsDF$V2)) names(argsL) <- argsdf="" nbsp="" span="" style="font-style: italic;">## Arg1 default if(is.null(args$arg1)) { ## do something } ## Arg2 default if(is.null(args$arg2)) { ## do something } ## Arg3 default if(is.null(args$arg3)) { ## do something }

## ... your code here ...

Created by Pretty R at inside-R.org

It is some work, but I find it pretty neat and use it for quite a while now. I do wonder what others have come up for this task. I hope I did not miss some very general solution.

2013年6月27日星期四

distinguish effects of geographic and ecological isolation on differentiation

Populations can be genetically isolated by both geographic distance and by differences in their ecology or environment that decrease the rate of successful migration. Empirical studies often seek to investigate the relationship between genetic differentiation and some ecological variable(s) while accounting for geographic distance, but common approaches to this problem (such as the partial Mantel test) have a number of drawbacks. In this article, we present a Bayesian method that enables users to quantify the relative contributions of geographic distance and ecological distance to genetic differentiation between sampled populations or individuals. We model the allele frequencies in a set of populations at a set of unlinked loci as spatially correlated Gaussian processes, in which the covariance structure is a decreasing function of both geographic and ecological distance. Parameters of the model are estimated using a Markov chain Monte Carlo algorithm. We call this method Bayesian Estimation of Differentiation in Alleles by Spatial Structure and Local Ecology (BEDASSLE), and have implemented it in a user-friendly format in the statistical platform R. We demonstrate its utility with a simulation study and empirical applications to human and teosinte datasets.

http://arxiv.org/abs/1302.3274

Piping with samtools, bwa and bedtools

http://www.biostars.org/p/43677/

In this tutorial I will introduce some concepts related to unix piping. Piping is a very useful feature to avoid creation of intermediate use once files. It is assumed that bedtools, samtools, and bwa are installed.

Lets begin with a typical command to do paired end mapping with bwa: (./ means look in current directory only)

#-t 4 is for using 4 threads/cores
bwa aln -t 4 ./hg19.fasta ./s1_1.fastq > ./s1_1.sai
bwa aln -t 4 ./hg19.fasta ./s1_2.fastq > ./s1_2.sai
bwa sampe ./hg19.fasta ./s1_1.sai ./s1_2.sai ./s1_1.fastq ./s1_2.fastq > s1.sam

Supposed we wish to compress sam to bam, sort, remove duplicates, and create a bed file.

samtools view -Shu s1.sam > s1.bam
samtools sort s1.bam s1_sorted
samtools rmdup -s s1_sorted.bam s1_sorted_nodup.bam
bamToBed -i s1_sorted_nodup.bam > s1_sorted_nodup.bed

This workflow above creates many files that are only used once (such as s1.bam) and we can use the unix pipe utility to reduce the number intermediate files created. The pipe function is the character | and what it does is take the output from one command and sets it as input for next command (assuming next command knows how to deal with this input).

For example, when you type in head myfile.txt the command 'head' will read first 10 lines of myfile.txt and output it (by default) to the display (stdout) so you will see the first 10 lines printed on your screen. However, if you were to do something like this:

head myfile.txt | wc -l

you will instead get 10 printed out on your screen. What the pipe command did was take the output of head command and used it as input for the wordcount (wc) command. wc -l command will count the number of lines passed in (which is 10 since head by default prints the first 10 lines). An equivalent command would be

head myfile.txt > myfile_top10_lines.txt
wc -l myfile_top10_lines.txt

By using the pipe command we are able to eliminate the intermediate file that was generated (myfile_top10_lines.txt).

Below is an example of how pipe command can be used in samtools

bwa aln -t 4 ./hg19.fasta ./s1_1.fastq > ./s1_1.sai
bwa aln -t 4 ./hg19.fasta ./s1_2.fastq > ./s1_2.sai
bwa sampe ./hg19.fasta ./s1_1.sai ./s1_2.sai ./s1_1.fastq ./s1_2.fastq | \
samtools view -Shu - | \
samtools sort - - | \
samtools rmdup -s - - | \
tee s1_sorted_nodup.bam | \
bamToBed > s1_sorted_nodup.bed

A couple of things I wanted to point out here:

The \ symbol after the pipe symbol tells the terminal that I want to keep writing command on the next line, that way you can have multiline commands instead of super-long command on one line.
The tee command takes the information being piped and creates a file while passing the data along the pipe, this command is useful if you want to save an intermediate file but still need to do additional processing. (Think of a T shaped drain pipe)
The - symbol in samtools it to tell the samtools program to take the input from pipe
It is not recommended you run this command all at once like I showed above if you have a giant bam file because sort will take forever. Running things in a pipe can be more space efficient since you are not storing intermediate files but can create issues/bottlenecks elsewhere. If problems arise, run the pieces without piping to troubleshoot.
The two bwa aln commands can be run in parallel (for example: on another machine) before the sampe command

Other common uses of pipes that I have not shown above it to pipe output into snp calling, you can find some examples here: (NOTE: pileupcommand is depreciated, the new way to call SNPs using samtools using mpileup command)

samtools merge - *.bam | tee merged.bam | samtools rmdup - - | tee rmdup.bam | samtools mpileup - uf ./hg19.fasta - | bcftools view -bvcg - | gzip > var.raw.bcf.gz

Compressing files with gzip can save a lot of space, but makes these files non-human readable (if you use the head or less command it will give you gibberish). To view these files you can use zless and zcat command.

zless var.raw.bcf.gz #similar to gzip -cd var.raw.bcf.gz | less
zcat var.raw.bcf.gz | head

In case you are wondering if there is a way to run bwa sampe without storing intermediate sai files the command below will do this:

bwa sampe ./hg19.fasta <(bwa aln -t 4 ./hg19.fasta ./s1_1.fastq) <(bwa aln -t 4 ./hg19.fasta ./s1_2.fastq) ./s1_1.fastq ./s1_2.fastq | samtools view -Shb /dev/stdin > s1.bam

Notice here that I used the -b command instead of -u command in samtools view. -u is useful for piping since it spits out an uncompressed bam so you don't have to waste CPU cycles compressing/decompressing at every step, however, if you have a file you want to store, it would be better to compress it with -b. Additionally, I used /dev/stdin instead of using the - symbol, they both mean the same thing (afaik) and I wanted to point it out that /dev/stdin might be a more explicit way of showing things. (ie samtools view -Shu /dev/stdin )

The Elements of Bioinformatics

http://elements.eaglegenomics.com/

2013年6月16日星期日

X11 forwarding over SSH: run remote graphical app and display locally

http://linuxcommando.blogspot.ca/2013/05/x11-forwarding-over-ssh-run-remote.html

In the modern networked environment, we often wish to run an application on a remote host while we are comfortably logged in on our local computer.

Assuming both machines are Linux-based, and the application runs on the graphical X desktop, the following approaches come to mind:

VNC
X11 forwarding over SSH

This article focuses only on X11 forwarding. X11 forwarding over SSH enables you to run a remote X app and display it locally, with traffic between the 2 hosts encrypted by SSH.

For X11 forwarding over SSH to work, both the SSH client and SSH server must be properly configured.

X11 forwarding must be enabled on The SSH server side. This is the machine where the application resides. To enable the feature, make sure the X11 configuration file /etc/ssh/sshd_config on the server contains this line:

X11Forwarding yes

If you edit the said file, you need to restart the sshd daemon for the change to take effect.

On Debian or Ubuntu systems, you restart the SSH daemon like this:

$ sudo service ssh restart
[ ok ] Restarting OpenBSD Secure Shell server: sshd.
$

On the ssh client side, you need to run SSH command with the proper parameters. For instance, suppose you want to run the xclockapplication on the remote SSH server and have it displayed back on the local client.

$ ssh -fX peter@192.168.1.112 xclock 
peter@192.168.1.112's password: 
$

The -X parameter allows an one-off X11 forwarding session.

The -f parameter instructs the SSH client to go to the background just before xclock is run.

If you want to permanently enable X11 forwarding for an user, insert this line in the user's own ~/.ssh/config file on the local host.

ForwardX11 yes

With X11 forwarding permanently enabled for the client, you can leave out the -X parameter:

$ ssh -f peter@192.168.1.112 xclock 
peter@192.168.1.112's password: 
$

If X11 forwarding is not enabled on the SSH server, any attempt to tunnel X11 will fail with the following error message:

$ ssh -X peter@192.168.1.112 xclock 
peter@192.168.1.112's password: 
X11 forwarding request failed on channel 0
Error: Can't open display: 
$

If X11 forwarding is properly enabled on the server side, you will see a nice looking clock displayed on your local screen.

2013年5月26日星期日

Framework for evaluating variant detection methods: comparison of aligners and callers

I’ll show a 12 way comparison between 2 different aligners (novoalign and bwa mem), 2 different post-alignment preparation methods (GATK best practices and the Marth lab’s gkno pipeline), and 3 different variant callers (GATK UnifiedGenotyper, GATK HaplotypeCaller, and FreeBayes). This allows comparison of openly available methods (bwa mem, gkno preparation, and FreeBayes) with those that require licensing (novoalign, GATK’s variant callers). I’ll also describe bcbio-nextgen, the fully automated open-source pipeline used for variant calling and evaluation, which allows others to easily bring this methodology into their own work and extend this analysis.

2013年5月25日星期六

How does your bash prompt look like

http://kevin-gattaca.blogspot.ca/

http://www.cyberciti.biz/tips/howto-linux-unix-bash-shell-setup-prompt.html

http://www.maketecheasier.com/8-useful-and-interesting-bash-prompts/2009/09/04

2013年5月22日星期三

Learn Perl in about 2 hours 30 minutes

http://qntm.org/files/perl/perl.html

2013年5月18日星期六

Multiple populations of artemisinin-resistant Plasmodium falciparum in Cambodia

http://www.nature.com/ng/journal/vaop/ncurrent/full/ng.2624.html

We describe an analysis of genome variation in 825 P. falciparum samples from Asia and Africa that identifies an unusual pattern of parasite population structure at the epicenter of artemisinin resistance in western Cambodia. Within this relatively small geographic area, we have discovered several distinct but apparently sympatric parasite subpopulations with extremely high levels of genetic differentiation. Of particular interest are three subpopulations, all associated with clinical resistance to artemisinin, which have skewed allele frequency spectra and high levels of haplotype homozygosity, indicative of founder effects and recent population expansion. We provide a catalog of SNPs that show high levels of differentiation in the artemisinin-resistant subpopulations, including codon variants in transporter proteins and DNA mismatch repair proteins. These data provide a population-level genetic framework for investigating the biological origins of artemisinin resistance and for defining molecular markers to assist in its elimination.