Genome Spot

You have probably seen word clouds before - but have you tried with gene expression data? I used the following bash script to process a DESeq spreadsheet from a a previous RNA-seq  post . The script extracts the gene name and p-value of the genes with differential expression. I used awk to separate the up and down regulated genes into different files. The score used to inform the font size is the exponent of the p-value. So this works best when there are a lot of statistically significant genes p-values. The data looks like this: AC011899.9:60 CAMK1D:54 GNG4:44 CGNL1:41 APCDD1:37 HSD11B1:33 Now go to the "advanced" tab of the Wordle page  and paste in the data. Experiment, with the colours, layouts and be sure to increase the "maximum words" to get a real appreciation for the number of changes in your experiment. Here is an example I made using both the up and down regulated gene sets showing the effect of azacytidine on AML3 cells. The result is pretty

In the previous post we looked at how choice of RNA-seq aligner influenced results. In this post, we'll use simulated reads to empirically determine the accuracy of OLego, STAR, SubJunc and SubRead. I made a pretty basic bash script (see below) to generate fasta format reads at uniform intervals along the length of transcripts without errors. The user can readily change the read length and the interval between read initiation. The script created 11.3 million 100 bp reads in 9 minutes so it is OK in terms of speed. It is designed to take Ensembl cDNA files as an input and has been tested on human  and  Arabidopsis . Here is what the generated reads look like >AT2G01210.1:gene:AT2G01210_197863 CGTCAGCTTTCGTTCTGGGGAAGAGCGGAATCGGAATTGTCTACAAAGTG >AT2G01210.1:gene:AT2G01210_197864 GTTCTAGAGAACGGGCTCACACTGGCCGTACGGAGATTGGGTGAAGGAGG >AT2G01210.1:gene:AT2G01210_197865 GTCTCAGAGATTCAAGGAGTTTCAGACAGAAGTTGAAGCCATAGGGAAAC >AT2G01210.1:gene:AT2G01210_197866 TAAAACATCCTAA

If you like science, you will love Forty-Four and Two . A new blog that delves into genetics and molecular biology for the wider audience. His latest post is an impressive poster describing some data we published together earlier this year in CMLS . Keep up the good work!

In the previous post , I compared the speed of several RNA aligners. In this post, we'll take a closer look at the results generated by these different aligners; Subread/Subjunc , STAR , HPG aligner and Olego . In addition, we will use BWA MEM as an example of what happens when you use a DNA aligner to analyse RNA-seq. For the test, we've been using 101 bp  Arabidopsis mRNA-seq data from GEO accession  GSE42968 . For simplicity, I use only the forward read in the paired-end dataset. I used fastx-toolkit to remove low quality bases from the 3' end (base qual threshold of 20). One of the simplest ways to assess the quality of an alignment is to determine the proportion of reads that are mapped to the genome and the proportion that map to exons. After the aligners did their work, I used featureCounts to quantify the number of aligned reads with a mapping quality >10. Here is the data for the first sample in the series, SRR634969 which contained 14.5 million reads.

In this post I compare the speed of Subread/Subjunc , STAR , HPG aligner and Olego . I will use real RNA-seq data from GEO accession  GSE42968  and align to the Arabidopsis thaliana genome . All code used is found at the bottom of this post. The benchmarking was performed on a standard 8 core workstation with 8 GB RAM. Index the Arabidopsis thaliana genome The first step in any RNA-seq analysis is to prepare the genome for alignment. Subread was the quickest. Time required to index the Arabidopsis genome (119.67 Mbp). Align RNA-seq reads I down sampled the sequence (SRA accesion:SRR634969) file to generate files containing 1M, 2M, 5M, 10M and 14M reads selected randomly. Fastq quality trimmer was used to trim low quality bases from the 3' end and filter low quality reads. Then I used Olego, Subread, Subjunc, STAR and HPG aligner to align the reads. (I wanted to test HPG aligner in both the BWT and SA modes but didn't have the 17 GB RAM to test SA algorithm sadly.)