Making coding T cell receptor sequences from V-J-CDR3

If, like me, you work on T cell receptors, occasionally you’re probably going to want to express a particular TCR in cells. However you're not always going to have the sequence, at either nucleotide or protein sequence level.

Not a problem, you can sort this out. You can look up all the relevant germline sequences from IMGT, trim away all the non-used bits, add in the non-templated stuff, then manually stitch it all together and have a look to see if it still makes what you were expecting. You can do all that... or you can just use the code I wrote.

StiTChR does it all: give it a V gene, a J gene, and a CDR3 amino acid sequence and it'll look up, trim and stitch together all the relevant TCR nucleotide sequences for you, back-translating the non-templated region using the most frequent codon per residue. It also translates it all, and will run a quick rudimentary alignment against a known partial protein sequence if you have one for a visual confirmation that it's made the right thing.

You can then take the alpha/beta TCR sequences it generates, bang them into an expression vector (typically split by a T2A sequence or something) and transduce your cells of interest.

I wrote this code to save me a bit of time in future, but hopefully it can do the same for some of you!

Bulk downloading proteome files from UniProt using Python

It's that time again, where the following has happened:

1. I want to do some niche bioinformatics related thing
2. I cobble together a quick script to do said thing
3. I throw that script up on the internet on the offchance it will save someone else the time of doing 2

It's a little shift of target and scale from a similar previous post (in which I used Python to extract specific DNA sequences from UCSC. This time I've been downloading a large number of proteome files from UniProt.

It's all explained in the docstring, but the basic idea is that you go on UniProt, search for the proteomes you want, and use their export tool to download tsv files containing the unique accession numbers with identify the data you're after. Then you simply run this script in the same directory; it takes those accessions, turns them in to URLs, downloads the FASTA data at that address and outputs it to new FASTA files on your computer, with separate files named after whatever the tsv files were named.

The best thing about this is you can download multiple different lists of accessions, and have them output to separate files. Say maybe you have a range of pathogens you're interesting in, each with multiple proteomes banked; this way you end up with one FASTA file for each, containing as many of their proteomes as you felt like including in your search.

Extracting antigen-specific TCRs from vdjdb on the command line

My previous lab and I recently published a review on the techniques and possibilities of analysing TCR repertoire data produced from high-throughput sequencing.

By and large I'm exceedingly happy with it - apart from a couple of missed references and one very unfortunate mix up regarding the accessibility, I'm very pleased with how it came out (and hope it will prove useful!).

One thing that I'm particularly pleased that we included (in spite of the lack of published descriptions yet) is the pair of manually curated TCR databases that have recently emerged: VDJdb and McPAS-TCR, in which you can find a small (but growing) host of TCRs of known specificity and/or disease association. We thought it was important to get these out there as soon as possible, as this is a rapidly changing field which is currently sorely needing for such efforts at standardisation and resource development.

With that in mind, I've been playing around with both of these, and thought I'd share some of the bare bones of the bash code I've been using to pull out sequences related to epitopes I'm interested in. Here is my quick vignette using VDJdb to pull out HIV-reactive TCR sequences - and even then just the fields of the database I'm interested in - using basically just default terminal commands.

Download specific DNA sequences from hg19 using Python

I've been working on a little side-project recently that involved needing to grab lots of different human DNA sequences based on their position, which lead me to discover the wonderful UCSC DAS server (from this informative Biostars thread).

Seeing as the rest of the project was written in Python, I knocked together a quick function to do just that. It's all nice and easy: just give it the chromosome number/letter*, and a numerical start and stop position, and the function returns the hg19 DNA sequence in that range.

I'm also trying to make a bit more use of GitHub (including knocking together a place for my publications), so I thought this was the perfect thing to make a gist from:

* Currently this function won't be able to grab anything from the unassigned chromosome contigs - just chromosomes 1-22, X, Y and mitochondrial (M) sequences.

Installing Trinity on Mac OS X via homebrew -- update

A couple of months ago I wrote a short post about how I managed to install the RNA assembler Trinity on Mac OS X (El Capitan), on the off-chance it would be useful to someone else.

This morning I received an email from my friend Mazlina, who I worked with in London, saying she had been trying to do just that and had coincidentally stumbled on my post*. However it hadn't worked out quite as easily for her as it had for me.

It turned out to be due to a Java version problem. While 1.8 was installed, brew --config claimed that only 1.6 was, which is insufficient for Trinity installation.

Here's how she solved it, quoting from her email:

"First, I did

brew doctor

and just cleared up whatever it told me to [...]

$ brew doctor

Your system is ready to brew.

So by running

brew cask search java

it lists down the available java versions (the one you have, 8, is just java), and I went with 9-beta because that was the only one I could download at that time. And when I ran 

$ brew cask install java9-beta

$ brew install trinity

worked like a charm.

Not sure if that's exciting enough to go on the blog, but it solves it anyway."

(For reference, as I told her, given the average excitement level of the blog this should fit right in!)

* I'm never really sure whether these posts are read or not (as the stats from Blogger are always inflated by scanning bots) apart from when people let me know they've seen them - so if you ever bump in to me at a conference or something and have found one of them useful or interesting please let me know, I love to hear about it!

Installing Trinity on Mac OS X

One of the inevitable joys of bioinformatic life is the installation of a variety of esoteric softwares on a variety of system. As I've just moved to a new position in a new institution, I get to go through this rigmarole again.

This time around I have an extra layer of faffery, as I am now for the first time using a Mac (having been on Ubuntu for the last ten years, and Windows in the distant recollections from before that). While the machine is gorgeous and responsive, I am still in the interminable murky phase where I don't know the intricacies and easy ways of doing things yet (and am still battling muscle memory for keyboard shortcuts!), which means that I'm back down the learning curve a little.

Anyway, as I've just discovered an incredibly easy way to install a very useful tool, I thought I'd share it.

I was installing the excellent RNA assembler Trinity on my iMac running OS X (El Capitan), or at least trying to, according to its website. However, despite attempts at using different (and newer) compilers, I kept running into this error, presumably reflecting my attempts at using alternative compilers failing:

clang: error: unsupported option '-fopenmp' trinity mac

Happily it turns out that Trinity is supported by the fantastic third party package manager homebrew, which I had coincidentally just installed anyway (you don't bundle wget in, what the heck Apple?).

Homebrew is easily installed following the details on their website, and then installing Trinity was as simple as this:

brew cask install java
brew install homebrew/science/trinity

Not only was this dead simple, but it automatically installed a number of other programs (as dependencies of Trinity) that were on my list to install anyway (e.g. trimmomatic and bamtools). It also installs everything directly to /usr/local/bin/, so there's no mucking about with your PATH required. Lovely.

NB: Whilst looking around for hints as to how to solve this problem, I did find this thread on SeqAnswers which suggests that you might need to take a little extra care when running Trinity on Mac systems as opposed to Linux. Something to bear in mind.

Count how many MiSeq reads derived from each surface of the flowcell

I recently had call to perform one of those tasks that I think others might, yet not be entirely sure how to go about it.

Specifically, in troubleshooting a MiSeq run's poor yield, I wanted to see whether there were significantly more reads derived from one of the flow cell surfaces (top or bottom) relative to the other. The reason I did this was my FWHM (full cluster width at half maximum, a measure of the focus during imaging) was noticeably higher for that surface.

I mean, I have no idea if ~3 is that much worse than ~2.8-2.9, but there's no harm in checking right?

This is very easily achieved as all of the information required to work it out is contained within the FASTQ reads themselves, in tile section of the identifier line of each each.

Therefore with a quick bit of basic bash we can find out exactly how many reads derived from each surface.

# Get all index reads (as the shortest) in one file
zcat *I1*z > I1.fq

# Extract the identifier lines with sed
 # and grep for those with a '1' at the right position
 # This indicated they derived from the top surface
sed '2~4d;3~4d;4~4d' I1.fq | grep ^.............................1 -c

# Do the same for '2', i.e. the bottom surface
sed '2~4d;3~4d;4~4d' I1.fq | grep ^.............................2 -c

And there you have it. Simple, quick and effective.

(As it turned out I have almost equal numbers derived from both surfaces, so it wasn't to blame in my case, but this might be useful for other situations!)

Heterogeneity in the polymerase chain reaction

I've touched briefly on some of the insights I made writing my thesis in a previous blog post. The other thing I've been doing a lot of over the last year or so is writing and contributing to papers. I've been thinking that it might be nice to write a few little blog posts on these, to give a little background information on the papers themselves, and maybe (in the theme of this blog) share a little insight into the processes that went into making them.

The paper I'll cover in this piece was published in Scientific Reports in October. I won't go into great detail on this one, not least because I'm only a (actually, the) middle author on it: this was primarily the excellent work of my friends and colleagues Katharine Best and Theres Oakes, who performed the bulk of the analysis and wet-lab work respectively (although I also did a little of both). Also, our supervisor Benny Chain summarised the findings of the article itself on his own blog, which covers the principles very succinctly.

Instead, I thought I'd write this blog to share that piece of information that I always wonder about when I read a paper: what made them look at this, what put them on this path? This is where I think I made my major contribution to this paper, as (based on my recollections) it began with observations made during my PhD.

My PhD primarily dealt with the development and application of deep-sequencing protocols for measuring T-cell receptor (TCR) repertoires (which, when I started, there were not many published protocols for). As a part of optimising our library preparation strategies, I thought that we might use random nucleotide sequences in our PCR products – which were originally added to increase diversity, overcoming a limitation in the Illumina sequencing technology – to act as unique molecular barcodes. Basically, adding random sequences to our target DNA before amplification uniquely labels each molecule. Then, in the final data we can infer that any matching sequences that share the same barcode are probably just PCR duplicates, if we have enough random barcodes*, meaning that sequence was less prevalent in the original sample than one might think based on raw read counts. Not only does this provide better quantitative data, but by looking to see whether different sequences share a barcode we can find likely erroneous sequences produced during PCR or sequencing, improving the qualitative aspects of the data as well. Therefore we thought (and still do!) that we were on to a good thing.

(Please note that we are not saying that we invented this, just that we have done it: it has of course been done before, both in RNA-seq (e.g. Fu et al, 2011 and Shiroguchi et al, 2012) at large and in variable antigen receptor sequencing (Weinstein et al, 2009), but it certainly wasn't widespread at the time; indeed there's really only one other lab I know of even now that's doing it (Shugay et al, 2014).)

However, in writing the scripts to 'collapse' the data (i.e. remove artificial sequence inflation due to PCR amplification, and throw out erroneous sequences) I noticed that the degree to which different TCR sequences were amplified followed an interesting distribution:

Here I've plotted the raw, uncollapsed frequency of a given TCR sequence (i.e. the number of reads containing that TCR, here slightly inaccurately labelled 'clonal frequency') against that value divided by the number of random barcodes it associated with, giving a 'duplication rate' (not great axis labels I agree, but this is pulling the plots straight out of a lab meeting I gave three years ago). The two plots show the same data, with a shortened X axis on the right to show the bulk of the spread better.

We can see that above a given frequency – in this case about 500, although it varies – we observe a 'duplication rate' around 70. This means that above a certain size, sequences are generally amplified at a rate proportional to their prevalence (give or take the odd outlier), or that for every input molecule of cDNA it gets amplified and observed seventy times. This is the scenario we'd generally imagine for PCR. However, below that variable threshold there is a very different, very noisy picture, where the amount to which a sequence is found to be amplified and observed is not related to the collapsed prevalence. This was the bait on the hook that lead our lab down this path.

Now, everyone knows PCR doesn't behave like it does in the diagrams, like it should. That's what everyone always says (usually as they stick another gel picture containing mysteriously sized bands into their labbooks). However, people have rarely looked at what's actually going on. There's a bit of special PCR magic that goes on, and a million different target and reaction variables that might affect things: you just optimise your reaction until your output looks like what you'd expect. It's only with the relatively recent advances in DNA sequencing technology that we can actually look at exactly what molecules are being made in the reaction that we can start to get actual data showing how just un-like the schematics the reaction can in fact behave.

This is exactly what Katharine's paper chases up, applying the same unique molecular barcoding strategy to TCR sequences from both polyclonal and monoclonal** T-cells. I won't go into the details, because hey, you can just read the paper (which says it much better), but the important thing is that this variability is not due to the standard things you might expect like CG content, DNA motifs or amplicon length, because it happens even for identical sequences. It doesn't matter how well you control your reactions, the noise in the system breeds variability. This makes unique molecular barcoding hugely important, at least if you want accurate relative quantitation of DNA species across the entire dynamic range of your data.

* Theoretically about 16.7 million in our case, or 4¹², as we use twelve random nucleotides in our barcodes.

** Although it's worth saying that while the line used, KT-2, is monoclonal, that doesn't mean the TCR repertoire is exactly as clean as you'd expect. T-cell receptor expression in T-cell lines is another thing that isn't simple as the textbook diagrams pretend.