samtools upgrade for exome & target region SNP calling

off the samtools helplist, a new feature that's very helpful for exome / target region NGS studies

BTW, samtools is now able to compute read depth and call SNPs in regions specified by an input BED file. For both "samtools mpileup" (not pileup) and "bcftools view", you may provide the BED via the "-l" option. If the input has two numeric columns, it is parsed as a BED (region list); if has one numeric column, parsed as a position list file, so the "-l" option is backward compatible with old versions. It is also possible to retrieve alignments overlapping a BED file with 

samtools view -L in.bed 

. For mpileup, using "-l" to call SNPs in target/exome regions can be much faster than doing whole-genome calling and then filtering.

Heng

Convert SAM / BAM to fasta / fastq

Probably one of the most freq FAQ

latest thread in biostar
http://biostar.stackexchange.com/questions/6993/convert-bam-file-to-fasta-file

Samtofastq using Picard
http://picard.sourceforge.net/command-line-overview.shtml#SamToFastq

Samtools and awk to make fasta from sam

samtools view filename.bam | awk '{OFS="\t"; print ">"$1"\n"$10}' - > filename.fasta

Biopython and pysam (code contributed by Brad Chapman)

http://biostar.stackexchange.com/questions/6993/convert-bam-file-to-fasta-file/6994#6994 

A practical comparison of de novo genome assembly software tools for next-generation sequencing technologies.

    A practical comparison of de novo genome assembly software tools for next-generation sequencing technologies.
    Zhang W, Chen J, Yang Y, Tang Y, Shang J, Shen B.
    PLoS One. 2011 Mar 14;6(3):e17915.
    PMID: 21423806 [PubMed - in process]

Should be interesting..

de novo assembly of Illumina CEO genome in 11.5 h - new ver of Ray

Kevin:You can't ignore an email with that subject header.. but 512 compute cores? Shall have a chat with my HPC vendor.. 
Also am waiting for public release of Cortex http://sourceforge.net/projects/cortexassembler/
Strange that courses that teach the software are available but the software ain't ... 
http://www.ebi.ac.uk/training/onsite/NGS_120510.html


Velvet and Curtain seems promising for de novo assembly as well.

Ray 1.3.0 is now available online.
http://sourceforge.net/projects/denovoassembler/files/Ray-1.3.0.tar.bz2

The most important change is the correction of a major bug that caused
parallel infinite loop on the human genome.

This, along concepts incorporated in Ray 1.2.4, allowed Ray to assemble
the genome of Illumina's CEO in 11.5 hours using 512 compute cores (see
below for the link).

What's new?

1.3.0

2011-03-22

   * Vertices with less than 1 of coverage are ignored during the
computation of seeds and during the computation of extensions.
   * Computation of library outer distances relies on the virtual
communicator.
   * Expiry positions are used to toss away reads that are out-of-range
   * When only one choice is given during the extension and some reads
are in-range, then the sole choice is picked up.
   * Fixed a bug for empty reads.
   * A read is not added in the active set if it is marked on a
repeated vertex and its mate was not encountered yet.
   * Grouped messages in the extension of seeds.
   * Reads marked on repeated vertices are cached during the extension.
   * Paths are cached in the computation of fusions.
   * Fixed an infinite loop in the extension of seeds.
   * When fetching read markers for a vertex, send a list of mates to
meet if the vertex is repeated in order to reduce the communication.
   * Updated the Instruction Manual
   * Added a version of the logo without text.


I fixed a bug that caused an infinite loop. Now Ray can assemble large
genomes. See my blog post for more detail about that.
http://dskernel.blogspot.com/2011/03/de-novo-assembly-of-illumina-ceo-genome.html


Version 1.2.4 of Ray incorporated also new concepts that I will present
at RECOMB-Seq 2011.

The talk is available online:
http://boisvert.info/dropbox/recomb-seq-2011-talk.pdf


Sébastien Boisvert

Cheat Sheets Galore-bioinformatics, biology, linux,perl, python, R

started with Keith's post here
http://omicsomics.blogspot.com/2011/03/whats-on-your-cheat-sheet.html

and a thread at
http://biostar.stackexchange.com/questions/6683/bioinformatics-cheat-sheet


I have soooo many of them! *this is going to be a long post
Vim
http://www.viemu.com/vi-vim-cheat-sheet.gif
Python
http://www.addedbytes.com/cheat-sheets/python-cheat-sheet/
R (pdf)
http://cran.r-project.org/doc/contrib/Short-refcard.pdf
Perl 
Hmmmm where did that go to?
AWK one liners
Sed examples
Linux common tasks


I have these too 

• IUPAC ambiguity codes for nucleotides:
• Amino acid single letter codes.

FAQ-SNPs missing when called with more samples

Question
Using mpileup called with 2 different samples, Im getting this
particular SNP which has a gd coverage (DP=55) only in one particular
sample. This is good.

However when mpileup was called with 10 samples, the SNP got lost. Im
just trying to figure out if the SNP got 'drowned' out by the other 9
samples which doesnt have the SNP and hence, wasnt called. Is this how
mpileup works

Excellent answer by Heng Li 

With more samples, you gain power on SNPs shared between samples, but lose power on singleton SNPs. Here is a way of thinking of this. Suppose we have 1% false positive rate (FPR) for one sample. If we call SNPs from 10 samples separately and then combine the calls, the FPR would be around 5% (not 10% because more SNPs are found given 10 samples). To retain a low FPR on singletons we have to be more stringent. Nonetheless, with more samples, we can usually get overall better calls than calling SNPs in each sample separately because information between samples is used more effectively.

Common numbers / statistics for Uniquely mapped reads?

Was asked if there was a commonly reported numbers for
uniquely mapped reads (which is troublesome to define with bowtie)
vs
total mapped reads

Not sure also if the numbers differ for applications
e.g.
WGS
exome reseq

human vs other organisms.
Got this figure from a 2009 paper. Not sure if anyone collates data like this
http://bioinformatics.oxfordjournals.org/content/25/7/969.full.pdf