Genome Spot

I've been a fan of GNU parallel for a long time. Initially I was sceptical about using it, preferring to write huge for loops but over time I've grown to love it. The beauty of GNU parallel is that it spawns a specified number of jobs in parallel and then submits more jobs as others are completed. This means that you get maximum usage out of the CPUs without overloading the system. There are many excuses for not using it , but perhaps the only valid one is that you have Sun Grid Engine or another job scheduler or manager in place. GNU parallel is particularly useful when used with functions . Functions are subroutines that may be repeated many times to complete a piece of work. In bash, here is a simple example , which declares a function consisting of a chain of piped commands, and then executes 4 jobs in parallel, until all of *files.txt have been processed. #!/bin/bash my_func2() { INPUT=$1 VAR1=bar cmd1 $INPUT $VAR1 | cmd2 | cmd3 > ${1}.out } export -...

CpG shores are the regions immediately flanking and up to 2 kbp away from CpG islands. These regions are interesting because methylation they are variably methylated in cancer and development (see reference ). Image courtesy of 7-themes.com To identify CpG shores, you first need a list of coordinates of CpG islands. If you have a linux computer with Emboss tools installed, you can determine CpG island positions with the cpgplot command . (if you don't have Emboss installed you can get CpG island coordinates from the UCSC table browser  just be wary of chromosome naming ie, "chr1" vs "1") cpgreport -score 17 -outfile genome.fa.cpg -outfeat /dev/stdout genome.fa | head Will produce a file like this one: ##gff-version 3 ##sequence-region 1 1 59999934 #!Date 2015-04-10 #!Type DNA #!Source-version EMBOSS 6.6.0.0 1 cpgreport sequence_feature 10469 11240 1299 + . ID=1.1 1 cpgreport sequence_feature 11284 11301 37 + . ID=1.2 1 cpgrepo...

Repetitive sequences pose significant challenges for the analysis of short read sequence data. Any sequence read that is shorter than the repeat unit is going to fail at being uniquely placed in the genome. At a given read length, genome can therefore be divided into mappable and unmappable compartments. For analyses such as gene expression and chromatin profiling, we generally ignore those unmappable regions and ignore reads that map to multiple positions. But what if we were actually interested in the abundance of repeat sequences in the dataset? How could it be quantified with short read sequencing? Why is that even relevant? The "Why" About 50% of the human genome is comprised of repetitive DNA ( Ref ). While the function of these elements in many cases remains unknown, some have been well defined. Alu elements have been shown to be a birthplace of enhancers ( Ref ). Endogenous retroviruses have been domesticated to boost antiviral response in B cells after exposure t...

Quality trimming and adapter clipping paired end reads is a tricky business. For paired-end alignments to work, the order of the sequences in both files (forward and reverse) needs to be retained. While there are plenty of read trimmers available open-source (think FASTX-toolkit, SeqTK, CutAdapt, etc), I haven't found many that: Retain pairing info Run parallel and are fast Are easy to set up and run Have good docs Then I fell in love with Skewer . It smashed through 14.5 million gzipped read pairs, doing adapter clipping and quality trimming in two minutes on my 8-core workstation. It auto-detects quality encoding so you can safely analyse any Illumina data. Awesome! $ skewer -t 8 -q 20 SRR634969.sra_1.fastq.gz SRR634969.sra_2.fastq.gz Parameters used: -- 3' end adapter sequence (-x): AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -- paired 3' end adapter sequence (-y): AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA -- maximum error ratio allowed (-r): 0.100 -- maximum indel er...

HISAT is a brand new RNA-seq aligner which promises great speed with a low memory footprint. The Nature Methods paper is worth a look. So I thought I'd give it a test run with some simulated data to check its accuracy compared to other aligners. I generated synthetic 100bp reads based on Arabidopsis cDNAs and then incorporated mutations with msbar . I then aligned these reads to the Arabidopsis genome with default settings and counted ( featureCounts ) the number of correctly and incorrectly assigned reads with a mapping quality of 20. Here is the result for unmutated (perfect) 100 bp reads. Table 1. Alignment of simulated 100 bp cDNA sequences to the Arabidopsis genome. HISAT shows similar accuracy to TopHat2 but was not as accurate as STAR. Now here are the results of the mutation experiment. Figure 1. Accuracy of mapping mutated 100bp reads. Left hand side graphs show the correct mapping rates and right hand side shows incorrect mapping rates. Top panels ...

When analyzing ChIP-seq data, we sometimes like to focus on promoters as a proxy of gene activation especially for histone marks that are normally associated with gene activation such as H3K9 acetylation or H3K4 tri-methylation. FeatureCounts has emerged as a competitor to HTSeq and BedTools MultiCov for counting reads across features (ie, exons, genes, promoters). FeatureCounts is great for RNA-seq because it can natively read GTF annotation files, but can't read BED format (that we use a lot in ChIP-seq analysis). In order to make featureCounts work, we need to extract the TSS coordinates and convert to a BED-like format that it can read ( SAF format ). In the below script, I extract the positions of the TSSs using a grep search followed by stripping only the neccessary information from the GTF, then using awk and BedTools to expand the region  around the TSS by 3kbp. #!/bin/bash #Generate an SAF file for TSS intervals from a GTF #Specify some parameters GTF=Homo_sapien...

Little known UNIX features to avoid writing temporary files in your data pipelines explained by Vince Buffalo in his digital notebook . Introducing named pipes and process substitution.