Nick Crawford | Evolution and more...

Designing qPCR primers from just a GTF/GFF file and a genome sequence

By Nick On April 3, 2013 · Add Comment · In Blog

I recently had to design qPCR primers for some genes. I had a genome and an annotated GTF file derived from Cufflinks. Since I wanted the primers to span introns, to prevent the amplification of genomic DNA, I needed a both fasta file of coding sequence to use as input in to Primer3 as well as some associated information about where the introns were spliced out so I could ensure the primers I design spanned introns. So, what I did was use the cufflinks gffread command to convert the GTF file to create set of fasta transcripts annotated with the positions of each of the exons.

Here’s how to do it yourself. If you don’t have cufflinks installed you can download precompiled binaries from here, or if you’re using OS X you can use homebrew + the science tap to install it.

Once you have cufflinks installed type gffread –h to make sure everything is copacetic. Assuming that prints a bunch documentation to the terminal you can then run:

grep 'GeneID' my.gtf | gffread –g my.fasta -W -x GeneID.fasta

You’ll need to know how your gene or transcript is labeled in the GTF file (e.g. ‘GeneID’). The ‘-x’ flag ensures that the results only include coding sequence, the ‘-W’ adds the exon coordinates to the fasta header. Then run head GeneID.fasta to check that sequences were added to the output file.

One could probably further automate the next steps, but basically you just need to copy the first couple of exons/segs from the file and add brackets around the exon splice site you want the primers to span. You can use the ‘segs’ in the fasta header to determine where the splice junctions are.

>ENS0123412 gene=bactin loc:3(-)22053-225439 segs:1-333,334-490,491-518
ATGGCTCAGAGAGATGCTGACAAATACCTCTATGTGGATAGAAATCTCATCAACAACCCTCTTGCTCAGG
CCGATTGGGCAGCTAAGAAACTGGTGTGGGTCCCATCAGAAAAGAATGGCTTTGAGCCTGCTAGCTTAAA
AGAGGAAGTAGGAGATGAAGCCATTGTGGAGCTTGCAGAGAACGGGAAGAAAGTGCGAGTAAACAAAGAT
GATATCCAAAAGATGAACCCGCCTAAGTTCTCTAAAGTGGAAGACATGGCTGAATTGACCTGCCTGAATG
AGGCCTCTGTGTTGCACAACTTAAAGGAACGATACTACTCGGGGCTT[ATCTATACCTACT]CAGGCCTA
CTGTGTGGTCATAAATCCCTACAAGAACTTGCCCATCTACTCAGAAGAGATTGTGGAAATGTATAAGGGC
AAAAAGAGACACGAGATGCCCCCTCACATCTATGCCATTACAGACACAGCCTACAGGAGTATGATGCAAG

Lastly, use the modified fasta as input into Primer3. You’ll need to adjust the PCR product size to range from 70-200 base-pairs, but the default TM of 60 should be fine.

The final result should look something like this:

PRIMER PICKING RESULTS FOR ENS0123412 gene=bactin loc:3(-)22053-225439 segs:1-333,334-490,491-518

No mispriming library specified
Using 1-based sequence positions
OLIGO            start  len      tm     gc%   any    3' seq 
LEFT PRIMER        275   20   59.83   50.00  8.00  0.00 TGAATGAGGCCTCTGTGTTG
RIGHT PRIMER       450   20   60.03   55.00  5.00  3.00 TAGATGTGAGGGGGCATCTC
SEQUENCE SIZE: 488
INCLUDED REGION SIZE: 488

PRODUCT SIZE: 176, PAIR ANY COMPL: 4.00, PAIR 3' COMPL: 2.00
TARGETS (start, len)*: 328,13

    1 ATGGCTCAGAGAGATGCTGACAAATACCTCTATGTGGATAGAAATCTCATCAACAACCCT
                                                                  

   61 CTTGCTCAGGCCGATTGGGCAGCTAAGAAACTGGTGTGGGTCCCATCAGAAAAGAATGGC
                                                                  

  121 TTTGAGCCTGCTAGCTTAAAAGAGGAAGTAGGAGATGAAGCCATTGTGGAGCTTGCAGAG
                                                                  

  181 AACGGGAAGAAAGTGCGAGTAAACAAAGATGATATCCAAAAGATGAACCCGCCTAAGTTC
                                                                  

  241 TCTAAAGTGGAAGACATGGCTGAATTGACCTGCCTGAATGAGGCCTCTGTGTTGCACAAC
                                        >>>>>>>>>>>>>>>>>>>>      

  301 TTAAAGGAACGATACTACTCGGGGCTTATCTATACCTACTCAGGCCTACTGTGTGGTCAT
                                 *************                    

  361 AAATCCCTACAAGAACTTGCCCATCTACTCAGAAGAGATTGTGGAAATGTATAAGGGCAA
                                                                  

  421 AAAGAGACACGAGATGCCCCCTCACATCTATGCCATTACAGACACAGCCTACAGGAGTAT
                <<<<<<<<<<<<<<<<<<<<                              

  481 GATGCAAG

Filtering contigs/chromosomes from a multi-fasta file

By Nick On July 31, 2012 · 2 Comments · In Blog, Software

A colleague needed to remove some individual fastas from a multi-fasta file. Googling didn’t reveal a canned way to do it so I hacked up this script.

8.29.12 – As Jason Gallant pointed out, if your fasta is very small you don’t need to index your fasta file. Just use the simple biopython code he mentions in the comments.

#!/usr/bin/env python # encoding: utf-8 """ removeFaFromFasta.py Created by Nick Crawford on 2012-31-07. The author may be contacted at ngcrawford@gmail.com Requires samtools and pysam to be installed. http://samtools.sourceforge.net/ http://code.google.com/p/pysam/ """ import os import pysam import argparse def get_args(): """Parse sys.argv""" parser = argparse.ArgumentParser() parser.add_argument('-fa','--fasta', required=True, help='Path to fasta file to filter.') parser.add_argument('-o','--output', required=True, help='Path to filtered fasta file.') parser.add_argument('-l', '--filter-list', required=True, nargs='+', help='List of of chrm/contig names to remove.') args = parser.parse_args() return args def splitIterator(text, size): "Iterator that splits string into list of substrings." assert size > 0, "size should be > 0" for start in xrange(0, len(text), size): yield text[start:start + size] def main(args): # Setup files fasta = args.fasta filtered_fa = args.output faidx = args.fasta + '.fai' # Make sure .fai exists try: os.path.exists(faidx) except: print "You need to index the fasta file with samtools.\n {} does not exist.".format(faidx) # Filter names for contings bad_names = args.filter_list chrm_names = (line.strip().split()[0] for line in open(faidx,'rU')) filtered_chrm_names = (cn for cn in chrm_names if cn not in bad_names) # Write contigs/chrms to output fasta fasta = pysam.Fastafile(fasta) filtered_fa = open(filtered_fa,'w') for name in filtered_chrm_names: print "Processed:", name chrm = fasta.fetch(name) filtered_fa.write(">" + name + "\n") # split lines on 80 characters [filtered_fa.write(chars + "\n" ) for chars in splitIterator(chrm, 80)] # Clean up open files filtered_fa.close() if __name__ == '__main__': args = get_args() main(args)

Python: Adding Read Group (@RG) tags to BAM or SAM files

By Nick On April 17, 2012 · 1 Comment · In Blog, Software

The SAM specification now requires @RG tags to be included in all SAM/BAM alignments. If you are using GATK you have probably noticed that it will not run without them. Since @RG tags weren’t standard until recently, I’ve written a script to add them in post hoc. You’ll need to install pysam and python2.7 to get it to work.

import os import sys import glob import pysam import argparse import multiprocessing def get_args(): '''Parse sys.argv''' parser = argparse.ArgumentParser() parser.add_argument('--cores', type=int, help='the number of avalible processor cores') parser.add_argument('-i','--input-dir', required=True, help='The input directory containing the bam files.') parser.add_argument('-CN',type=str, help="Name of sequencing center producing the read. \ GATK not required.") parser.add_argument('-DS',type=str, help="Description. GATK Not Required. ") parser.add_argument('-DT',type=str, required=True, help="Date the run was produced (ISO8601 date or date/time). \ GATK Not Required. ") parser.add_argument('-PI',type=int, required=True, help="Predicted median insert size. GATK Not Required.") parser.add_argument('-PL',type=str, required=True, choices = ['CAPILLARY', 'LS454', 'ILLUMINA', 'SOLID', 'HELICOS', 'IONTORRENT', 'PACBIO'], help="Platform/technology used to produce the reads.") args = parser.parse_args() return args def addRG2Header(filename, sam_info, args): """Add read group info to a header.""" # CREATE TEMPLATE # Read group. Unordered multiple @RG lines are allowed. RG_template = { 'ID': '', # Read group identifier. e.g., Illumina flowcell + lane name and number 'CN': '', # GATK Not Required. Name of sequencing center producing the read. 'DS': '', # GATK Not Required. Description 'DT': '', # GATK Not Required. Date the run was produced (ISO8601 date YYYY-MM-DD or YYYYMMDD) 'PI': '', # GATK Not Required. Predicted median insert size. 'PU': '', # GATK Not Required. Platform unit (e.g. flowcell-barcode.lane for Illumina or slide for SOLiD). 'SM': '', # Sample. Use pool name where a pool is being sequenced. 'PL': 'ILLUMINA'} # Platform/technology used to produce the reads. samfile = pysam.Samfile(filename, 'r') new_header = samfile.header.copy() samfile.close() # ADD INFO TO TEMPLATE RG_template = RG_template.copy() RG_template['ID'] = sam_info["sample_name"] if args.CN: RG_template['CN'] = args.CN.upper() if args.DS: RG_template['DS'] = args.DS RG_template['DT'] = args.DT RG_template['LB'] = sam_info["sample_name"] RG_template['SM'] = sam_info["sample_name"] RG_template['DS'] = "{0}.{1}".format(sam_info['sample_name'], sam_info['locality']) RG_template['PI'] = args.PI RG_template['PU'] = '{0}.{1}'.format(sam_info['flowcell_id'], sam_info['lane']) new_header['RG'] = [RG_template] return new_header def add_RGs_2_BAMs_runner(data): """Generates the correct @RG header and adds a RG field to a bam file.""" # Make Bam of Sam if it doesn't exist filename, new_RG_header = data if filename.endswith('sam'): sam_handle = pysam.Samfile(filename) bam_name = os.path.splitext(filename)[0]+".bam" bam_handle = pysam.Samfile( bam_name, "wb", template = sam_handle ) for s in sam_handle: bam_handle.write(s) filename = bam_name # Massage paths and make outputfiles pysam.sort(filename, os.path.splitext(filename)[0]+".sorted") pysam.index(os.path.splitext(filename)[0]+".sorted.bam") filename = os.path.splitext(filename)[0]+".sorted.bam" path, filename = os.path.split(filename) name, ext = os.path.splitext(filename) new_name = name + '.wRG.' + 'bam' outfile_name = os.path.join(path,new_name) outfile = pysam.Samfile( outfile_name, 'wb', header = new_RG_header ) # Step 2: Process Samfile adding Read Group to Each Read samfile = pysam.Samfile(os.path.join(path, filename)) samfile.fetch() for count, read in enumerate(samfile.fetch()): name = read.qname read_group = os.path.split(filename)[1].split(".")[0] new_tags = read.tags new_tags.append(('RG', read_group)) read.tags = new_tags outfile.write(read) outfile.close() # Step 3: Make index of read group enabled samfile pysam.index(outfile_name) sys.stdout.write(".") sys.stdout.flush() return def parseFileName(filepath): path, filename = os.path.split(filepath) sample_name, locality, inline_tag, third_read_tag = filename.split(".")[:-1] sam_handle = pysam.Samfile(filepath,'r') sam_line = sam_handle.next() read_info = sam_line.qname sam_handle.close() instrument, run_id, flowcell_id, lane = read_info.split(":")[:4] info = {'sample_name': sample_name, 'locality': locality, 'inline_tag': inline_tag, 'third_read_tag': third_read_tag, 'instrument':instrument, 'run_id': run_id, 'flowcell_id': flowcell_id, 'lane': lane} return info def add_RGs_2_BAMs(pool, args): data_for_map = [] print 'Making RG headers.' for count, filename in enumerate(glob.glob(os.path.join(args.input_dir,'*'))): if filename.endswith('sam') or filename.endswith('bam'): sam_info = parseFileName(filename) new_RG_header = addRG2Header(filename, sam_info, args) data_for_map.append([filename, new_RG_header]) sys.stdout.write("\nAdding RGs and making BAMs") sys.stdout.flush() pool.map(add_RGs_2_BAMs_runner, data_for_map) return def main(): args = get_args() cores = args.cores pool = multiprocessing.Pool(args.cores) add_RGs_2_BAMs(pool, args) if __name__ == '__main__': main()

Python: Multiprocessing large files

By Nick On March 29, 2012 · Add Comment · In Blog

I been working with a lot of very large files and it has become increasing obvious that using a single processor core is a major bottleneck to getting my data processed in a timely fashion. A MapReduce style algorithm seemed like the way to go, but I had a hard time finding a useful example. After a bit of hacking about I came up with the following code.

The basic algorithmic idea is to first read in a large chunk of lines from the file. These are then partitioned out to the available cores and processed independently. The new set of lines are then written to an output file or in this example just printed to the screen. Normally this would be tricky code to write, but python 2.7′s wonderful multiprocessing module handles all the synchronization for you.

#!/usr/bin/env python # encoding: utf-8 import multiprocessing from textwrap import dedent from itertools import izip_longest def process_chunk(d): """Replace this with your own function that processes data one line at a time""" d = d.strip() + ' processed' return d def grouper(n, iterable, padvalue=None): """grouper(3, 'abcdefg', 'x') --> ('a','b','c'), ('d','e','f'), ('g','x','x')""" return izip_longest(*[iter(iterable)]*n, fillvalue=padvalue) if __name__ == '__main__': # test data test_data = """\ 1 some test garbage 2 some test garbage 3 some test garbage 4 some test garbage 5 some test garbage 6 some test garbage 7 some test garbage 8 some test garbage 9 some test garbage 10 some test garbage 11 some test garbage 12 some test garbage 13 some test garbage 14 some test garbage 15 some test garbage 16 some test garbage 17 some test garbage 18 some test garbage 19 some test garbage 20 some test garbage""" test_data = dedent(test_data) test_data = test_data.split("\n") # Create pool (p) p = multiprocessing.Pool(4) # Use 'grouper' to split test data into # groups you can process without using a # ton of RAM. You'll probably want to # increase the chunk size considerably # to something like 1000 lines per core. # The idea is that you replace 'test_data' # with a file-handle # e.g., testdata = open(file.txt,'rU') # And, you'd write to a file instead of # printing to the stout for chunk in grouper(10, test_data): results = p.map(process_chunk, chunk) for r in results: print r # replace with outfile.write()

Python: Interleave Paired-End Reads

By Nick On March 28, 2012 · Add Comment · In Blog

Here’s a simple script for interleaving paired-end fastq files. You’ll need to do this if you want to create input files for velvet. Unlike the velvet’s shuffleSequences_fastq.pl perl script, this script handles gzipped input and output. It requires python 2.7.

#!/usr/bin/env python # encoding: utf-8 import gzip import argparse def interface(): args = argparse.ArgumentParser() args.add_argument('-l', '--left-input', help='The first input file in fastq format.') args.add_argument('-r', '--right-input', help='The first input file in fastq format.') args.add_argument('-o', '--output', help='The output file in fastq format.') args = args.parse_args() return args def process_reads(args): if args.left_input.endswith('.gz') or args.right_input.endswith('.gz'): left = gzip.open(args.left_input,'rb') right = gzip.open(args.right_input,'rb') fout = gzip.open(args.output,'wb') else: left = open(args.left_input,'rU') right = open(args.right_input,'rU') fout = open(args.output,'wb') # USING A WHILE LOOP MAKE THIS SUPER FAST # Details here: # http://effbot.org/zone/readline-performance.htm while 1: # process the first file left_line = left.readline() if not left_line: break fout.write(left_line) left_line = left.readline() fout.write(left_line) left_line = left.readline() fout.write(left_line) left_line = left.readline() fout.write(left_line) # process the second file right_line = right.readline() fout.write(right_line) right_line = right.readline() fout.write(right_line) right_line = right.readline() fout.write(right_line) right_line = right.readline() fout.write(right_line) left.close() right.close() fout.close() return 0 if __name__ == '__main__': args = interface() process_reads(args)

Bowtie2 output as BAM

By Nick On March 14, 2012 · 3 Comments · In Software

Bowtie2 is a short read aligner that is optimized for aligning longer reads of lengths of 50 bp or greater. I’ve been playing around with it and was initially puzzled by the fact that it only outputs SAM formated alignments. Then I realized you can pipe the output straight into samtools which will do the compression to BAM for you.

$ bowtie2 \
-p 4 \
-x /genome/index \
-1 pair2.fastq \
-2 pair2.fastq \
-U unpaired.fastq \
--very-sensitive \
-X 1000 \
-I 200 \
| samtools view -bS - > output.bam

Installing ABySS

By Nick On December 5, 2010 · Add Comment · In Blog

To install ABySS on an system running an older version of gcc and use the following commands.

>./configure –enable-maxk=96 –disable-openmp \

CPPFLAGS=-I<path to google-sparsehash install>/include \

–prefix=<ABySS install directory>/

>make AM_CXXFLAGS=-Wall

>make install

Word by word in Terminal

By Nick On September 27, 2010 · Add Comment · In Blog

One of the annoying things in OS terminal is that if you want to traverse word by word in a line of text you need to type ‘esc-b’ and ‘exc-f’. This post on macromates, the textmate blog, explains how to reset these keys. Enjoy.

MultiMarkdown

By Nick On June 9, 2010 · Add Comment · In Blog, Software

Since I started using github in a serious way back in January I’ve begun writing my documentation in the markdown format that displays so nicely on github. Markdown is essentially a parsing tool and a simple text syntax that allows the easy conversion of human ‘readable text’ to html. It’s intuitive, it took less than 5 minutes to pick up, and saves me a ton of time not writing HTML. However, its ease of use is tempered, a bit, by a lack of features. Although it is easy to create headers, lists, and code bocks – simple HTML stuff – it doesn’t include the option to create tables, formated mathematical formulas, citations and bibliographies. Since I’m a scientist who wants to produce documents with these sorts of features, this is annoying.

Luckily, the markdown syntax has recently been extended, in a project called MultiMarkdown, to include many of the aforementioned features. Multimarkdown essentially merges the markdown syntax with LaTeX which, if you haven’t heard of it, is a rather inscrutable, but extremely powerful text formatting language. It’s popular in the CS and physics disciplines. LaTeX produces beautiful documents, but it’s easy to spend a week or more adjusting the formatting and reading the API trying to figure out some of the more complicated features. Multimarkdown looks like it will do much of the more basic LaTeX formatting, but without the headache.

Vertebrate Zoology: Bi302

By Nick On January 22, 2010 · Add Comment · In Courses

Welcome to Vertebrate Zoology lab (Bi302 Lab). I’ll be using this portion of my website to post notes/slides and to answer your questions. Please make use of the comments. More to come.