use Entrez Direct to access NCBI database

I was reading Applied Bioinformatics 2014 lecture 3 and learned that one can use Entrez Direct to access NCBI database (Pubmed, nucleotide, protein sequence etc).

After installing Entrez Direct, I played around with it:

Commonly-used fields for PubMed queries include:

  [AFFL]  Affiliation       [FILT]  Filter              [MESH]  MeSH Terms
  [ALL]   All Fields        [JOUR]  Journal             [PTYP]  Publication Type
  [AUTH]  Author            [LANG]  Language            [WORD]  Text Word
  [FAUT]  Author - First    [MAJR]  MeSH Major Topic    [TITL]  Title
  [LAUT]  Author - Last     [SUBH]  MeSH Subheading     [TIAB]  Title/Abstract 

[PDAT] Date - Publication [UID] UID

Filters that limit search results to subsets of PubMed include:

  humans [MESH]                has abstract [FILT]
  pharmacokinetics [MESH]      historical article [FILT]
  chemically induced [SUBH]    loprovflybase [FILT]
  all child [FILT]             randomized controlled trial [FILT]
  english [FILT]               clinical trial, phase ii [PTYP]
  free full text [FILT]        review [PTYP]

Sequence databases are indexed with a different set of search fields, including:

  [ACCN]  Accession       [GENE]  Gene Name            [PROT]  Protein Name
  [ALL]   All Fields      [JOUR]  Journal              [SQID]  SeqID String
  [AUTH]  Author          [KYWD]  Keyword              [SLEN]  Sequence Length
  [GPRJ]  BioProject      [MLWT]  Molecular Weight     [SUBS]  Substance Name
  [ECNO]  EC/RN Number    [ORGN]  Organism             [WORD]  Text Word
  [FKEY]  Feature Key     [PACC]  Primary Accession    [TITL]  Title
  [FILT]  Filter          [PROP]  Properties           [UID]   UID

and a sample query in the protein database is:

  "alcohol dehydrogenase [PROT] NOT (bacteria [ORGN] OR fungi [ORGN])"

Please refer to the documents for more examples http://www.ncbi.nlm.nih.gov/books/NBK179288/

Install Inkscape on Mac

It has been a while since I wrote my last post. I am now in China and the Internet connection is very bad... I will start my postdoc in Dr.Role Verhaak's lab at MD Anderson Cancer Center. Dr. Verhaak's lab studies genomic alternations of brain tumor by analyzing whole exome sequencing, RNA-seq, whole genome sequencing, methylation and copy-number data. Yes, I am going to do a postdoc on computational biology. For sure, my computational skills would be strengthened in Dr.Verhaak's lab. Moreover, I will not give up my bench skills. I will use experiments to validate the functions of computational predictions.

I am writing a proposal and want to make a figure with Inkscape. It is a very good drawing software to deal with vector based figures, for bitmap based images, Gimp is the right one. On a linux machine (I have a ubuntu machine), inkscape can be installed by:

$sudo apt-get install inkscape

I only have a mac machine now at hand, so I have to install inkscape on my mac. Installation on mac is different from that on linux. Mac OS needs Xquartz to be installed first. I followed instructions here http://stackoverflow.com/questions/21049815/mavericks-trying-to-install-xquartz-x11-for-inkscape-image-not-found and here https://www.youtube.com/watch?v=7kvSc_PYokM

And now it is ready to use!  I am excited to start my new job at MD Anderson and I should have an exciting 2015. I will update more frequently on this blog as I learn new stuff in my new lab.

csvkit to manipulate csv at command line, Rio to interact with R at command line

I came across csvkit long time ago, but I just began to use it and found it is very handy. It is a python utility to manipulate csv files at command line. you can install it by pip:

pip install csvkit

I played around with the iris data.
csvlook gives you a database view of the data by adding grids.
$ cat iris.csv | head | csvlook
|---------------+-------------+--------------+-------------+--------------|
|  sepal_length | sepal_width | petal_length | petal_width | species      |
|---------------+-------------+--------------+-------------+--------------|
|  5.1          | 3.5         | 1.4          | 0.2         | Iris-setosa  |
|  4.9          | 3.0         | 1.4          | 0.2         | Iris-setosa  |
|  4.7          | 3.2         | 1.3          | 0.2         | Iris-setosa  |
|  4.6          | 3.1         | 1.5          | 0.2         | Iris-setosa  |
|  5.0          | 3.6         | 1.4          | 0.2         | Iris-setosa  |
|  5.4          | 3.9         | 1.7          | 0.4         | Iris-setosa  |
|  4.6          | 3.4         | 1.4          | 0.3         | Iris-setosa  |
|  5.0          | 3.4         | 1.5          | 0.2         | Iris-setosa  |
|  4.4          | 2.9         | 1.4          | 0.2         | Iris-setosa  |
|---------------+-------------+--------------+-------------+--------------|

csvstat can output summary statistics for specific columns with -c flag

-c can accept both the column names and the column numbers as arguments

$ cat iris.csv | csvstat -c sepal_length,sepal_width

  1. sepal_length

<type 'float'>

Nulls: False

Min: 4.3

Max: 7.9

Sum: 876.5

Mean: 5.84333333333

Median: 5.8

Standard Deviation: 0.825301291785

Unique values: 35

5 most frequent values:

5.0: 10

6.3: 9

5.1: 9

6.7: 8

5.7: 8

  2. sepal_width

<type 'float'>

Nulls: False

Min: 2.0

Max: 4.4

Sum: 458.6

Mean: 3.05733333333

Median: 3.0

Standard Deviation: 0.434410967735

Unique values: 23

5 most frequent values:

3.0: 26

2.8: 14

3.2: 13

3.4: 12

3.1: 11

Row count: 150

you can also output only the max, min, mean by adding --max, --min etc

$ cat iris.csv | csvstat -c sepal_length,sepal_width --max

  1. sepal_length: 7.9

  2. sepal_width: 4.4

csvcut can cut out specific columns and reorder the column while the unix cut command can

not reorder the column, one needs to use awk to achieve that.

$ cat iris.csv | csvcut -c5,4,2 | csvlook | head

|------------------+-------------+--------------|

|  species         | petal_width | sepal_width  |

|------------------+-------------+--------------|

|  Iris-setosa     | 0.2         | 3.5          |

|  Iris-setosa     | 0.2         | 3.0          |

|  Iris-setosa     | 0.2         | 3.2          |

|  Iris-setosa     | 0.2         | 3.1          |

|  Iris-setosa     | 0.2         | 3.6          |

|  Iris-setosa     | 0.4         | 3.9          |

|  Iris-setosa     | 0.3         | 3.4          |

csvgrep can grep based on a specific column:

$ cat iris.csv | csvgrep -c species -m setosa | wc -l

51

other commands like csvsort, csvjoin, csvsql are also useful, please refer to the link above to see more examples.

I want to plot some nice figure at command line, but I am too lazy to fire Rstudio. by using Rio, you can leverage the power of R at the command line. you can then visualize it by using the display command from imagemagic

-g to load the ggplot2 package 

-e to execute the command

update on 11/19/2014, Rio has changed the flag to load dplry and tidyr
https://github.com/jeroenjanssens/data-science-at-the-command-line/commit/66fc1fda604f2d420de9490934259093a76fb105

boxplot of sepal_length for different species :

histogram for sepal_length for all species:

 you can change the color and fill of the bars:

you can also plot histogram for different species with a facet:

scatter plot of sepal_width for different species:

command line is awesome! I love it! It brings me more power!

Test significance of overlapping genes from 2 gene sets

The question is very straightforward:
I have three microarray or RNA-seq data sets (control, treatment1, treatment2, each should have several replicates, otherwise you can not get a pvalue for differential expression!) , I then get two sets of differentially up-regulated genes ( treatment1 vs control, treatment2 vs control).  There are some overlapping genes which are both up-regulated by treatment1 and treatment2. How significant the overlapping is?

The easiest way is to use an online program http://nemates.org/MA/progs/overlap_stats.html
The program calculates the p value by using Hypergeometric_distribution.

see below for a python and an R script for this problem.

read  posts here http://stats.stackexchange.com/questions/16247/calculating-the-probability-of-gene-list-overlap-between-an-rna-seq-and-a-chip-c
http://blog.nextgenetics.net/?e=94
https://www.biostars.org/p/90662/

list1=3000, list2=400, Total gene number (population)=15000,
overlapping between list1 and list2 =100
From the stackexchange link:
"The p-value you want is the probability of getting 100 or more white balls in a sample of size 400 from an urn with 3000 white balls and 12000 black balls. Here are four ways to calculate it.

sum(dhyper(100:400, 3000, 12000, 400))
1 - sum(dhyper(0:99, 3000, 12000, 400))
phyper(99, 3000, 12000, 400, lower.tail=FALSE)
1-phyper(99, 3000, 12000, 400)

These give 0.0078.

dhyper(x, m, n, k) gives the probability of drawing exactly x. In the first line, we sum up the probabilities for 100 – 400; in the second line, we take 1 minus the sum of the probabilities of 0 – 99.

phyper(x, m, n, k) gives the probability of getting x or fewer, so phyper(x, m, n, k) is the same as sum(dhyper(0:x, m, n, k)).

The lower.tail=FALSE is a bit confusing. phyper(x, m, n, k, lower.tail=FALSE) is the same as 1-phyper(x, m, n, k), and so is the probability of x+1 or more "

If I run with the python script, I got the same p-value:

tommy@tommy-ThinkPad-T420 ~/scripts_general_use
$ ./hyper_geometric.py 15000 3000 400 100

total number in population: 15000
total number with condition in population: 3000
number in subset: 400
number with condition in subset: 100

p-value <= 100: 0.996049535673

p-value >= 100: 0.00784035207977

The bottom line is that we need more statistics!!

convert bam file to bigwig file and visualize in UCSC genome browser in a Box (GBiB)

I just came across UCSC genome browser in a box Gbib, and I want to try it out.

what is Gbib?
from the link above:

"Genome Browser in a Box (GBiB) is a "virtual machine" of the entire UCSC Genome Browser website that is designed to run on most PCs (Windows, Mac OSX or Linux). GBiB allows you to access much of the UCSC Genome Browser's functionality from the comfort of your own computer. It is particularly directed at individuals who want to use the Genome Browser's functionality to view protected data. With the anticipation that the majority of protected data use will focus on the human genomes (primarily GRCh37/hg19, and transitioning to GRCh38/hg38), the current version of GBiB has been optimized for use with the hg19 assembly. Many other recent genome assemblies can also be viewed, but without mirroring of additional data, they may be slow."

Download the compressed file from here https://genome-store.ucsc.edu/ and follow the installation instructions. Open virtual box and start the virtual machine and then direct your web browser to 127.0.0.1:1234, you should see the UCSC genome browser page.

It looks the same as the UCSC genome browser page https://genome.ucsc.edu/. You can then go to the browser and do what you usually do in the genome browser. when I searched a gene, and errors came out:

• nextRow failed
• mySQL error 1194: Table 'rmsk' is marked as crashed and should be repaired (profile=db, host=localhost, db=hg19)

I asked the UCSC mail list, and I typed the command below in the console:
 sudo mysqlcheck --all-databases --auto-repair --fast
it works like a charm now!

I have some bam files and I want to visualize them in the UCSC genome browser. I would first convert them to bigwig http://genome.ucsc.edu/goldenpath/help/bigWig.html and then upload to some open http https or ftp locations.  see my old post:
hosting bigwig by Dropbox for UCSC

However, with GBiB one only needs to share a local folder containing the bigwig files with virtual box.

"Loading Local Big Data Tracks

Your computer can share directories with GBiB so that you can load big files without the need to upload them to a web server. The big file formats are compressed, indexed, binary files, and include bigBed,bigWig, BAM, and VCF files. Normally, one would have to place these types files onto a publicly accessible web server to upload them as a custom track. However, GBiB acts as its own web server, meaning that you can share local files with GBiB for easy uploading as a custom track."

I will start from converting bam to bigwig.
several posts you may want to have a look:
https://www.biostars.org/p/64495/#64680
https://www.biostars.org/p/2699/
http://onetipperday.blogspot.com/2012/07/three-ways-to-convert-bambed-file-to.html

I used genomeCoverageBed from bedtools to convert bam to bedgraph first, and then used UCSC utilities bedClip and bedGraphToBigWig http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/ to do the conversion.

bedtools requires the bam files be sorted. so samtools sort was used to sort the bam file first.
I used a script from Tao Liu (the author of famous MACS ChIP-seq calling software) to do the conversion after I got the bedgraph files. I then put the bigwig files in a folder and shared it with virtual box. Follow the instructions http://genome.ucsc.edu/goldenPath/help/gbib#LocalTracks. I successfully visualized a bunch of bigwig files!

speed up grep

I had a list of gene names in txt file. There are around 500 genes with one gene name in one line, and I want to filter the gtf file from ensemble Homo_sapines.GRCh37.74gtf.gz

the gtf file contains 2244857 lines. I used grep to do it, but it takes very long (~1 hour).

what I used:

zcat Homo_sapines.GRCh37.74gtf.gz | grep -f gene_names.txt -w > my_genes.gtf

I searched on line, and found several posts in stackoverflow to speed up grep:
http://stackoverflow.com/questions/14602963/faster-grep-function-for-big-27gb-files
http://stackoverflow.com/questions/13913014/grepping-a-huge-file-80gb-any-way-to-speed-it-up
http://stackoverflow.com/questions/9066609/fastest-possible-grep

options to speed up:

1) Prefix your grep command with LC_ALL=C to use the C locale instead of UTF-8.

2) Use fgrep because you're searching for a fixed string, not a regular expression

I then used:

zcat Homo_sapines.GRCh37.74gtf.gz | LC_ALL=C fgrep -f gene_names.txt -w > my_genes.gtf

It runs much faster!

everyone will write a blog post on dplyr

dply is a nice R package to manipulate big data. There are several good tutorials I came across:

http://www.dataschool.io/dplyr-tutorial-for-faster-data-manipulation-in-r/
http://stat545-ubc.github.io/block009_dplyr-intro.html
http://stat545-ubc.github.io/bit001_dplyr-cheatsheet.html
http://gettinggeneticsdone.blogspot.com/2014/08/do-your-data-janitor-work-like-boss.html

I followed the example in the early edition of Bioinformatics Data Skills on the dplyr part and put a gist below. The example used dplyr to manipulate and summarize the human annotation gff file downloaded here http://useast.ensembl.org/info/data/ftp/index.html.

I downloaded the ensemble 74 build which is different from the 75 build in the example of the book.
I believe the gencode v19 annotation file is based on ensemble 74.

I also included comments in the R code describing how to do the same job with linux command line.

make a soft link

I got to know the ln linux command when I was reading The linux command line. I did not realize why I need this command to make links to the same data or file.

One situation I faced is that I have a lot of bam files  with long names in one folder, and when I execute a program in another folder, I have to pass the program a full path to the bam files.

while you think the TAB auto-completion helps most of the time,  it is still too much to specify the full path + long names (ugly). So I want to make a copy of the bam files in the current working director to save me from typing the full path. This is where ln command comes to rescue.

The other scenario is when you have multiple executable files in one folder ~/myprogram/bin, you want to execute the program from anywhere, you could link it to /usr/local/bin, or you can just copy the executable to /usr/local/bin or you can add ~/myprogram/bin to $PATH. see a link here http://unix.stackexchange.com/questions/116508/adding-to-path-vs-linking-from-bin

Let me demonstrate the usage of ln 

look at the man page

man ln

# I have two files in the playground directory and make two directories 

mkdir dir1 dir2
ls

output: tommy@tommy-ThinkPad-T420[playground] ls                              [ 2:32PM]
dir1  dir2  genes_h.txt  genes.txt

I then move the two txt files into dir1 and redirect to dir2 and make soft links of files in dir1 to dir2

mv *txt dir1
cd dir2
ln -s ../dir1/*txt .
ls

if we go to dir1 and make a softlink from there

cd ../dir1
ln -s ./*txt  ../dir2

The links will not work (too many levels of symbolic levels error message if you want to use the link ), instead you need to specify the full path

ln -s $PWD/*txt  ../dir2

see links here http://www.thegeekstuff.com/2010/10/linux-ln-command-examples/
and http://www.thegeekstuff.com/2010/10/linux-ln-command-examples/

converting gene ids using bioconductor with biomaRt and annotation packages

I had a post using  mygene to convert gene ids. Bioconductor can do the same job.
I put a gist on github.

For more examples see posts from Dave Tang:
http://davetang.org/muse/2013/12/16/bioconductor-annotation-packages/
http://davetang.org/muse/2013/05/23/using-the-bioconductor-annotation-packages/
http://davetang.org/muse/2013/11/25/thoughts-converting-gene-identifiers/

mapping gene ids with mygene

Mapping gene ids is one of the routine jobs for bioinformatics. I was aware of several ways to do it including Biomart.

Update on 10/30/14, a mygene bioconductor package is online http://bioconductor.org/packages/release/bioc/html/mygene.html

Recently I got to know mygene, a python wrapper for the mygene.info services to map gene ids.
I found it very handy to convert gene ids.  see a gist below.

To use it,  cat input.txt | python geneSymbol2Entrez.py > output.txt
or python geneSymbol2Entrez.py input.txt > output.txt  where input.txt contains one gene name in each line. pretty neat!