Genome Spot

The genome sequencing field is undoubtedly dominated by Illumina, and now they have two new sequencing platforms to offer ( link ). The HiSeq X Ten is aimed at large sequencing centres and Illumina claims that it can sequence genomes for $1000. Very similar to the Hiseq2500 in most regards, it doesn't offer Rapid Run mode (2 lanes) but has 4x faster run times in High Capacity mode (8 lanes). Only catch is the price, at USD$10M for a minimum order of 10 instruments, this isn't for labs without huge resources. The other instrument released was the NextSeq500 (info here ). It runs a 4 lane flow cell capable of generating 400M reads. NextSeq is aimed at labs which do RNA-Seq and similar protocols but can't afford the price tag of their big brother HiSeq machines. Given the history of Illuminas platforms, I'd expect the NextSeq to undergo further incremental increases in output as the platform matures. Priced at about $250,000 it is placed in between the MiSeq and the Hi...

MSigDB contains thousands of gene sets which have been mined from a range of genome wide studies and these are a valuable resource for gene ontology and pathway analysis. You probably know that many gene sets are curated by KEGG, REACTOME and BIOCARTA - as well as dry lab scientists who specialise in analysing these data sets and curating gene sets. What you may not know is that if you follow a few basic guidelines, you can start generating your own custom gene sets and these can become a valuable resource for running gene ontology and Gene Set Enrichment Analysis (as per the graphic) For instance within our lab, we have extensively used ENCODE ChIP-Seq data to help us to analyse our mRNA-seq data and this has provided a huge leg-up in generating hypothesis and designing follow-up experiments. For this example, I want to show an overview of how I made biomarker gene sets for a bunch of cell types analysed by ENCODE. Biomarker gene sets can useful in array or mRNA-Seq analysis tha...

Here's  a selection of one liners that you might find useful in data analysis. I'll update this post regularly and hope you can share some of your own gems. Compress an archive: tar cf backup.tar.bz2 --use-compress-prog=pbzip2  directory Sync a folder of data to a backup rsync -azhv /scratch/project1 /backup/project1/ Extract a tar archive: tar xvf backup.tar   Head a compressed file: bzcat file.bz2 | head   pbzip2 -dc file.bz2 | head zcat file.gz | head  Uniq on one column (field 2): awk '!arr[$2]++' file.txt Explode a file on the first field: awk '{print > $1}' file1.txt Sum a column: awk '{ sum+=$1} END {print sum}' file.txt Put spaces between every three characters sed 's/\(.\{3\}\)/\1 /g' file.txt Join 2 files based on field 1. Both files need to be properly sorted (use sort -k 1b,1)  join -1 1 -2 1 file1.txt file2.txt Join a bunch of files by field 1. Individual files don't need to be sorted but the ...

In my previous post , I discussed an efficient approach to align deep sequencing reads to a genome, which could be used for most applications (mRNA-Seq, ChIP-Seq, exome-Seq, etc). MicroRNA datasets are rather a bit different because they often contain millions of reads from only a few thousand RNA species, meaning there is a lot of redundancy. We can speed up this analysis considerably if we collapse the sequences and remove the redundancy before genome alignment. You can see in the below script that several steps are done without writing a single intermediate file: Adapter clipping Sequence collapsing BWA alignment Bam conversion Sorting bam file for FQ in *fq do fastx_clipper -a TGGAATTCTCGGGTGCCAAGG -l 18 -i $FQ \ | fastx_collapser \ | bwa aln -t 30 $REF - \ | bwa samse $REF - ${FQ} \ | samtools view -uSh - \ | samtools sort - ${FQ}.sort & done wait for bam in *bam ; do samtools index $bam & done ; wait for bam in *bam ; do samtools flagstat $bam > ${bam}.sta...

We know that sequencing platforms are generating more and more output, and more of this data is finding itself in online databanks. The volume of data is a blessing because it is so information rich, but it is also increasingly becoming a problem - we need more machine time to process it. I was recently aligning dozens of mouse encode RNA-seq data sets while we had a server fault. I was left with a mess of half finished jobs and intermediate files and basically had to start the whole process again. This prompted me to make an alignment pipeline that would be more resistant to unexpected interruptions. To do this, I took inspiration from a post at SEQanswers which piped the output from the BWA aligner and generated bam files without any intermediate files. But I wanted to take this a bit further, by going from compressed fastq to bam file without writing any intermediate files, doing quality trimming on the fly and saving a lot of storage in the process. So here are two examples, o...

Pathway and gene set analysis is a common procedure for interpretation of RNA-seq or other genome-wide expression assays. Most of the time, we use GSEA to tell us whether our gene sets of interest are up- or down-regulated. We can use gene sets from KEGG , Reactome , GO ,  MSigDB  and other sources, but you can also generate your own gene sets. The format used for GSEA is gmt . I'm going to take you through two examples of generating custom gene sets: Generate gene sets from published data sets using GEO2R Let's say you're interested in the transcription factor STAT1. I found a dataset in GEO called "Knockdown of STAT1 in SCC61 tumor xenografts leads to alterations in the expression of energy metabolic pathways", which has a paper in BMC Med.  Most uploaded array data sets can be reanalysed with GEO2R, which runs the array analysis tool Limma but this is embedded in the webpage and has a GUI which makes it very accessible for biologists. Click this link ...

In the previous post , I mentioned that in order to do counting, we need some regions (exons) to count over. We generated a BED file based on GENCODE exons which looks like this: chr9 33248310 33248419 ENSG00000107262.11_BAG1 chr9 33252477 33255306 ENSG00000107262.11_BAG1 chr9 33255862 33255925 ENSG00000107262.11_BAG1 Now we can actually count how many reads from each of our samples maps to these exons. To do that there are a few different programs: htseq-count samtools view GATK countreads Bams2mx BedTools  multicov For this demo, I will outline use of BedTools for RNA-Seq. I think BedTools nulticov is a good option for 3 main reasons: It can handle multiple bam files at the same time, producing a count table/matrix which is pretty much ready for statistical analysis. It is aware of mapQ scores of alignments, meaning that if the program comes across non-unique alignments, you can set it to ignore or include those reads. The BedTools package is well ...