Why don't you see all GO terms in GSEA-MSigDB?

 This is why...

This procedure has been modified from that described previously for MSigDB v5.2. First, for each GO term we got the corresponding human genes from the gene2go file. Next, we have applied the path rule. Gene products are associated with the most specific GO terms possible. All parent terms up to the root automatically apply to the gene product. Thus, the parent GO term gene sets should include all genes associated with the children GO terms. Then we removed sets with fewer than 5 or more than 2,000 Gene IDs. Finally, we resolved redundancies as follows. We computed Jaccard coefficients for each pair of sets, and marked a pair as highly similar if its Jaccard coefficient was greater than 0.85. We then clustered highly similar sets into "chunks" using the hclust function from the R stats package according to their GO terms and applied two rounds of filtering for every "chunk". First, we kept the largest set in the "chunk" and discarded the smaller sets. This left "chunks" of highly similar sets of identical sizes, which we further pruned by preferentially keeping the more general set (i.e., the set closest to the root of the GO ontology tree).

Ref: https://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/MSigDB_v5.2_Release_Notes

A bug related to R factor

Note a bug in my code today. Sometimes you need to put a certain level (e.g. healthy control) in the first position for your covariance. 

Here is my old code:

dds[[variable]]=factor(dds[[variable]])

levels(dds[[variable]])= union(variable_REF, levels(dds[[variable]])))

Note that this can cause problem. For example, you have two levels: HC and AD in your diagnosis. By default, setting the covariance to factor() will make the levels sorting in alphabetic order (e.g. AD then HC). In the 2nd line, if you force to change the levels to c("HC", "AD"), you are actually doing something like c(AD = "HC", HC = "AD"). This is not what you want. 

What you wanted is actually this:

dds[[variable]]=factor(dds[[variable]], levels= union(variable_REF, unique(dds[[variable]]))) 

* Too bad that blogger does not adapt to markdown language yet. 

compression level in samtools and gzip

samtools fastq can convert bam to fastq format, e.g. samtools fastq input.bam -o output.fastq

The output file will be automatically compressed if the file names have a .gz or .bgzf extension, e.g. 

samtools fastq input.bam -o output1.fastq.gz

Alternatively, you can also pipe the stdout to compressor explicitly, e.g.

samtools fastq input.bam | gzip > output2.fastq.gz

Interestingly, I noticed that output2.fastq.gz is significantly smaller than output1.fastq.gz, even though the uncompressed file content is the same. 

Actually, this is because of the different default compression ratios used in samtools and gzip. 

In samtools fastq, its default compression level is 1 (out of [0..9]) while gzip's default compression level is 6 (out of [1..9]). 

Conda - A must-have for bioinformatician

As a bioinformatician, when you are given a new system (Mac, HPC, or any Linux environment), the first thing to do is probably to configure your work environment, install required tools, etc. 

Here is what I did for the new Mac:

1. This step is optional. Apple changed to default shell to zsh since the last OS. If you like bash, you can change the default zsh shell from bash:

chsh -s /bin/bash

2. Install Miniconda. Miniconda is a free minimal installer for conda, incl. only the basic ones (e.g. conda, python, zlib, etc). If you are under Linux, download the one for Linux. 

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh

bash Miniconda3-latest-MacOSX-x86_64.sh

3. Install commonly used tools with bioconda:

# update bash (the one coming with MacOS is too old)

conda install -c conda-forge bash

# samtools etc.

conda install -c bioconda samtools plink plink2

# UCSC Kent Utilities etc.

conda install -c bioconda ucsc-bedgraphtobigwig ucsc-wigtobigwig ucsc-liftover ucsc-bigwigsummary

# GNU utils etc.

conda install -c conda-forge coreutils gawk wget gzip datamash parallel make readline sed

# R and R package etc.

conda install -c r r-base

conda install -c r r-essentials

conda install -c r r-tidyverse

# Google cloud gsutil

conda install -c conda-forge gsutil 

Life is much easier after conda!

Explaination of "samtools flagstat" output

==> CRR119891.sam.flagstat <==

192013910 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 secondary

38830640 + 0 supplementary

0 + 0 duplicates

191681231 + 0 mapped (99.83% : N/A)

153183270 + 0 paired in sequencing

76591635 + 0 read1

76591635 + 0 read2

118230182 + 0 properly paired (77.18% : N/A)

152555312 + 0 with itself and mate mapped

295279 + 0 singletons (0.19% : N/A)

10160332 + 0 with mate mapped to a different chr

5532624 + 0 with mate mapped to a different chr (mapQ>=5)

Here are the notes from our careful check:

1. What's secondary, supplementary, and singleton? I was using "bwa mem" with default parameters. In that case, it only reports one alignment if the read can be mapped to multiple locations with equal quality. See discussion here: https://www.biostars.org/p/304614/#304620. That's why you don't see the secondary alignment in above example. Supplementary alignments are those alignments that are part of a chimeric alignment. For example, in a circRNA case, one read is broke into two pieces, mapped to different locations in back-splicing order. Then one of the two fragment will be supplementary (the other is not). Singleton is a mapped read whose mate is unmapped. So its flag has 0x1 and 0x8 but 0x4 is unset. Note that a singleton's mate (which is unmapped) can have a flag of 0110.  

2. What's the total number of reads? The number in the output is NOT for reads, but rather for alignment in the SAM files. The samtools flagstat only check the FLAG, not the read ID. So, in the above example, the total number of reads (R1+R2) should be 192013910 - 38830640. (If there is secondary or duplicate, it should also be excluded). It's also indicated by "paired in sequencing", 153183270 (R1 +R2). This can also be calculated by checking the read ID column e.g. cat CRR119891.sam | cut -f1 | sort -u | wc -l, then times 2 for paired-end reads.

3. How many reads are unmapped? The total alignment is sum of "mapped" + "singletons" + unmapped. So, the number of unmapped reads should be calcualted as 192013910-191681231-295279 = 37400 in above case. The unmapped reads can also be counted by checking the RNAME column in the SAM file (e.g. cat CRR119891.sam | awk '$3=="*"' | wc -l). 

4. Mappability? The above mapped % is based on "mapped" / "total" alignment, not based on reads. If for reads, it should be calculated based on above two numbers. 

Be careful of "sort -k1,1 -k2,2n -u"

When you attempted to sort and extract the unique genomic regions using "sort -k1,1 -k2,2n -u", you might make a mistake by missing the region with the same chr and start, but different end position.

The right way should be  "sort -k1,1 -k2,2n -k3,3n -u" or  "sort -k1,1 -k2,2n | sort -u"

try to be consistent

Worthy to note the lesson I learned today: try to be consistent when sorting with different chromosome order.

I got the same error message as the post: https://github.com/arq5x/bedtools/issues/109

In short, the chromosome order in my sorted BAM file header is like this chr10, chr11, chr12... chr19, chr1, chr20, ... chr22, chr2, chr3 ....chrM, chrX, chrY. The order is determined by the aligner which generates @SQ headers based on how the reference sequences are arranged in the FASTA file from which you built its aligner indices. That means, the order in @SQ is same as the order in your genome.fa.fai (index) file.

Educational papers in Bioinformatics from Nature Methods and Nature Biotechnology

Nature Methods and Nature Biotechnology have this super-great monthly column in bioinformatics for years. For example, Nature Methods has both the "Points of Significance" covering topics in statistics and "Points of View" in data visualization. And Nature Biotechnology's Primer column covers a collection of topics in Bioinformatics, such as "What's a hidden Markov model?". 

Today I've compiled all the articles from above three columns together and made a PDF for you. Here is the link. 

My 15 practical tips for a bioinformatician

Tips below are based on the lessons I learnt from making mistakes during my years of research. It's purely personal opinion. Order doesn't mean anything. If you think I should include something else, please comment below.

1. Always set a seed number when you run tools with random option, e.g. bedtools shuffle, random etc.; You (or your boss in a day) want your work reproducible. 
2. Set your own temporary folder (via --tmp, $TMPDIR etc., depending on your program). By default, many tools, e.g. sort, use the system /tmp as temporary folder, which may have limited quote that is not enough for your big NGS data. 
3. Always use a valid name for your variables, column and rows of data frame. Otherwise, it can bring up unexpected problem, e.g. a '-' in the name will be transferred to '.' in R unless you specify check.names=F.
4. Always make a README file for the folder of your data; For why and how, read this: http://stackoverflow.com/questions/2304863/how-to-write-a-good-readme
5. Always comment your code properly, for yourself and for others, as you very likely will read your ugly code again.
6. Always backup your code timely, using github, svn, Time Machine, or simply copy/paste whatever.
7. Always clean up the intermediate or unnecessary data, as you can easily shock your boss and yourself by generating so much data (and perhaps most of them are useless). 
8. Don't save into *.sam if you can use *.bam. Always zip your fastq (and other large plain files) as much as you. This applies to other file format if you can use the compressed one. As you cannot imagine how much data (aka "digital garbage") you will generate soon.
9. Using parallel as much as you can, e.g. using "LC_ALL=C sort --parallel=24 --buffer-size=5G" for sorting (https://www.biostars.org/p/66927/), as multi-core CPU/GPU is common nowaday.
10. When a project is completed, remember to clean up your project folder, incl. removing the unnecessary code/data/intermediate files, and burn a CD or host it somewhere in cloud for the project. You never know when you, your boss or your collaborators will need the data again;
11. Make your code sustainable as possible as you can. Remember the 3 major features of OOP: Inheritance, Encapsulation, Polymorphism. (URL)
12. When you learn some tips from others by Google search, remember to list the URL for future reference and also for acknowledging others' credit. This applies to this post, of course :)
13. Keep learning, otherwise you will be out soon. Just like the rapid development of NGS techniques, computing skills are also evolving quickly. Always catch up with the new skills/tools/techniques.
14. When you learn some tips from others, remember to share something you learned to the community as well, as that's how the community grows healthily. 
15. Last but not least, stand up and move around after sitting for 1-2 hours. This is especially important for us bioinformaticians who usually sit in front of computer for hours. Only good health can last your career long. More reading: https://www.washingtonpost.com/news/wonk/wp/2015/06/02/medical-researchers-have-figured-out-how-much-time-is-okay-to-spend-sitting-each-day/

PLINK2 vs. SNAP for linkage disequilibrium (LD) calculation

Among other ways to calculate SNPs in linkage disequilibrium (LD), two methods have been used in many literatures: PLINK2 and SNAP

With PLINK2, SNPs in LD can be calculated using parameters like:
PLINK2 -r2 dprime --ld-window-kb 1000 --ld-window 10 --ld-window-r2 0.8
SNAP is a web-based service with pre-calculated data: https://www.broadinstitute.org/mpg/snap. You can set the population, distance limit, r2 threshold.

Note that PLINK2 has one more control: --ld-window

In its manual (https://www.cog-genomics.org/plink2/ld), it says "By default, when a limited window report is requested, every pair of variants with at least (10-1) variants between them, or more than 1000 kilobases apart".

That means PLINK2 limits SNPs with two kinds of radius: both the genomic distance (--ld-window-kb) and the number of SNPs apart (--ld-window), while SNAP only has the distance limit. Both of them have the r2 threshold.

How to extract the gap region in human genome?

Just notice that I should avoid the gap region, esp. when we generate a random background as your null distribution using tools such as bedtools shuffle.

Short answer: go below UCSC Table Browser link and choose to save as a bed file
 http://genome.ucsc.edu/cgi-bin/hgTables?clade=mammal&org=Human&db=hg19&hgta_group=allTables&hgta_track=hg19&hgta_table=gap&hgta_regionType=genome&hgta_outputType=primaryTable

As below table shown, 8.28% of hg19 assembly are simply gap.

Gap (gap) Summary Statistics

	item count 457 item bases 239,845,127 (8.28%) item total 239,845,127 (8.28%) smallest item 47 average item 524,825 biggest item 30,000,000

altColor in UCSC track setting

Many track types allow setting a color range that varies from color to altColor. For instance the CpG Island tracks use the altColor setting to display the weaker islands, while the stronger ones are rendered in color. If altColor is not specified, the system will use a color halfway between that specified in the color tag and white instead.

Be aware that wiggles with negative values are drawn in altColor not color as positive values are.

Hierarchical structure of UCSC Genome Browser track hub

UCSC Genome Browser tracks can be organized into groups by using the container multiWig, compositeTrack on, and superTrack on lines. Supertracks can contain composite tracks and container multiWigs, but not vice versa. With supertracks, composite tracks, and container multiWigs, children will inherit the settings from their parents, but can override their parent settings within their own stanzas.

Here is their hierarchical relationship:

superTrack on

|==== child tracks (bam, bigBed, bigWig or vcfTabix)

|==== container multiWig

     |---- child track (bigWig)

|==== compositeTrack on

     |==== view

         |---- child track with subGroups setting (bam, bigBed, bigWig or vcfTabix, but not mix)

         |---- child track with subGroups setting (bam, bigBed, bigWig or vcfTabix, but not mix) 

superTrack also allow to mix other supported types of hub tracks: bam, bigBed, bigWig or vcfTabix.

Simple script to generate random size-matched background regions

Here is a simple script I wrote for generating a random background regions with matched size distribution, using bedtools random and a input.bed as background. Not fancy, but works!

# ===============================================================
# Script to generate a random size-matched background regions 
# Author: Xianjun
# Date: Jun 3, 2014
# Usage:
# toGenerateRandomRegions.sh input.bed > output.random.bed
# or
# cat input.bed | toGenerateRandomRegions.sh -
# ===============================================================

bedfile=$1
# other possible options can be the genome (e.g. hg19, mm9 etc.) or an inclusion or exclusion regions (e.g. exons)

# download hg.genome
mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e "select chrom, size from hg19.chromInfo"  > hg19.genome

cat $bedfile | while read chr start end rest
do
    let l=$end-$start;
    bedtools random -g hg19.genome -n 1 -l $l
done

a bigWigSummary bug

Write down a bigwigsummary bug I found today. It's found when I attempted to get the max value in a region using -type=max and dataPoints=1:

$ bigWigSummary test.bw -type=max chr12 54070173 54072173 1

13.3672

$ bigWigSummary test.bw -type=max chr12 54070173 54072173 10

0.944904 1.02475 0.568405 0.741671 1.43119 1.08896 0.705965 0.542034 0.380971 0.591934

As you see, if I use dataPoints=1, the max value is 13.3672 and when I use dataPoints=10 the max value is 1.43119. So there must be something wrong, since the max value should not change no matter how many data points we check. Visualizing the bigwig in UCSC Genome Browser shows that 1.43119 is correct for this case. Interestingly, 13.3672 is indeed a peak summit, but not for this region, rather a region upstream. I don’t know why bigWigSummary take the summit from region outside. This only happened when dataPoints=1. 

I put the data below. You can download to test:

https://www.dropbox.com/sh/q963gxj9e1kj8ly/AAAG-LqCT_WQi6dgMm0lhoDIa?dl=0

People reported similar error for bigwigSummary, for example

bigWigSummary outputs different values as bigWigAverageOverBed: UCSC team explained "The reason for this is that the summary levels have some rounding error and some border conditions when extracting data over relatively small regions." They suggested to use bigWigAverageOverBed if you want the highest level of accuracy. But bigWigAverageOverBed won't output the max and also it's nothing with high or low accuracy, but rather a bug. 

Tao Liu also reported another bug for bigwig when it's converted from wig in compressed manner (by default), and suggested to fix it by using -unc when converting wig to bigwig. I've tried to use -unc when converting from bedGraph to bigwig, the bug is still there. 

Still looking for workout, and also report to UCSC:

https://groups.google.com/a/soe.ucsc.edu/forum/#!topic/genome/pWjcov-xQyQ

Update: One workout I found is, to use intersectBed and groupBy in bedtools on bedGraph file (rather than bigwig). Here is pseudocode:

intersectBed -a regions.bed -b signal.bedGraph -wo -sorted | groupBy -g 1,2,3 -c 7 -o max > regions.maxSignal.bed

Update2:  use the -minMax option from the latest version of bigWigAverageOverBed (from v304). See reply from UCSC group:

Thank you for contacting us. One of our engineers was able to reproduce this and says the bigWigSummary program uses the bigWig summary levels, which are lower resolution versions of the data, so when you ask for a single value in a particular range, the range that is used may include bases that are before and after the range.
If you want only the values exactly within the range, you can use (the latest version from v304) bigWigAverageOverBed like so:
echo "chr12 54070173 54072173 one" | bigWigAverageOverBed myTest.bw stdin stdout -minMax | cut -f 8

Trimmed mean and median in AWK

This can be easily done in R, but sometime you want to get it in scripting language like awk or perl in order to process the big data line by line.

Here is the code snip:

# for median
# Note: thanks to anonymous reply below, put c, j in arguments to define that they are local variables.

function median(v, c, j) 

{ 

    c=asort(v,j); 

    if (c % 2) return j[(c+1)/2]; 

    else return (j[c/2+1]+j[c/2])/2; 

}

# for trimmed mean (where p is the percentage of data to be trimmed)
function trimmedMean(v, p, c, j) 

{ 

    c=asort(v,j); 

    a=int(c*p);

    for(i=a+1;i<=(c-a);i++) s+=j[i];

    return s/(c-2*a); 

}

To use it, for example if we want to generate a merged bigwig track for a list of samples, we can take median value of all samples at each genomic position, here is it:

unionBedGraphs -i `ls *normalized.bedGraph` | awk 'function median(v) {c=asort(v,j); if (c % 2) return j[(c+1)/2]; else return (j[c/2+1]+j[c/2])/2.0; } {OFS="\t"; n=1; for(i=4;i<=NF;i++) S[n++]=$i; print $1,$2,$3, median(S)}' > AllSamples.median.normalized.bedGraph
bedGraphToBigWig AllSamples.median.normalized.bedGraph ChromInfo.txt AllSamples.median.normalized.bw

Confusing dbSNP ID

There have been enough confusion in ID or naming system in the popular biological databases. I mentioned one such example in an earlier post, where two genes (on totally different chromosomes) use the same symbol and the same gene has different Ensembl Gene ID. Other confusions include the Entrez ID vs. gene symbol vs. Ensembl ID. People who do GO analysis using tools like DAVID must have painful experience on this (btw, I don't understand why most GO analysis tools convert to Entrez ID at the end).

And now, I got another example of confusing ID system in biology or bioinformatics.

You expect each SNP should have their unique ID in the dbSNP database, right? But they are not. Look at this example:

chr1    13837   13838   -      rs7164031, rs79531918

chr1    13837   13838   +     rs200683566, rs28391190, rs28428499, rs71252448, rs79817774

At the same location, multiple SNPs ID are assigned, for both strands. 

Here is what I got from NCBI User service for the explanation:

This is expected for many reasons:

- short probes could have multiple mapping locations on the genome

- certain genes could have duplicate/pseudogene/paralogs

- variations found in repeat regions would be difficult to map to a unique 

location

An rsID represents a cluster of reported variations submitted to dbSNP.

In ideal situation, they can be mapped to a unique location in the genome. 

There are cases where such unique mapping is NOT attainable.

But I am still not clear of how a cluster of variants are defined. [UPDATE to add]

p.s. code snip to merge SNPs with the same location into single line:
grep single snp137.bed | sort -k1,1 -k2,2n -k6,6 | bedtools groupby -g 1,2,3,6 -c 4 -o collapse

about UCSC Genome Browser track

Two points:

1. The custom track data may be compressed by any of the following programs: gzip (.gz), compress (.Z), or bzip2 (.bz2). But not for bigwig and bam.

2. In a track hub Db configuration file, up to 9 subgroup types can be defined for a composite, such as:

subGroup1 <gTag1> <gTitle1> <mTag1a=mTitle1a> [mTag1b=mTitle1b…]
subGroup2 <gTag2> <gTitle2> <mTag2a=mTitle2a> [mTag2b= mTitle2b…]
...
subGroup9 <gTag9> <gTitle9> <mTag9a=mTitle9a> [mTag9b= mTitle9b…]

But these is no such limitation (I guess so, not test yet) for the tag/title pairs in each subGroup. For example, ENCODE data trackDb put all TFs in one subGroup:
http://ftp.ebi.ac.uk/pub/databases/ensembl/encode/integration_data_jan2011/hg19/trackDb.txt

One question: How to share tracks but secure the data files in the track hub?

The current directory hierarchy for a hub is like:

myHub/ - directory containing track hub files

     hub.txt -  a short description of hub properties
     genomes.txt - list of genome assemblies included in the hub data
     hg19/ - directory of data for the hg19 (GRCh37) human assembly
          trackDb.txt - display properties for tracks in this directory
          dnase.html - description text for a DNase track 
          dnaseLiver.bigWig - wiggle plot of DNase in liver
          dnaseLiver.bigBed - regions of active DNase
          dnaseLung.bigWig - wiggle plot of DNase in lung
          dnaseLung.bigWig - regions of active DNase
          ...
          rnaSeq.html - description text for an RNAseq track
          rnaSeqLiver.bigWig - wiggle plot of RNAseq data in liver
          rnaSeqLiver.bigBed - intron/exon lists for liver
          rnaSeqLung.bigWig - wiggle plot of RNAseq data in lung
          rnaSeqLung.bigBed - intron/exon lists for lung
     hg18/ - directory of data for the hg18 (Build 36) human assembly
          trackDb.txt - display properties for tracks in this directory
          dnase.html - description text for a DNase track 
          dnaseLiver.bigWig - wiggle plot of DNase data in liver
          dnaseLiver.bigBed - regions of active DNase
          dnaseLung.bigWig - wiggle plot of DNase data in lung
          dnaseLung.bigWig - regions of active DNase
          ...
          rnaSeq.html - description text for an RNAseq track
          rnaSeqLiver.bigWig - wiggle plot of RNAseq data in liver
          rnaSeqLiver.bigBed - intron/exon lists for liver
          rnaSeqLung.bigWig - wiggle plot of RNAseq data in lung
          rnaSeqLung.bigBed - intron/exon lists for lung

The UCSC webpage also indicates that "unlisted hubs are in no way secure." But this is definitely a unsolved problem. Maybe the only solution is to set up your own local mirror?

How to use the "textHistogram -aveCol" and textHist2 in the UCSC executables

The Jim Kent's executable tools are general quite friendly to use. But for some options of some tools, it's a bit too terse. For example, the "-aveCol" option in textHistogram says

-aveCol=N - A second column to average over. The averages
will be output in place of counts of primary column.

Also, for textHist2, it says

textHist2 - Make two dimensional histogram table out

of a list of 2-D points, one per line.

usage:

   textHist2 input

options:

   -xBins=N - number of bins in x dimension

   -yBins=N - number of bins in y dimension

   -xBinSize=N - size of bins in x dimension

   -yBinSize=N - size of bins in x dimension

   -xMin=N - minimum x number to record

   -yMin=N - minimum y number to record

   -ps=output.ps - make PostScript output

   -psSize=N - Size in points (1/72th of inch)

   -labelStep=N - How many bins to skip between labels

   -margin=N - Margin in points for PostScript output

   -log    - Logarithmic output (only works with ps now)

   -postScale=N (default 1.000000) - What to scale by after normalization

What do this mean? For example, I have userFile with 2 columns like,
1 8
2 9
3 6
5 3

Here is reply from Brooke Rhead in UCSC Genome Bioinformatics Group:
The -aveCol option prints the average values of the items in each bin instead of the number of items in each bin. (The default bin size is 1.) So, if you run textHistogram on your example input without the -aveCol option and the default bin size, you get the number of items in each bin, which is 1 item for all bins except #4:

$ textHistogram userFile
1 ************************************************************ 1
2 ************************************************************ 1
3 ************************************************************ 1
4 0
5 ************************************************************ 1

If you instead use -aveCol=2, you get the average of the *value* in column 2. Since bin size is still equal to 1, this amounts to printing the value in column 2 in each bin:

$ textHistogram -aveCol=2 userFile
1 ***************************************************** 8.000000
2 ************************************************************ 9.00000
3 **************************************** 6.000000
4 0.000000
5 ******************** 3.000000

If you specify a smaller bin size, the option might seem more useful:

$ textHistogram -binSize=2 -aveCol=2 userFile
0 ************************************************************ 8.00000
2 ******************************************************** 7.500000
4 *********************** 3.000000

The bins are taken from the first column: 0-1, 2-3, and 4-5, and the values are from the second column. So, for instance, the bin labeled "2" above contains the average of the values from this part of your file:
2 9
3 6
. . . bin "2" contains the average of 9 and 6, or 7.5.


The textHist2 program makes more sense if you specify the number of bins and bin sizes to use. For example, if we specify the following options on your same input file:

1 8
2 9
3 6
5 3

we get:

$ textHist2 -xBins=6 -yBins=10 userFile
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 1 0 0
0 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0

The first column in the input file specifies the x-axis in the two dimensional histogram, and the second column specifies the y-axis. The numbering of the rows and columns starts at 0 in both cases. The first line of your input:

1 8

Causes a 1 to appear in the second column from the right and the 9th row from the bottom, or column 1 and row 8, if you start counting from zero. The "2 9" line causes the 1 to appear in column 2, row 9, and so on.

Basic knowledge for a bioinformatician

Very often, esp. when I was interviewed for a job or talk with a knowledgable guy like Xiaopeng, I feel there are full of "holes" in the mass body of my knowledge. How awkward it is! Guess I am not the only one who feels the same.

As a senior-in-age-but-not-senior-in-knowledge bioinformatian, I would seriously recommend who will like to work in this field to have basic knowledge in the following subjects I can think of:
1. probability and statistics (not everyone know the difference between them)
2. machine learning (the 4-elements circle: data + algorithm + model + criteria)
3. programming design (knowing how to write script does not mean you know how to program; a good programer should learn the concept of how to write code in a inheritable manner).
4. algorithm and data structure (many know some algorithm, but to truly understand it is not a easy task. Binindex is a good example of using the concept of binary tree to store/query genomic coordinate in a super fast way.)
5. know how to appreciate a scientific work. (A paper can be good in way of (i) data sources (2) method and/or (3) idea. For sure it's also important to tell good paper from junk papers. I feel it's so important to enhance the sensitivity of 'smelling' a paper)

I found this nice reading list from Hendrik's page (http://www.liacs.nl/~hoogeboo/mcb/nature_primer.html)

How to apply de Bruijn graphs to genome assembly (Phillip E C Compeau, Pavel A Pevzner & Glenn Tesler) November 2011, Vol 29, No 11; pp 987 - 991 doi: 10.1038/nbt.2023 (?) Analyzing 'omics data using hierarchical models (Hongkai Ji & X Shirley Liu) April 2010, Vol 28, No 4; pp 337 - 340 doi: 10.1038/nbt.1619 (?) What is flux balance analysis? (Jeffrey D Orth, Ines Thiele & Bernhard Ø Palsson) March 2010, Vol 28, No 3; pp 245 - 248 doi: 10.1038/nbt.1614 (?) How does multiple testing correction work? (William S Noble) December 2009, Vol 27, No 12 ; pp 1135 - 1137 doi: 10.1038/nbt1209-1135 (?) How to visually interpret biological data using networks (Daniele Merico, David Gfeller & Gary D Bader) October 2009, Vol 27 No 10 ; pp 921 - 924 doi: 10.1038/nbt.1567 (?) How to map billions of short reads onto genomes (Cole Trapnell & Steven L Salzberg) May 2009, Vol 27, No 5; pp 455 - 457 doi: 10.1038/nbt0509-455 (?) SNP imputation in association studies (Eran Halperin & Dietrich A Stephan) April 2009, Vol 27, No 4; pp 349 - 351 doi: 10.1038/nbt0409-349 (?) Maximizing power in association studies (Eran Halperin & Dietrich A Stephan) March 2009, Vol 27, No 3; pp 255 - 256 doi: 10.1038/nbt0309-255 (?) Understanding genome browsing (Melissa S Cline & W James Kent) February 2009, Vol 27, No 2; pp 153 - 155 doi: 10.1038/nbt0209-153 (?) What are decision trees? (Carl Kingsford & Steven L Salzberg) September 2008, Volume 26, No 9; pp 1011 - 1013 doi: 10.1038/nbt0908-1011 (?) What is the expectation maximization algorithm? (Chuong B Do & Serafim Batzoglou) August 2008, Volume 26 No 8; pp 897 - 899 doi: 10.1038/nbt1406 (?) What is principal component analysis?  (Markus Ringnér) March 2008, Volume 26, No 3; pp 303 - 304 doi: 10.1038/nbt0308-303 (?) What are artificial neural networks? (Anders Krogh) February 2008, Volume 26, No 2; pp 195 - 197 doi: 10.1038/nbt1386 (?)

How does eukaryotic gene prediction work?  (Michael R Brent) August 2007, Volume 25, No 8; pp 883 - 885 doi: 10.1038/nbt0807-883 (?) How do shotgun proteomics algorithms identify proteins? (Edward M Marcotte) July 2007, Volume 25, No 7; pp 755 - 757 doi: 10.1038/nbt0707-755 (?) What is a support vector machine? (William S Noble) December 2006, Volume 24, No 12; pp 1565 - 1567 doi: 10.1038/nbt1206-1565 (?) How does DNA sequence motif discovery work? (Patrik D'haeseleer) August 2006, Volume 24, No 8; pp 959 - 961 doi: 10.1038/nbt0806-959 (?) What are DNA sequence motifs? (Patrik D'haeseleer) April 2006, Volume 24, No 4; pp 423 - 425 doi: 10.1038/nbt0406-423 (?) Inference in Bayesian networks (Chris J Needham, James R Bradford, Andrew J Bulpitt & David R Westhead) January 2006, Volume 24, No 1; pp 51 - 53 doi: 10.1038/nbt0106-51 (?) How does gene expression clustering work?  (Patrik D'haeseleer) December 2005, Volume 23, No 12; pp 1499 - 1501 doi: 10.1038/nbt1205-1499 (?) How do RNA folding algorithms work?  (Sean R Eddy) November 2004, Volume 22, No 11; pp 1457 - 1458 doi: 10.1038/nbt1104-1457 (?) What is a hidden Markov model?  (Sean R Eddy) October 2004, Volume 22, No 10; pp 1315 - 1316 doi: 10.1038/nbt1004-1315 (?) What is Bayesian statistics?  (Sean R Eddy) September 2004, Volume 22, No 9; pp 1177 - 1178 doi: 10.1038/nbt0904-1177 (?) Where did the BLOSUM62 alignment score matrix come from?  (Sean R Eddy) August 2004, Volume 22, No 8; pp 1035 - 1036 doi: 10.1038/nbt0804-1035 (?) What is dynamic programming?  (Sean R Eddy) July 2004, Volume 22, No 7; pp 909 - 910 doi: 10.1038/nbt0704-909 (?)

Getting Started in ... Getting Started in Gene Orthology and Functional Analysis (Fang G, Bhardwaj N, Robilotto R, Gerstein MB) PLoS Comput Biol (2010) 6(3): e1000703; doi: 10.1371/journal.pcbi.1000703 (?) Getting Started in Structural Phylogenomics (Sjölander K ) PLoS Comput Biol (2010) 6(1): e1000621 ; doi: 10.1371/journal.pcbi.1000621 (?) Getting Started in Gene Expression Microarray Analysis (Slonim DK, Yanai I) PLoS Comput Biol (2009) 5(10): e1000543; doi: 10.1371/journal.pcbi.1000543 (?) Getting Started in Text Mining: Part Two.  (Rzhetsky A, Seringhaus M, Gerstein MB) PLoS Comput Biol (2009) 5(7): e1000411. ; doi: 10.1371/journal.pcbi.1000411 (?) Getting Started in Computational Mass Spectrometry-Based Proteomics.  (Vitek O) PLoS Comput Biol (2009) 5(5): e1000366. ; doi: 10.1371/journal.pcbi.1000366 (?)

Getting Started in Computational Immunology. (Kleinstein SH ) PLoS Comput Biol (2008) 4(8): e1000128; doi: 10.1371/journal.pcbi.1000128 (?) Getting Started in Biological Pathway Construction and Analysis.  (Viswanathan GA, Seto J, Patil S, Nudelman G, Sealfon SC ) PLoS Comput Biol (2008) 4(2): e16; doi: 10.1371/journal.pcbi.0040016 (?) Getting Started in Text Mining (Cohen KB, Hunter L) PLoS Comput Biol (2008) 4(1): e20; doi: 10.1371/journal.pcbi.0040020 (?) Getting Started in Probabilistic Graphical Models.  (Airoldi EM ) PLoS Comput Biol (2007) 3(12): e252. ; doi: 10.1371/journal.pcbi.0030252 (?) Getting Started in Tiling Microarray Analysis (Liu XS) PLoS Comput Biol (2007) 3(10): e183; doi: 10.1371/journal.pcbi.0030183 (?) Ten Simple Rules Also the Ten Simple Rules series of editorials has a separate page at the PLoS journal. A link is now all you need to read about 'Ten Simple Rules for Getting Published' or '...for a Good Poster Presentation', etc. On the Process of Becoming a Great Scientist (Giddings MC) PLoS Comput Biol (2008) 4(2): e33; doi: 10.1371/journal.pcbi.0040033 (?)