Rumenation

Conveniently skipping over the background of my current analysis for the time being -- to avoid a lengthy technical post -- last night I sent a batch job to our compute cluster to count how many times features of some interest (identified via BLAST) appeared over all reads (partial genomic sequences) in our assembly.

This morning I inspected the results of the htseq-count joblot and noticed a consistent error across my log files.

[msn@bert limpet_swissprot_bac]$ grep "^Warning" htseq-count.sh.e* | cut -f2- -d':' | sort -u Warning: 354465166 reads with missing mate encountered.

This is somewhat unsettling. 354465166 is a large number of things to have missing. We can confirm its largeness by adding some commas: 354,465,166.

But is this really a large number in context? The log files have a final line of output (assuming success) for the number of pairs encountered during parsing of the alignment BAM:

[msn@bert limpet_swissprot_bac]$ tail -q -n1 htseq-count.sh.e* | sort -u 363027622 SAM alignment pairs processed.

Confusingly, htseq-count believes to have encountered 363027622 pairs, whereas counting the number of reads (or records) with samtools view gives a different value:

[msn@bert limpet_vh_51_bt_bam]$ samtools view -c Limpet-Magda.sort.bam 365107681

So are they reads or pairs and why are they different? Well the source shows that line is written to stderr via the following:

sys.stderr.write( "%d SAM alignment record%s processed.\n" % ( i, "s" if not pe_mode else " pairs" ) )

htseq-count merely inserts the word "pairs" if it has detected the records encountered are from paired-end sequencing (pe_mode), even though the input BAM really does just represent one read per line. This seems a bit confusing if not outright misleading.

For the disagreement on how many records are there, my guess is htseq-count has ignored these ~2M reads on some other basis such as quality.

Now given the figures and the ambiguity of reads vs. pairs apparently solved, we have 365107681 total reads processed - 354465166 reads missing mate leaving a mere 10,642,515 records for us to play with. At that rate it would appear we've "lost" 97% of our data!

What is going on? bowtie2; the alignment software I used to map the assembled contigs back to the original raw reads estimated that 42.53% of the reads aligned discordantly, which may account for a considerable number of reads not having a mate. Yet not over 97%, surely?

I thought perhaps htseq-count intended to say 354465166 pairs missing mate or equivalently, 177232583 reads missing mate. Which would leave us with 48.54% reads missing a mate - which could seem almost reasonable given our prior knowledge about the rate of discordant reads that make up the assembly.

But at this point I have to point the finger at myself and wonder what I might have done incorrectly. I'll keep you posted.

**

Update I think the mess was caused when I used samtools sort to sort the alignment BAM without realising that by default alignments are organised by position (thus contig) rather than by readname. This would mean that our paired-end data would not be stored adjacently in the resultant file, which is what htseq-count expects by default!

The solution is either to sort by readname with samtools sort -n or to specify to htseq-count that the alignment file is sorted by position with htseq-count -r pos.

I opted for the former as the behaviour of htseq-count for the latter is rather non-optimal. As mate pairs are not stored adjacently, reads are stored in a large hash until the relevant mate is found -- adding significant overhead.

Having spent quite some time ruminating over the contents of this inaugural blog post to a point beyond that of analysis paralysis, I have managed to avoid writing entirely.

Thus I'll remove this pedestal of preclusion and begin with an uninspiringly titled First Post.

**

Hello, I'm Sam Nicholls; one of the latest PhD students to appear on the roster at Aberystwyth University, enrolled jointly between the Department of Computer Science and the Institute of Biological, Environmental and Rural Sciences (colloquially known as IBERS).

This inter-departmental, inter-institute collaboration has me well placed for a project that sits on the boundaries of both computing and biology and whilst doubling the supervisor count and the administrative paperwork, I have two different coffee rooms to hide in and get my name on two distinct doors.

The project, provisionally entitled Metagenomics of the Rumen Microbiome -- which I'll get to in a moment -- is still being defined. Initially this was somewhat concerning but I quickly garnered from my admittedly short interactions with other PhD students that this is very much the norm at such an early stage. Phew.

I hope to provide a better overview of the project in a future post but suffice to say for now we're looking for "interesting" micro-organisms that live in the rumen: the first stop of the four-chambered stomach that grass and feed will take on their tour through the digestive system of a cow (or other ruminant animals like sheep, goats and giraffes).

The rumen is a complex and diverse ecosystem inhabited by bacteria (Earth's first life, abundant simple single celled micro-organisms found everywhere in their millions), archaea (micro-organisms that were initially thought to thrive only in extreme environments -- with an evolution influenced by an early Earth -- but more recently discovered in non-extreme environments too, play a similar role to bacteria ecologically), fungi (a broad class including single-celled yeasts, moulds and complex multi-cellular forms such as mushrooms), protozoa (single celled micro-organisms that feed from their direct surroundings and have the capacity for controlled movement with a tendency to thrive in moist environments) and viruses (the smallest of the microbes, which can only survive inside the cells of other living things). An unsuspecting cow yields anything from 20 to 50 billion bacteria in the volume of a raindrop alone.

But why is what goes on in a cow's stomach of any interest at all? These rumen residents are adapted to efficiently break down biomass for their host's continuing survival. The cell walls of grasses and other green plants for example, feature cellulose -- long chains of glucose molecules which are difficult to digest -- as a major structural component. Yet species of bacteria found in the rumen produce enzymes (catalysts for biochemical reactions) capable of breaking up these long chains in to smaller chains that are more easily further broken apart and used for energy. If we could isolate such enzymes we may be able to leverage them to improve the efficiency of biofuels like cellulosic ethanol.

This diverse and densely populated environment necessitates strong competition from its denizens to fight for their own survival. Under times of stress and competition for resources, microbes release potent antimicrobial peptides (very short proteins) to destroy their rivals.

Antibiotics are essential medicines for the treatment of bacterial infections in humans and animals. Bacteria presenting resistance to available antibiotics pose a serious problem in the world of medicine and the appearance of highly resistant and life threatening "superbugs" like MRSA are becoming increasingly commonplace in hospitals the world over. Discoveries of antimicrobials have dwindled over the past decade despite a clear medical need to synthesize new drugs to fight resistant bacteria. Yet in our unsuspecting cow's paunch, a puzzle of peptides that could potentially lead to humanity's next generation of antibiotics is ripe for exploration.

Another interesting possibility for analysis is working to lower methane emissions from cattle. Livestock are a significant global contributor to greenhouse gases, more so than transportation. The average cow will belch out around 500 litres of methane per day -- a greenhouse gas which traps almost 30 times more heat than carbon dioxide -- due to its biomass diet. By investigating the effects of various crops and feed on the rumen and the resulting interactions between the various microbes within, it might be possible to make a more environmentally friendly cow.

Finally, a methodology (both in terms of knowledge and software) for understanding the processes of such a complex biome could be applied and scaled to other microbiosystems such as that of our own gut, a world we know quite little about despite its considerable influence on our health and wellbeing.

This of course is all quite far off and currently I'm just finding my feet and when not occupied by compulsory soft-skills training, trying to absorb as much new biology terminology and techniques as possible -- especially having had such little prior experience with biology coming from my computer science and statistics background.

I hope to go in to more detail shortly about what I've been doing so far and why trying to answer these questions pose a difficult computational problem -- hopefully validating why we can make a PhD of it -- as it would certainly make for a full post of its own.