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.

Extracting sample names from Illumina FASTQ filenames

Another post in the vein of "minor analysis task and my solution": how to extract the sample name or identifier (i.e. what you called a given file in the samplesheet) from a list of Illumina FASTQ file names.

So you might have something that looks like this (which I've cobbled together from a couple of runs from different machines that I pulled off Basespace):
 
322_S18_R2_001.fastq.gz 
1c8_P6_act_S11_R2_001.fastq.gz
5c4_P8_res_S1_R1_001.fastq.gz

In this exercise, I wanted the first part of the IDs, which are what we labelled the actual samples with. In an ideal world this should be very easy: we could just use sed or cut to delete everything after the first instance of a character used in the text that Illumina appends to the sample name. However different people (or the same people at different times) use different naming conventions, which can often include underscores, the character that separates the string fields that get added.

Therefore what we need to do is not delete from say the first pattern, but from the nth pattern from the end, which turned out to be a slightly trickier problem (at least to one such as myself that doesn't find sed to be the most intuitive program in the world). Here's what I came up with:

sed 's/\(.*\)_\(.*_\)\(.*\)_/\1|\2/;s/|.*//' FilenameList.txt

What this does is find the third underscore to last (just before S[XX]_R[X].fastq.gz string) and replace it with a pipe character ('|'), which should never appear in a filename, before finding this pipe character and deleting everything after it. This produces the following output:
 
322
1c8_P6_act
5c4_P8_res

This feels a bit hacky to me, as really I feel like I should be able to do this in one substitution, but the first part of the command actually selects the latter half (so if anyone knows how to tidy this up I'd love to hear it!).

It's also worth pointing out that the output of different machines or demultiplexing settings will differ in the number of underscores, and so the command will need to change slightly. For instance, the above example was run on NextSeq data demultiplexed using bcl2fastq --no-lane-splitting. For an example of a file containing lane information ('L001' etc), here's the equivalent on some MiSeq filenames I pulled up.

# This 
324bAB_S8_L001_R1_001.fastq.gz 
076-PBMCs-alpha_S5_L001_R1_001.fastq.gz 

# After this
sed 's/\(.*\)_\(.*_\)\(.*\)_\(.*\)_/\1|\2/;s/|.*//' FilenameList.txt   

# Becomes this 
324bAB
076-PBMCs-alpha

The reason I was actually doing this was because I had a number of different NextSeq runs, some of which contained samples sequenced from the same original libraries (to increase the depth), so I needed to get the names of all samples sequenced in each run. I had each run in separate directories, and then ran this in the folder above:

for i in */
do echo $i
ls $i | grep R1 | grep fastq.gz | sed 's/\(.*\)_\(.*_\)\(.*\)_/\1|\
2/;s/|.*//' > ${i%?}.fl
done

This goes through all the directories, finds the R1 FASTQ files (as there are other files in there, and this gives me only one entry per sample) and strips their sample names out as above, and writes these to a textfile with the suffix '.fl'.

Retreiving GenBank FASTA records with Python

I recently needed to download a relatively large number of GenBank FASTA records from GenBank (as I wanted to have a play with the paired alpha-beta TCR chains sequenced from single cells in this lovely paper.

Unfortunately what is supposedly probably the easiest way - using NCBI's Batch Entrez function - did not work for me. Instead I just got lots of "the following records can't be retrieved/might be obsolete" and "No items found" warnings.

However there are a number of other ways, but I especially liked Pierre Lindenbaum's answer on https://www.biostars.org/p/3297/, which basically uses bash to pull the record directly from the NCBI website. Nice and simple, and it also gives the URL format required.

As I want to do some post-download processing in Python, I thought I may as well download the records directly using the urllib2 module. In case it may be of use to others (and for the next time I need to do this myself!) I thought I'd throw the basic code up here. The only requirement that this code has is to already have a file containing a list of GenBank accession IDs: in this case I just generated this in Excel for the ranges needed, but it could be easily coded directly in Python.

import urllib2

with open('GenbankAccessions') as infile, open('Sequences.fasta', "w") as outfile:
  for accession in infile:
    url = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=" + accession.rstrip() + "&rettype=fasta"
    entry = urllib2.urlopen(url)
    for line in entry:
      outfile.write(line)

This will write all of the FASTA records of the accessions listed in the 'GenbankAccessions' file to the new out file 'Sequences.fasta'.

The key to finding TCR sequences in RNA-seq data

I had previously written a short blog post touching on how I'd tried to mine some PBMC RNA-seq data (from the ENCODE project) for rearranged T-cell receptor genes, to try and open up this huge resource for TCR repertoire analysis. However, I hadn't gotten very far, on account of finding very few TCR sequences per file.

That sets the background for an extremely pleasant surprise this morning, when I found that Scott Brown, Lisa Raeburn and Robert Holt from Vancouver (the latter of whom being notable for producing one of the very earliest high-throughput sequencing TCR repertoire papers) had published a very nice paper doing just that!

This is a lovely example of different groups seeing the same problem and coming up with different takes. I saw an extremely low rate of return when TCR-mining in RNA-seq data from heterogeneous cell types, and gave up on it as a search for needles in a haystack. The Holt group saw the same problem, and simply searched more haystacks!

This paper tidily exemplifies the re-purposing of biological datasets to allow us to ask new biological questions (something that I consider a practical and moral necessity, given the complexity of such data and the time and costs involved in their generation).

Moreover, they do some really nice tricks, like estimating TCR transcript proportions in other data sets based on constant region usage, investigate TCR diversity relative to CD3 expression, testing on simulated RNA-seq data sets as a control, looked for public or known-specificity receptors and inferred possible alpha-beta pairs by checking all each sample's possible combinations for their presence in at least one other sample (somewhat akin to Harlan Robins' pairSEQ approach).

All in all, a very nice paper indeed, and I hope we see more of this kind of data re-purposing in the field at large. Such approaches could certainly be adapted for immunoglobulin genes. I also wonder if, given whole-genome sequencing data from mixed blood cell populations, we might even be able to do a similar analysis on rearranged variable antigen receptors from gDNA.

TCR trivia part one - translating T cell receptor sequences

Here's the first in what's likely to a be pretty niche series of posts – but hopefully useful to some – on various bits of TCR trivia I've gleaned during my PhD. This installment: translating TCR sequences.

Whether you're doing high- or low-throughput DNA TCR sequencing, it's highly likely that at some point you'll want to translate your DNA into amino acid sequences.

There's usually a few main goals to bear in mind when translating; to check the frame, and (potential) functionality of the sequence (i.e. lacking premature stop codons), before defining the hypervariable CDR3 region.

Assuming your reads are indel free, sequencing variable antigen receptors will always involve a bit more thought about checking the reading frame than regular amplicon sequencing, due to the non-templated addition and deletion of nucleotides that occurs during V(D)J recombination.

So how do can we tell whether our final recombined sequence is in frame or not?

This will mostly depend on how your sequencing protocols work.

If, say you've somehow managed to sequence an entire TCR mRNA sequence, then it's easy; you can just take the sequence from the start codon of the leader region to the stop codon of the constant region and if it divides exactly by three then it's in frame. Assuming there's no stop codons in between and the CDR3 checks out, chances are good your sequence encodes a productively rearranged TCR chain.

Given today's technology, that's unlikely to be the case; whether you're cloning into plasmids and Sanger sequencing or throwing libraries into some next-gen machine, chances are good that the sequence you'll be working with is actually just a small window around the recombination site.

What's likely is that you have amplified from a constant or J region primer at one end, to a V or RACE/template-switching primer at the other. Sorting out the frame is then just a matter of finding sequences on either side of the rearrangement that you know should be in frame, and seeing whether they are.

So far so self-explanatory. Here's the first of the fiddlier bits of TCR trivia I want to impart, the thing I need to remind myself of every time I tweak my TCR translating scripts: in fully-spliced TCR message the last nucleotide of the J region makes up the first nucleotide of the first codon of the constant region.

This means that if you go from the first nucleotide of the V (or the leader sequence, all functional leader sequences being divisible by three) to either the second to last nucleotide of the J, or the second nucleotide of the C (if you've started with mRNA instead of gDNA) then presto, you'll be in frame!

This feature also produces another noteworthy feature, at least in one chain; while every TRBJ ends in the same nucleotide, there are functional TRAJ genes ending in each base. This means that the first residue of the constant region of the alpha chain can be one of four different amino acids – which can be a bit off-putting if you don't know this and you're looking at what's supposed to be constant.

Once you've translated your TCR the next step is to find the CDR3, which should surely be the easiest bit; just run from the second conserved cysteine to the phenylalanine in the FGXG motif, right, just as IMGT says? Simple, or so it might seem.

Both sides of the CDR3 offer their own complications. First off, finding the second conserved cysteine isn't as easy as just picking the second C from the 5', or the 5'-most one upstream of the rearrangement, as some Vs have more than two germline cysteines, and it's feasible that new ones might be generated in the rearrangement. At some point you're just going to have to do an alignment of all the different Vs, and record the position of all the conserved cysteines.

The FGXG is also trickier than it might seem, by merit of the fact that there are J genes in both the alpha and gamma loci* that (while supposedly functional) lack the motif, being either FXXG or XGXG. If you're only looking for CDR3s matching the C to FGXG pattern, then there's going to be whole J genes which never seem to get used by your analysis!

There you have it, a couple of little tips that are worth bearing in mind when translating TCR sequence**.

Addendum - it's actually a little bit more complicated. Naturally.

* If you're interested (or don't believe me), they are: TRAJ16, TRAJ33, TRAJ38, TRGJP1 and TRGJP2.

** Note that everything written here is based on human TCR sequence as that's what I work with, but most of it will probably apply to other species as well.

Biopython SeqIO tips

The Biopython libraries are of tremendous use in analysing biological data. I use them in the vast majority of my bespoke fastq analysis, as I'm sure many do.

However, there's a couple of tasks I regularly find myself wanting to do, but that I could not easily find solutions for. In case anyone finds themselves in similar need, here are the solutions I found, maybe save you a bit of time.

Output quality scores

This is the main one; it's always annoyed me that there doesn't seem to be an easy way to output the quality score in SeqIO. Sure, I know with the letter_annotations option I can output the actual score in numbers, but sometimes I want to output the actual ASCII characters (such as if I want to take a subsection of a fastq record, both of nucleotide and quality strings).

Here's how I get around it; turn the whole record into a string, split that on its new lines, and just take the fourth:

 str(record.format("fastq")).split("\n")[3]  

Generate new fastq records in situ

Sometimes it's necessary to generate complete fastq records on the fly, as opposed to reading them in from an existing file.

I've found a couple of ways of doing this. The first comes buried in amongst some unrelated Biopython features, and looks like this:

 from Bio.Alphabet import generic_dna  
 from Bio.Seq import Seq  
 from Bio.SeqRecord import SeqRecord  
 new_record = SeqRecord(Seq("AAAAAA", generic_dna), id="New Record", description="")  
 new_record.letter_annotations["phred_quality"] = [40,40,40,40,40,40]   

The second way came from an answer on Biostars, and works like this:

 from Bio import SeqIO  
 from StringIO import StringIO  
 fastq_string = "@%s\n%s\n+\n%s\n" % ("New Record", "AAAAAA", "IIIIII")  
 new_record = SeqIO.read(StringIO(fastq_string), "fastq-illumina")  

They both work, with a couple of differences.

According to a quick test, the latter is faster, and obviously requires loading fewer modules.

Which to use might also depend on what format you have your qualities in; if you only have integers, then the first might be more tempting, whereas if you're making your new fastq from pieces of other existing records then the second is probably the way to go.

Generating random amino acid/DNA/RNA sequences in Python

This is my first post in what I hope to be a recurrent theme, that of bioinformatic practice.

I am by training a wet-lab scientist; my degrees were in molecular genetics, and I'm now doing my PhD in immunology. I've always had a keen interest in bioinformatics, and always tried to do as much as I can, but I've not been formally trained (and I've never claimed to be the most natural mathmatician/statistician).

Accordingly, whenever I'm trying to do some new bit of coding or analysis, I usually like to look around in manuals, papers, blogs and forums to see what and how people are doing things, and see if I can prevent some arduous wheel-reinventing.

I've also noticed that a lot of the time some little snippet of info or code that would have proved very useful to me is omitted by more practised programmers, to whom it would have been so obvious it goes unsaid or forgotten.

In that vein, I plan to post some of the bioinformatic tricks and codes I find and use as I go, in the hopes that I can pay it forward to others who find myself in my position, that is (usually), knowing exactly what I want to do, just not knowing off-hand the names of the functions to do it.

So, kicking off, here's something that I've been meaning to code for a little while; a way to generate random biological sequences (be they protein or nucleic acid sequences), written in Python, based on this biostars post.

This example produces random peptide sequences. For my current purposes, I also wanted these peptides to be of variable length, around a given normal distribution (that I'd measured empirically from the data set I was making these random sequences to act as a control for).

Note that I also wanted to plot the distribution of the newly created peptides; if you didn't want that functionality you can just delete every line with a "# plot" comment.

import random
import matplotlib.pyplot as plt # plot

def random_AA_seq(length):
    return ''.join(random.choice('ACDEFGHIKLMNPQRSTVWY') for i in range(length))

list_size = 3000

lengths = []

for j in range(list_size):

    a = int(random.normalvariate(17, 2))
    lengths.append(a)  # plot
    print random_AA_seq(a)

plt.hist(lengths, bins=14) # plot
plt.show() # plot

You can adjust the distribution of lengths by changing the values fed to normalvariate; the former is the mean, the latter the standard deviation (I deduced the values of 17 and 2 here from the sample I was trying to match).

Then of course these is easily adapted to DNA or RNA by swapping 'ACDEFGHIKLMNPQRSTVWY' for 'ACGT' or 'ACGU'.

Also, thanks go to David Craft, whose blog post on Syntax Highlighter came in useful in this post!