Dropbox is having sync problem on my Snow Leopad laptop. There is a triangle but the uploading is stuck.
I deleted ./dropboxcache and the problem persists.
Quit and restart Dropbox by 'open Dropbox.app': Starting at 19:52:25 -> uploading up-traingle 19:54:10 -> still stuck.
19:57:23, $ mv Dropbox.app Dropbox.app.old
Download Dropbox.1.6.13.dmg
20:00:00, Reinstall Dropbox
Icon appears at up-panel with a spinning double-arrow circle, "indexing, connecting, uploading 291 files .."
20:02: it shows "Uploading 291 files ... "
By 20:10, it is still showing uploading 291 files
Maybe it is the 2G encrypted drive that makes the sync looks stuck. But the previous icon did not show a spining cirlc, it was flashing square). In any case, maybe it is working now. I decided to remove this 2G encryptic disk.
20:28, Uploading 479 files, indexing 53 files.
20:31, uploading 523 files, connecting ...
20:35, "Unlink" the laptop on Dropbox website. (This seems to solves the problem)
20:37, open Drobpx.app
Requesting link .. ...
20:39 chose previous setting
20:39, indexing 591 files, downloading file list
20:42, Uploading 2694 files, indexing 14,771 files, downloading file list
20:44, Uploading 4156 files, indexing 11,582 files, ...
20:56, Uploading 241 files, indexing 8,590 files ...
Around 22:07, all of the files are finally synced.
21:13, I found file '.dropbox' in the $Home directory. Maybe I should examine this directory next time.
I realizing the slow Spelman network also plays a role here. The uploading speed is 1,288KB/sec. Because the encrypted disk is a 2G file, updating it takes very long time.
I spent 1.5 hours on fixing Dropbox sync problem.
Thursday, January 10, 2013
Tuesday, January 8, 2013
A look at identical twins from the systems biology perspective
After I listened to a BBC podcast about Poincaré and bifurcation in dynamic systems, I had an interesting question. A key point of dynamic systems is the so-called butterfly effect - small changes in initial conditions can lead to large changes in the future. I thought about the identical twins in biology. They always look physically similar, presumably because structural development has a lot of constraints, its layout was set during the early stages of development, and there was not enough time for changes (noises) to propagate. The personality differences between the identical twins could be understood by the butterfly effect because neurological and emotional development took decades and allow environmental input signals to be amplified in the systems. I do not have expert understanding of both developmental biology and dynamics systems, and I wonder what other factors and principles are behind this. This could be an interesting story for discussion in an introductory course on systems biology or developmental biology.
After I posted the initial paragraph, a friend reminded me of Stuart Kauffman's argument that complex systems ought to adapt to the border between chaos and order during evolution. Kauffman's argument indeed makes sense. In fact, if I can speak freely, bio-diversity might have root in the butterfly effect. Although many species have similar sets of orthologous genes, their phenotypic diversity could be explained by the different trajectories caused by the changes in the initial conditions such as the divergent genotypes. A good example maybe chimpanzees and humans: Small genotypic differences (initial conditions) are amplified to drastically different phenotypes in non-structural traits.
References:
After I posted the initial paragraph, a friend reminded me of Stuart Kauffman's argument that complex systems ought to adapt to the border between chaos and order during evolution. Kauffman's argument indeed makes sense. In fact, if I can speak freely, bio-diversity might have root in the butterfly effect. Although many species have similar sets of orthologous genes, their phenotypic diversity could be explained by the different trajectories caused by the changes in the initial conditions such as the divergent genotypes. A good example maybe chimpanzees and humans: Small genotypic differences (initial conditions) are amplified to drastically different phenotypes in non-structural traits.
References:
- BBC podcast on the story of Henri Poincare and butterfly effect
- Stuart Kauffman, The origins of order and self organization selection
Friday, January 4, 2013
Compare tree topologies using the ape package in R
The following R code read a newick file with multiple trees, and calculate the topological distance between all pairs.
require(ape)
help(package=ape)
trees = read.tree("treesall.nwk") #list of 10
trees[[1]]
trees[[2]]
dist.topo( trees[[1]], trees[[10]])
d.tree = matrix( nrow=10, ncol=10 )
for ( i in 1:10 ) {
for (j in i:10) {
print ( c(i, j))
d.tree[i,j] = dist.topo( trees[[i]], trees[[j]] )
}
}
d.tree
An output from the 10 tree for Bacillus essential genes are:
> d.tree
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 2 2 4 6 2 4 4 4 4
[2,] NA 0 4 6 6 4 2 6 6 6
[3,] NA NA 0 2 4 4 6 6 6 6
[4,] NA NA NA 0 2 4 6 6 6 8
[5,] NA NA NA NA 0 6 6 8 8 10
[6,] NA NA NA NA NA 0 2 4 4 6
[7,] NA NA NA NA NA NA 0 6 6 8
[8,] NA NA NA NA NA NA NA 0 2 8
[9,] NA NA NA NA NA NA NA NA 0 8
[10,] NA NA NA NA NA NA NA NA NA 0
require(ape)
help(package=ape)
trees = read.tree("treesall.nwk") #list of 10
trees[[1]]
trees[[2]]
dist.topo( trees[[1]], trees[[10]])
d.tree = matrix( nrow=10, ncol=10 )
for ( i in 1:10 ) {
for (j in i:10) {
print ( c(i, j))
d.tree[i,j] = dist.topo( trees[[i]], trees[[j]] )
}
}
d.tree
An output from the 10 tree for Bacillus essential genes are:
> d.tree
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 2 2 4 6 2 4 4 4 4
[2,] NA 0 4 6 6 4 2 6 6 6
[3,] NA NA 0 2 4 4 6 6 6 6
[4,] NA NA NA 0 2 4 6 6 6 8
[5,] NA NA NA NA 0 6 6 8 8 10
[6,] NA NA NA NA NA 0 2 4 4 6
[7,] NA NA NA NA NA NA 0 6 6 8
[8,] NA NA NA NA NA NA NA 0 2 8
[9,] NA NA NA NA NA NA NA NA 0 8
[10,] NA NA NA NA NA NA NA NA NA 0
CONSEL result interpretation and reformat using perl
From "program.pdf" in CONSEL package,
2.4 catpv
The p-values calculated by consel are stored in the pv ¯le, and its contents are
shown by the command catpv. For example,
catpv foo
reads foo.pv and shows the contents like below.
# reading foo.pv
# rank item obs au np | bp kh sh wkh wsh |
# 1 1 -2.7 0.789 0.575 | 0.579 0.639 0.944 0.639 0.948 |
# 2 3 2.7 0.516 0.318 | 0.312 0.361 0.799 0.361 0.791 |
# 3 2 7.4 0.114 0.037 | 0.036 0.122 0.575 0.122 0.422 |
# 4 5 17.6 0.076 0.014 | 0.013 0.044 0.178 0.044 0.210 |
'au' The p-value of the approximately unbiased test. This is the main result of
consel, while rest of the p-values in the table are supplementary. It is derived
from the theory of the signed distance and the curvature of the boundary
as explained later. Let us denote it AUi for the item i. We may reject the
possibility that item i has the largest expected value of Yi when AUi < 0.05
at the significance level 0.05.
In its foo example, the trees are ranked based on Y_i test statistic. Difference of Y_i, \delta Y_i is shown in the report. I noticed that au p-value are often the highest for the top ranked tree. When a tree's au pvalue is less than 0.05, this tree can be considered as rejected. Y_i statistic and p-value orders can give different orders, and this can be see in the foo example.
The following perl code can convert CONSEL report to a tab-delimited file for R.
use strict; use warnings;
use lib '/Users/hongqin/lib/perl'; use Util;
my $debug= 2 ;
my $inclusterfile = '_test.csv';
my $sourcedir = '~coat/ortholog.analysis/essetial.mcl.rbh.090208/bacillus.phylip.I3.1/wkdir';
open (IN,"<$inclusterfile"); my @lines = <IN>; close (IN); shift @lines;
my @clusterdirs = ();
foreach my $line (@lines) {
chomp $line;
my ($num, $bg, @res) = split( /\t/, $line);
print "[$num][$bg]------";
push (@clusterdirs, $bg);
}
open ( OUT, ">_conselreport.34essenGene.2012Jan4.tab" );
print OUT "BG\trank\titem\tobs\tau\tnp\tbp\tpp\tkh\tsh\twkh\twsh\n";
my $count = 0;
foreach my $mydir ( @clusterdirs ) {
chomp $mydir;
if ( !( -e "wkdir/$mydir/$mydir.consel.report.txt" ) ) {
print "Error: $mydir has no consel report";
} else {
$count ++; print "Cluster:$count I am now working on [$mydir]::\n";
open (IN, "<wkdir/$mydir/$mydir.consel.report.txt"); @lines = <IN>; close (IN);
shift @lines; shift @lines; shift @lines; #remove the first three lines
pop @lines; #remove the last empty line
foreach my $line (@lines) {
$line =~ s/^#\s+//o;
$line =~ s/\|//g;
chomp $line;
my @els = split( /\s+/, $line);
print OUT $mydir;
foreach my $el (@els) { print OUT "\t$el"; }
print OUT "\n";
} }
}
close (OUT);
exit;
2.4 catpv
The p-values calculated by consel are stored in the pv ¯le, and its contents are
shown by the command catpv. For example,
catpv foo
reads foo.pv and shows the contents like below.
# reading foo.pv
# rank item obs au np | bp kh sh wkh wsh |
# 1 1 -2.7 0.789 0.575 | 0.579 0.639 0.944 0.639 0.948 |
# 2 3 2.7 0.516 0.318 | 0.312 0.361 0.799 0.361 0.791 |
# 3 2 7.4 0.114 0.037 | 0.036 0.122 0.575 0.122 0.422 |
# 4 5 17.6 0.076 0.014 | 0.013 0.044 0.178 0.044 0.210 |
'au' The p-value of the approximately unbiased test. This is the main result of
consel, while rest of the p-values in the table are supplementary. It is derived
from the theory of the signed distance and the curvature of the boundary
as explained later. Let us denote it AUi for the item i. We may reject the
possibility that item i has the largest expected value of Yi when AUi < 0.05
at the significance level 0.05.
In its foo example, the trees are ranked based on Y_i test statistic. Difference of Y_i, \delta Y_i is shown in the report. I noticed that au p-value are often the highest for the top ranked tree. When a tree's au pvalue is less than 0.05, this tree can be considered as rejected. Y_i statistic and p-value orders can give different orders, and this can be see in the foo example.
The following perl code can convert CONSEL report to a tab-delimited file for R.
use strict; use warnings;
use lib '/Users/hongqin/lib/perl'; use Util;
my $debug= 2 ;
my $inclusterfile = '_test.csv';
my $sourcedir = '~coat/ortholog.analysis/essetial.mcl.rbh.090208/bacillus.phylip.I3.1/wkdir';
open (IN,"<$inclusterfile"); my @lines = <IN>; close (IN); shift @lines;
my @clusterdirs = ();
foreach my $line (@lines) {
chomp $line;
my ($num, $bg, @res) = split( /\t/, $line);
print "[$num][$bg]------";
push (@clusterdirs, $bg);
}
open ( OUT, ">_conselreport.34essenGene.2012Jan4.tab" );
print OUT "BG\trank\titem\tobs\tau\tnp\tbp\tpp\tkh\tsh\twkh\twsh\n";
my $count = 0;
foreach my $mydir ( @clusterdirs ) {
chomp $mydir;
if ( !( -e "wkdir/$mydir/$mydir.consel.report.txt" ) ) {
print "Error: $mydir has no consel report";
} else {
$count ++; print "Cluster:$count I am now working on [$mydir]::\n";
open (IN, "<wkdir/$mydir/$mydir.consel.report.txt"); @lines = <IN>; close (IN);
shift @lines; shift @lines; shift @lines; #remove the first three lines
pop @lines; #remove the last empty line
foreach my $line (@lines) {
$line =~ s/^#\s+//o;
$line =~ s/\|//g;
chomp $line;
my @els = split( /\s+/, $line);
print OUT $mydir;
foreach my $el (@els) { print OUT "\t$el"; }
print OUT "\n";
} }
}
close (OUT);
exit;
Wednesday, January 2, 2013
First try to enumerate truncated trees for bacilli
import warnings
import re
from Bio import Phylo
stree = "((Bha,Bcl), (Bpu,(Bli,(Bam,(Bmo,Bsu)))), (Bwe,(Bce,(Ban,Bth))))"
stree2 = re.sub( r'[\(\)\s+]', '', stree)
taxa = re.split( ',', stree2)
outgrouptaxa = taxa[0:2]
bsutaxa = taxa[2:7]
bcetaxa = taxa[7:12]
bacillustaxa = bsutaxa + bcetaxa
#tmpfl = "/tmp/_tmpfile2013Jan2.nwk"
#tmpFH = open( tmpfl, 'w')
#tmpFH.write( stree + ";\n" )
#tmpFH.close()
#mastertree = Phylo.read(tmpfl, 'newick')
def lookup_by_names(tree):
names = {}
for clade in tree.get_terminals():
if clade.name:
if clade.name in names:
raise ValueError("Duplicate key: %s" % clade.name)
names[clade.name] = clade
return names
from StringIO import StringIO
#mastertree = Phylo.read(StringIO(stree), 'newick')
#names = lookup_by_names(mastertree)
#single leaf-node deletion
for node in bacillustaxa:
tree = Phylo.read(StringIO(stree), 'newick')
names = lookup_by_names(tree)
tree.prune(names[node])
Phylo.write( tree, 'sandbox/-' + node + '.nwk', 'newick')
#double leaf node deletion
for i in range(len(bacillustaxa)):
for j in range( (i+1), len(bacillustaxa)):
print i, j
tree = Phylo.read(StringIO(stree), 'newick')
names = lookup_by_names(tree)
tree.prune(names[bacillustaxa[i]])
tree.prune(names[bacillustaxa[j]])
Phylo.write( tree, 'sandbox/_' + bacillustaxa[i] + bacillustaxa[j]+'.nwk', 'newick')
$ pwd
~coat.protein/ortholog.analysis/coatpaml/python/
This is not bad for a one-hour coding on a first try. I should try to write multiple trees in a single file. I should try a generic loop or more likely, a recursive call.
import re
from Bio import Phylo
stree = "((Bha,Bcl), (Bpu,(Bli,(Bam,(Bmo,Bsu)))), (Bwe,(Bce,(Ban,Bth))))"
stree2 = re.sub( r'[\(\)\s+]', '', stree)
taxa = re.split( ',', stree2)
outgrouptaxa = taxa[0:2]
bsutaxa = taxa[2:7]
bcetaxa = taxa[7:12]
bacillustaxa = bsutaxa + bcetaxa
#tmpfl = "/tmp/_tmpfile2013Jan2.nwk"
#tmpFH = open( tmpfl, 'w')
#tmpFH.write( stree + ";\n" )
#tmpFH.close()
#mastertree = Phylo.read(tmpfl, 'newick')
def lookup_by_names(tree):
names = {}
for clade in tree.get_terminals():
if clade.name:
if clade.name in names:
raise ValueError("Duplicate key: %s" % clade.name)
names[clade.name] = clade
return names
from StringIO import StringIO
#mastertree = Phylo.read(StringIO(stree), 'newick')
#names = lookup_by_names(mastertree)
#single leaf-node deletion
for node in bacillustaxa:
tree = Phylo.read(StringIO(stree), 'newick')
names = lookup_by_names(tree)
tree.prune(names[node])
Phylo.write( tree, 'sandbox/-' + node + '.nwk', 'newick')
#double leaf node deletion
for i in range(len(bacillustaxa)):
for j in range( (i+1), len(bacillustaxa)):
print i, j
tree = Phylo.read(StringIO(stree), 'newick')
names = lookup_by_names(tree)
tree.prune(names[bacillustaxa[i]])
tree.prune(names[bacillustaxa[j]])
Phylo.write( tree, 'sandbox/_' + bacillustaxa[i] + bacillustaxa[j]+'.nwk', 'newick')
$ pwd
~coat.protein/ortholog.analysis/coatpaml/python/
This is not bad for a one-hour coding on a first try. I should try to write multiple trees in a single file. I should try a generic loop or more likely, a recursive call.
Monday, December 31, 2012
*** Git/Github repository for Bacillus spore coat project
$ cd /BacillusSporeCoat-master
$ git init
Initialized empty Git repository in /Users/hongqin/github/BacillusSporeCoat-master/.git/
$ git add _coat.profile.for.submission.Dec31.2012.xlsx
$ git commit -m "coat profile"
[master (root-commit) fe3211c] coat profile
1 file changed, 0 insertions(+), 0 deletions(-)
create mode 100644 _coat.profile.for.submission.Dec31.2012.xlsx
$ git init
Initialized empty Git repository in /Users/hongqin/github/BacillusSporeCoat-master/.git/
$ git add _coat.profile.for.submission.Dec31.2012.xlsx
$ git commit -m "coat profile"
[master (root-commit) fe3211c] coat profile
1 file changed, 0 insertions(+), 0 deletions(-)
create mode 100644 _coat.profile.for.submission.Dec31.2012.xlsx
$ git remote add origin https://github.com/hongqin/BacillusSporeCoat.git
$ git add README.md
$ git commit -m "coat profile 2"
[master 5aac99d] coat profile 2
1 file changed, 4 insertions(+)
create mode 100644 README.md
ace:BacillusSporeCoat-master hongqin$ git remote add origin https://github.com/hongqin/BacillusSporeCoat.git
fatal: remote origin already exists.
[master 5aac99d] coat profile 2
1 file changed, 4 insertions(+)
create mode 100644 README.md
ace:BacillusSporeCoat-master hongqin$ git remote add origin https://github.com/hongqin/BacillusSporeCoat.git
fatal: remote origin already exists.
$ git pull https://github.com/hongqin/BacillusSporeCoat.git -m "test"
$ git remote add origin https://github.com/hongqin/BacillusSporeCoat.git
fatal: remote origin already exists.
fatal: remote origin already exists.
$ git push origin master
Counting objects: 7, done.
Delta compression using up to 2 threads.
Compressing objects: 100% (6/6), done.
Writing objects: 100% (6/6), 21.62 KiB, done.
Total 6 (delta 1), reused 0 (delta 0)
To https://github.com/hongqin/BacillusSporeCoat.git
5d270ff..1b82f44 master -> master
Counting objects: 7, done.
Delta compression using up to 2 threads.
Compressing objects: 100% (6/6), done.
Writing objects: 100% (6/6), 21.62 KiB, done.
Total 6 (delta 1), reused 0 (delta 0)
To https://github.com/hongqin/BacillusSporeCoat.git
5d270ff..1b82f44 master -> master
This works.
Saturday, December 29, 2012
Install Git and GitHub for snow leopard
I have a snow leopard laptop.
From http://git-scm.com/download/mac, it automatically downloaded "git-1.8.0.2-3-intel-universal-snow-leopard.dmg". Double-click this file generated a "Git 1.8.0.2 Snow Leopard Intel Universal" disk, which contains four files: 1) uninstall.sh, 2) README, 3) 'setup git PATH for non-terminal programs.sh, 4) git-1.8.0.2-intel-universal-snow-leopard.pkg. After running the fourth file, a successful install message occurred.
I then found out that the free GUI for Git, GitHub, has stopped support for snow leopard. This leaves me the only option to follow the command-line version on https://help.github.com/articles/set-up-git.
$ git config --global user.name "hongqin"
From http://git-scm.com/download/mac, it automatically downloaded "git-1.8.0.2-3-intel-universal-snow-leopard.dmg". Double-click this file generated a "Git 1.8.0.2 Snow Leopard Intel Universal" disk, which contains four files: 1) uninstall.sh, 2) README, 3) 'setup git PATH for non-terminal programs.sh, 4) git-1.8.0.2-intel-universal-snow-leopard.pkg. After running the fourth file, a successful install message occurred.
I then found out that the free GUI for Git, GitHub, has stopped support for snow leopard. This leaves me the only option to follow the command-line version on https://help.github.com/articles/set-up-git.
$ git config --global user.name "hongqin"
$ git config --global user.email "myworkemail@gmail.com"
$ curl -s -O \
http://github-media-downloads.s3.amazonaws.com/osx/git-credential-osxkeychain
$ chmod u+x git-credential-osxkeychain
$ cd /opt/local/bin #I saw osxkeychain already here
$ git config --global credential.helper osxkeychain
By now, I have finished all the steps on https://help.github.com/articles/set-up-git.
I then proceed to created a 'Public' repository, named "sandbox". I chose a README option. After hitting the green "create repository" button, a new page appeared in about one second.
When I clicked "Clone in Mac", it took me back to GitHub Mac version for 10.7+. When I clicked "ZIP", it started download "sandbox-master.zip" into local Download. This zip file contains the README.md file and the sandbox-master directory structure. I then noticed that I was in the 'branch:master' mode.
To test how GitHub works locally, I created a hello world file in local sandbox-master.
$cd sandbox-master
$git init
Initialized empty Git repository in /Users/hongqin/github/sandbox-master/.git/
$ git add helloworld.txt
$ git commit -m "first local commit"
[master (root-commit) 6654d46] first local commit
1 file changed, 1 insertion(+)
create mode 100644 helloword.txt
To test how GitHub works locally, I created a hello world file in local sandbox-master.
$cd sandbox-master
$git init
Initialized empty Git repository in /Users/hongqin/github/sandbox-master/.git/
$ git add helloworld.txt
$ git commit -m "first local commit"
[master (root-commit) 6654d46] first local commit
1 file changed, 1 insertion(+)
create mode 100644 helloword.txt
$ git remote add origin https://github.com/hongqin/sandbox.git
$ git push origin master
Username for 'https://github.com': hongqin
Password for 'https://hongqin@github.com':
To https://github.com/hongqin/sandbox.git
! [rejected] master -> master (non-fast-forward)
error: failed to push some refs to 'https://github.com/hongqin/sandbox.git'
hint: Updates were rejected because the tip of your current branch is behind
hint: its remote counterpart. Merge the remote changes (e.g. 'git pull')
hint: before pushing again.
hint: See the 'Note about fast-forwards' in 'git push --help' for details.
This non-fast-forward error seems to be a common error:
https://help.github.com/articles/dealing-with-non-fast-forward-errors
Username for 'https://github.com': hongqin
Password for 'https://hongqin@github.com':
To https://github.com/hongqin/sandbox.git
! [rejected] master -> master (non-fast-forward)
error: failed to push some refs to 'https://github.com/hongqin/sandbox.git'
hint: Updates were rejected because the tip of your current branch is behind
hint: its remote counterpart. Merge the remote changes (e.g. 'git pull')
hint: before pushing again.
hint: See the 'Note about fast-forwards' in 'git push --help' for details.
This non-fast-forward error seems to be a common error:
https://help.github.com/articles/dealing-with-non-fast-forward-errors
$ git push --force origin master
Counting objects: 3, done.
Writing objects: 100% (3/3), 230 bytes, done.
Total 3 (delta 0), reused 0 (delta 0)
To https://github.com/hongqin/sandbox.git
+ c8d6cad...6654d46 master -> master (forced update)
Counting objects: 3, done.
Writing objects: 100% (3/3), 230 bytes, done.
Total 3 (delta 0), reused 0 (delta 0)
To https://github.com/hongqin/sandbox.git
+ c8d6cad...6654d46 master -> master (forced update)
I then refreshed GitHub sandbox web page, and found "helloworld.txt" there. But "README.md" disappeared.This may explain my 'non-fast-forward error'.
$ git add README.txt
$ git commit -m "second local commit"
[master bdcefa8] second local commit
1 file changed, 2 insertions(+)
create mode 100644 README.txt
$ git add README.txt
$ git commit -m "second local commit"
[master bdcefa8] second local commit
1 file changed, 2 insertions(+)
create mode 100644 README.txt
$ git remote add origin https://github.com/hongqin/sandbox.git
fatal: remote origin already exists.
$ git push origin master
Counting objects: 4, done.
Delta compression using up to 2 threads.
Compressing objects: 100% (2/2), done.
Writing objects: 100% (3/3), 288 bytes, done.
Total 3 (delta 0), reused 0 (delta 0)
To https://github.com/hongqin/sandbox.git
6654d46..bdcefa8 master -> master
Delta compression using up to 2 threads.
Compressing objects: 100% (2/2), done.
Writing objects: 100% (3/3), 288 bytes, done.
Total 3 (delta 0), reused 0 (delta 0)
To https://github.com/hongqin/sandbox.git
6654d46..bdcefa8 master -> master
I then refreshed my GitHub webpage and see both files.
Oddly, there is still "README.md" file, but I do not see this file in the downloaed ZIP archive. I leave it for now.
NPR food-for-thought blog on lactose tolerance.
NPR has a blog on the evolution of lactose tolerance: An Evolutionary Whodunit: How Did Humans Develop Lactose Tolerance? by Helen Thompson. This blog discuss a review on lactase persistence in Europe based on archaeological and genetic evidences, written by Leonardi, Gerbault, Thomas, and Burger. This blog highlighted the speculative combination of famine and diarrhea to explain the quick adaptive sweep of lactase persistence in European populations. I was surprised that lactase persistence was argued by Cordain to provide advantage against malaria in Africa and Southern Europe and rickets in Northern Europe.
I found many interesting ideas in the comments posted on this blog, often with entertaining exchanges. One reader linked another blog on evolution of lactose tolerance at darwinstudents.blogspot.com, that also gave some interesting speculations. I have to say that some of the discussions on lactose intolerance in Chinese population contain inaccurate facts.
I found many interesting ideas in the comments posted on this blog, often with entertaining exchanges. One reader linked another blog on evolution of lactose tolerance at darwinstudents.blogspot.com, that also gave some interesting speculations. I have to say that some of the discussions on lactose intolerance in Chinese population contain inaccurate facts.
Saturday, December 22, 2012
Read pairwise network data into igraph, then do permutation.
#continental US demo, permutation of pariwse interaction network, igraph usage
require(RCurl)
URL = "https://docs.google.com/spreadsheet/pub?key=0ArLJZixvTlU7dDI5c3dnTzdRX1dndzBORUk1UC0wYlE&single=true&gid=0&output=csv"
states = read.csv(textConnection(getURL(URL)), colClass = c("character", "character"))
g = graph.data.frame(states, directed=F)
g.degree = degree(g)
g.degree [g.degree == max(g.degree)] #TN and MO have 8 bordering states
g.shortestpath.m = shortest.paths(g)
str(g.shortestpath.m)
sorted.names = sort( rownames(g.shortestpath.m) )
gsm = g.shortestpath.m[, sorted.names]
gsm = gsm[sorted.names, ]
URL2 = "https://docs.google.com/spreadsheet/pub?key=0ArLJZixvTlU7dC1rTlRhTkw5UlY0eWZGT1JOVUVCWGc&single=true&gid=0&output=csv"
state.year.tb = read.csv(textConnection(getURL(URL2)), colClass = c("character", NA))
permute.pairs.wo.selfpairs = function( inpairs, ncycles=10, debug=1 ) {
if (ncycles >= 1 ) {
if(debug) {
print(paste('ncycles=', ncycles))
}
longids = c(as.character(inpairs[,1]), as.character(inpairs[,2]) )
longids = sample(longids)
len = length(inpairs[,1])
newpairs = data.frame( cbind( longids[1:len], longids[(len+1): (2*len)]) )
names(newpairs) = c('id1', 'id2')
newpairs$id1 = as.character( newpairs$id1)
newpairs$id2 = as.character( newpairs$id2)
newpairs$selfpairs = ifelse( newpairs$id1 == newpairs$id2, 1, 0 )
self.tb = newpairs[ newpairs$selfpairs==1, ]
nonself.tb = newpairs[newpairs$selfpairs==0, ]
if(debug) {
print(self.tb)
}
if( length(self.tb[,1])>0 ) {
if ( ncycles == 0) { return (c(NA,NA, NA) )
} else {
ncycles = ncycles - 1
selectedpairs = rbind(self.tb, nonself.tb[1: (length(self.tb[,1])*2), ] )
restpairs = nonself.tb[ (length(self.tb[,1])*2+1): length(nonself.tb[,1]), ]
return( rbind(restpairs, permute.pairs.wo.selfpairs(selectedpairs, ncycles)))
}
} else {
return (newpairs)
}
} else {
return( c(NA,NA,NA ))
}
}
#test
x = permute.pairs.wo.selfpairs(states)
x2 = permute.pairs.wo.selfpairs(states)
require(RCurl)
URL = "https://docs.google.com/spreadsheet/pub?key=0ArLJZixvTlU7dDI5c3dnTzdRX1dndzBORUk1UC0wYlE&single=true&gid=0&output=csv"
states = read.csv(textConnection(getURL(URL)), colClass = c("character", "character"))
g = graph.data.frame(states, directed=F)
g.degree = degree(g)
g.degree [g.degree == max(g.degree)] #TN and MO have 8 bordering states
g.shortestpath.m = shortest.paths(g)
str(g.shortestpath.m)
sorted.names = sort( rownames(g.shortestpath.m) )
gsm = g.shortestpath.m[, sorted.names]
gsm = gsm[sorted.names, ]
URL2 = "https://docs.google.com/spreadsheet/pub?key=0ArLJZixvTlU7dC1rTlRhTkw5UlY0eWZGT1JOVUVCWGc&single=true&gid=0&output=csv"
state.year.tb = read.csv(textConnection(getURL(URL2)), colClass = c("character", NA))
permute.pairs.wo.selfpairs = function( inpairs, ncycles=10, debug=1 ) {
if (ncycles >= 1 ) {
if(debug) {
print(paste('ncycles=', ncycles))
}
longids = c(as.character(inpairs[,1]), as.character(inpairs[,2]) )
longids = sample(longids)
len = length(inpairs[,1])
newpairs = data.frame( cbind( longids[1:len], longids[(len+1): (2*len)]) )
names(newpairs) = c('id1', 'id2')
newpairs$id1 = as.character( newpairs$id1)
newpairs$id2 = as.character( newpairs$id2)
newpairs$selfpairs = ifelse( newpairs$id1 == newpairs$id2, 1, 0 )
self.tb = newpairs[ newpairs$selfpairs==1, ]
nonself.tb = newpairs[newpairs$selfpairs==0, ]
if(debug) {
print(self.tb)
}
if( length(self.tb[,1])>0 ) {
if ( ncycles == 0) { return (c(NA,NA, NA) )
} else {
ncycles = ncycles - 1
selectedpairs = rbind(self.tb, nonself.tb[1: (length(self.tb[,1])*2), ] )
restpairs = nonself.tb[ (length(self.tb[,1])*2+1): length(nonself.tb[,1]), ]
return( rbind(restpairs, permute.pairs.wo.selfpairs(selectedpairs, ncycles)))
}
} else {
return (newpairs)
}
} else {
return( c(NA,NA,NA ))
}
}
#test
x = permute.pairs.wo.selfpairs(states)
x2 = permute.pairs.wo.selfpairs(states)
Thursday, December 20, 2012
The meaning of dN/dS (Ka/Ks) ratio in molecular evolution
Omega, Ka/Ks or dN/dS, is the ratio of nonsynonymous substitution per site to synonymous substitution per site. Ks (dS) can be viewed as a proxy of background mutation. The ratio of Ka/Ks, omega, indicate the extent of changes at the amino acid levels after normalized by silent mutational changes at DNA level, hence, is a proxy for positive selection pressure in coding genes. This definition assumes selection only at the protein level, not at the DNA or RNA level. It is clear that omega is only applicable in protein-coding genes.
Omega is designed to study divergence because its definition assumes fixed changes. When Omega is used in the population context, standing genetic variation would be treated as fixed ones, which can inflate the estimation of omega.
Kryazhimskiy and Plotkin (2008) showed that in population genetics context, dN/dS ratio is not a reliable index for selection pressure. When population scaled mutation rate is small enough, the population omega depends on both selection coefficient (s) and mutation rate (\mu). Based its figure 1, it can be seen that population omega is generally higher on divergent omega for negative selections. KP08 used a two-allele model to find an analytic solution, and then used simulation to study the behavior in a Markov chain model with 64 states. KP08 is mostly well written. I found Equation 3 was presented without explanations.
Omega is designed to study divergence because its definition assumes fixed changes. When Omega is used in the population context, standing genetic variation would be treated as fixed ones, which can inflate the estimation of omega.
Kryazhimskiy and Plotkin (2008) showed that in population genetics context, dN/dS ratio is not a reliable index for selection pressure. When population scaled mutation rate is small enough, the population omega depends on both selection coefficient (s) and mutation rate (\mu). Based its figure 1, it can be seen that population omega is generally higher on divergent omega for negative selections. KP08 used a two-allele model to find an analytic solution, and then used simulation to study the behavior in a Markov chain model with 64 states. KP08 is mostly well written. I found Equation 3 was presented without explanations.
Wednesday, December 19, 2012
Check and balance between adaptation and mutation in E. coli
Wielgoss et al, 2012 PNAS, Mutation rate dynamics in a bacterial population reflect tension between adaptation and genetic load. The authors sequenced one E. coli strain that has been evolved for 20 years in the Lenski lab. The author estimated that 20 years of E coli is about 50,000 generations. Twenty-two clones sampled at various time point during evolution were resequenced for this study. Previously, it was found that a frame-shift in mutT became 'established' in the population around 25K generations. After the fixation of this hypermutator mutT, mutation rate in two branches show reducted mutation rates (confirmed by fluctuation tests). Because the number of substitutions can be revealed by sequencing and the number of generations can be estimated, the mutations rates at different time points (hence, mutation dynamics) can be monitored (This is very nice thanks to the 20 years of in-vitro evolution experiment).
For reading, a few key concepts here are mutation rate, genetic load, hypermutators. Mutation rate are measured using synonymous substitutions. Genetic load is the normalized difference between maximal and mean fitness. The authors seem to use relative growth rate as compared to the ancestral clone to calculate the genetic load.
Specific terms include mutT, mutY, which are the specific genes under discussion.
I noticed the use of "established" instead of 'fixed' in the population, which indicates that mutT did not reach 100% of the population.
A question: How does the author establish the causal interaction between mutT and mutY? It seems that they begin with a candidate gene approach, and used association arguments. I was hoping for a regression based causal inference like those in quantitative genetics, but this kind of analysis may not be available here.
For reading, a few key concepts here are mutation rate, genetic load, hypermutators. Mutation rate are measured using synonymous substitutions. Genetic load is the normalized difference between maximal and mean fitness. The authors seem to use relative growth rate as compared to the ancestral clone to calculate the genetic load.
Specific terms include mutT, mutY, which are the specific genes under discussion.
I noticed the use of "established" instead of 'fixed' in the population, which indicates that mutT did not reach 100% of the population.
A question: How does the author establish the causal interaction between mutT and mutY? It seems that they begin with a candidate gene approach, and used association arguments. I was hoping for a regression based causal inference like those in quantitative genetics, but this kind of analysis may not be available here.
Tuesday, December 18, 2012
Dietary restriction in chronologicaled aged cells leads to extened replicative life span, Delaney et al, 2012 Exp Gen
Delaney et al, 2012 Experimental Gerontology
A paper from the Kaeberlein group ( http://www.ncbi.nlm.nih.gov/pubmed/23235143) showed the dietary restriction during chronologically aging can lead to extended replicative life span. This is a technically tricky experiment, because when chronologically aged cells start to lose viability, much more cells will be required to conduct the subsequent replicative life span assays. They author stated that they pick as much as 1000 cells to start the replicative lifespan assays.
One of my questions is a technique one: How did they carry out DR during chronological aging? It seems that starting cultures are SC with various glucose concentration. Five microliters of chronologically aged cultures were spread to YPD (2% glucose) plate for RLS assay, which is a standard practice. DR was achieved in 0.05% glucose in this study.
The authors double-stained yeast cells with 5nM SYTOX red and 17.5nM DiOC6, and monitor them by flow cytometry. DiOC6 is a mitochrondrial membrane potential dependent dye, and SYTOX Red is a used distinguish live and dead cells. SYTOX Red-negative cells are considered live cells.
This paper argues that both DR and buffering during chronological aging can delay the subsequent replicating aging process by maintiaing a low mitochondria membrane potential. Overall, the paper suggests that mitochondria is a link between post-mitotic and mitotic aging.
One question that I have is the role of mitochondria during chronological aging. The author conducted chronological aging in rich media. Although many cells are probably metabolically active during in the depleted media, many are also metabolically inactive. Is it possible that DR have differential effect on different cell sub-populations in different metabolic states?
For comparison, I often measure CLS in water, and many cells are metabolically much less active, if not domant in G0. DR cannot be studied in water-based chronological aging.
A paper from the Kaeberlein group ( http://www.ncbi.nlm.nih.gov/pubmed/23235143) showed the dietary restriction during chronologically aging can lead to extended replicative life span. This is a technically tricky experiment, because when chronologically aged cells start to lose viability, much more cells will be required to conduct the subsequent replicative life span assays. They author stated that they pick as much as 1000 cells to start the replicative lifespan assays.
One of my questions is a technique one: How did they carry out DR during chronological aging? It seems that starting cultures are SC with various glucose concentration. Five microliters of chronologically aged cultures were spread to YPD (2% glucose) plate for RLS assay, which is a standard practice. DR was achieved in 0.05% glucose in this study.
The authors double-stained yeast cells with 5nM SYTOX red and 17.5nM DiOC6, and monitor them by flow cytometry. DiOC6 is a mitochrondrial membrane potential dependent dye, and SYTOX Red is a used distinguish live and dead cells. SYTOX Red-negative cells are considered live cells.
This paper argues that both DR and buffering during chronological aging can delay the subsequent replicating aging process by maintiaing a low mitochondria membrane potential. Overall, the paper suggests that mitochondria is a link between post-mitotic and mitotic aging.
One question that I have is the role of mitochondria during chronological aging. The author conducted chronological aging in rich media. Although many cells are probably metabolically active during in the depleted media, many are also metabolically inactive. Is it possible that DR have differential effect on different cell sub-populations in different metabolic states?
For comparison, I often measure CLS in water, and many cells are metabolically much less active, if not domant in G0. DR cannot be studied in water-based chronological aging.
Monday, December 17, 2012
Generate pruned phylogeny using Biopython
Biopython Phylo
Its cookbook (http://biopython.org/wiki/Phylo_cookbook) offers some useful examples. Biopython is still tied to python 2.7.
The following code will read a Newick tree, generate a lookup dictionry for leaf-nodes, and then remove individual leaf node to generate a pruned tree.
from Bio import Phylo
#from Bio.Phylo import PhyloXML, NewickIO
def lookup_by_names(tree):
names = {}
for clade in tree.get_terminals():
if clade.name:
if clade.name in names:
raise ValueError("Duplicate key: %s" % clade.name)
names[clade.name] = clade
return names
EX_NEWICK = 'spec.nwk'
treeA = Phylo.read(EX_NEWICK, 'newick')
print(treeA)
names = lookup_by_names(treeA)
treeA.prune(names['Bcl']) #prune() takes an object from the same tree
treeA.count_terminals()
treeA.prune(names['Bha']) #prune() takes an object from the same tree
treeA.count_terminals()
treeA.prune(names['Bsu']) #prune() takes an object from the same tree
treeA.count_terminals()
treeA.prune(names['Bpu']) #prune() takes an object from the same tree
treeA.count_terminals()
Its cookbook (http://biopython.org/wiki/Phylo_cookbook) offers some useful examples. Biopython is still tied to python 2.7.
The following code will read a Newick tree, generate a lookup dictionry for leaf-nodes, and then remove individual leaf node to generate a pruned tree.
from Bio import Phylo
#from Bio.Phylo import PhyloXML, NewickIO
def lookup_by_names(tree):
names = {}
for clade in tree.get_terminals():
if clade.name:
if clade.name in names:
raise ValueError("Duplicate key: %s" % clade.name)
names[clade.name] = clade
return names
EX_NEWICK = 'spec.nwk'
treeA = Phylo.read(EX_NEWICK, 'newick')
print(treeA)
names = lookup_by_names(treeA)
treeA.prune(names['Bcl']) #prune() takes an object from the same tree
treeA.count_terminals()
treeA.prune(names['Bha']) #prune() takes an object from the same tree
treeA.count_terminals()
treeA.prune(names['Bsu']) #prune() takes an object from the same tree
treeA.count_terminals()
treeA.prune(names['Bpu']) #prune() takes an object from the same tree
treeA.count_terminals()
Saturday, December 15, 2012
DendroPy tutorial
Followed DendroPy's tutorial and I was able to prune a phylogenetic tree
in 50 minutes. DendroPy's example scripts can be found both on its website and in the examples/ folder in the downloaded source code package.
The interactive mode of Python is also very useful in this learning experience. It should be noted that DendroPy does not support Python3.0, which is a problem shared by most other bioinformatics related python projects, such as biopython.
import dendropy
dir(dendropy)
help(dendropy.Tree)
tree_str = "[&R] ((A, (B, (C, (D, E)))),(F, (G, H)));"
tree = dendropy.Tree.get_from_string(
tree_str,
"newick")
dir(tree)
print("Before:")
print(tree.as_string('newick'))
print(tree.as_ascii_plot())
tree.prune_taxa_with_labels(["A", "C", "G"])
print("After:")
print(tree.as_string('newick'))
print(tree.as_ascii_plot())
Later, I found that DendroPy cannot parse nodes with multiple labels, such as those output by codeml.
The interactive mode of Python is also very useful in this learning experience. It should be noted that DendroPy does not support Python3.0, which is a problem shared by most other bioinformatics related python projects, such as biopython.
import dendropy
dir(dendropy)
help(dendropy.Tree)
tree_str = "[&R] ((A, (B, (C, (D, E)))),(F, (G, H)));"
tree = dendropy.Tree.get_from_string(
tree_str,
"newick")
dir(tree)
print("Before:")
print(tree.as_string('newick'))
print(tree.as_ascii_plot())
tree.prune_taxa_with_labels(["A", "C", "G"])
print("After:")
print(tree.as_string('newick'))
print(tree.as_ascii_plot())
Later, I found that DendroPy cannot parse nodes with multiple labels, such as those output by codeml.
Monday, November 12, 2012
botseq_hong_variant_calling.mk
# Run freebayes on bam files
#
.DELETE_ON_ERROR:
ROOT_DIR?=/ndata/botlabseq/seqshare
include ${ROOT_DIR}/src/config.mk
include ${ROOT_DIR}/src/hong_analysis/botseq_hong_mapping.mk
REFERENCE_INDEX?=${ROOT_DIR}/reference_genomes/saccharomyces_cerevisiae_s288c_saccer3/saccharomyces_cerevisiae_s288c_saccer3.fasta
HONG_VCF_DIR=${ROOT_DIR}/hong_analysis/vcf
HONG_VCF_FILES=$(addprefix ${HONG_VCF_DIR}/, $(patsubst %.bam,%.vcf,$(notdir ${HONG_SORTEDBAM_FILES})))
HONG_VCF_GZ_FILES=$(addsuffix .gz, ${HONG_VCF_FILES})
HONG_VCF_INDEX_FILES=$(addsuffix .tbi, ${HONG_VCF_GZ_FILES})
.DEFAULT_GOAL:=hongvcall
hongvcall: ${HONG_VCF_INDEX_FILES}
${HONG_VCF_DIR}:
mkdir ${HONG_VCF_DIR}
.SECONDARY: ${HONG_VCF_FILES}
${HONG_VCF_FILES}: ${HONG_VCF_DIR}/%.vcf: ${HONG_MAPPED_BAM_DIR}/%.bam | ${HONG_VCF_DIR}
${FREEBAYES} --fasta ${REFERENCE_INDEX} --pooled $< > $@
${HONG_VCF_GZ_FILES}: ${HONG_VCF_DIR}/%.gz: ${HONG_VCF_DIR}/%
${BGZIP} $<
${HONG_VCF_INDEX_FILES}: ${HONG_VCF_DIR}/%.tbi: ${HONG_VCF_DIR}/%
${TABIX} $<
hongvctest:
@echo ${HONG_VCF_INDEX_FILES}
#
.DELETE_ON_ERROR:
ROOT_DIR?=/ndata/botlabseq/seqshare
include ${ROOT_DIR}/src/config.mk
include ${ROOT_DIR}/src/hong_analysis/botseq_hong_mapping.mk
REFERENCE_INDEX?=${ROOT_DIR}/reference_genomes/saccharomyces_cerevisiae_s288c_saccer3/saccharomyces_cerevisiae_s288c_saccer3.fasta
HONG_VCF_DIR=${ROOT_DIR}/hong_analysis/vcf
HONG_VCF_FILES=$(addprefix ${HONG_VCF_DIR}/, $(patsubst %.bam,%.vcf,$(notdir ${HONG_SORTEDBAM_FILES})))
HONG_VCF_GZ_FILES=$(addsuffix .gz, ${HONG_VCF_FILES})
HONG_VCF_INDEX_FILES=$(addsuffix .tbi, ${HONG_VCF_GZ_FILES})
.DEFAULT_GOAL:=hongvcall
hongvcall: ${HONG_VCF_INDEX_FILES}
${HONG_VCF_DIR}:
mkdir ${HONG_VCF_DIR}
.SECONDARY: ${HONG_VCF_FILES}
${HONG_VCF_FILES}: ${HONG_VCF_DIR}/%.vcf: ${HONG_MAPPED_BAM_DIR}/%.bam | ${HONG_VCF_DIR}
${FREEBAYES} --fasta ${REFERENCE_INDEX} --pooled $< > $@
${HONG_VCF_GZ_FILES}: ${HONG_VCF_DIR}/%.gz: ${HONG_VCF_DIR}/%
${BGZIP} $<
${HONG_VCF_INDEX_FILES}: ${HONG_VCF_DIR}/%.tbi: ${HONG_VCF_DIR}/%
${TABIX} $<
hongvctest:
@echo ${HONG_VCF_INDEX_FILES}
Subscribe to:
Posts (Atom)