Welcome to Scrimer¶

Special cases¶

If you’re in an environment where you’re not able to install virtualenv systemwide, we recommend using the technique described at http://eli.thegreenplace.net/2013/04/20/bootstrapping-virtualenv/.

If you’re in a grid environment, this can help with paths that differ on different nodes

virtualnev --relocatable ~/scrimer-env

Non-python dependecies¶

Apart from the Python modules, the Scrimer pipeline relies on other tools that should be installed in your PATH. Follow the installation instructions in each package.

For reference we recorded the commands used to install those dependencies in the scrimer virtual box image. If your system is Debian 7, the commands could work verbatim.

• bedtools [8] is a dependency of pybedtools, used for manipulating with gff and bed files
• samtools [12] is used for manipulating short read alignments, and for calling variants
• LASTZ [7] is used to find the longest isotigs
• tabix [9] creates compressed and indexed verisions of annotation files
• GMAP [11] produces a spliced mapping of your contigs to the reference genome
• smalt [13] maps short reads to consensus contigs to discover variants
• GNU parallel [14] is used throughout the pipeline to speed up some lengthy calculations [34]
• blat and isPcr [15] are used to check the designed primers
• Primer3 [16] is used to find the most optimal primes sequences
• cutadapt [6] is used to remove cDNA synthesis primers.

Additional tools can be installed to provide some more options.

• FastQC [20] can be used to check the tag cleaning process
• agrep [21] and tre-agrep [22] can be used to check the tag cleaning process
• sort-alt [10] provides alphanumeric sorting of chromosome names, rename sort to sort-alt after compiling
• IGV [23] is great for visualizing the data when checking the results
• newbler [24] is the best option for assembling 454 mRNA data [32] [33]
• MIRA [25] does well with 454 transcriptome assembly as well [32] [33]
• sim4db [26] can be used as alternative spliced mapper, part of the kmer suite, apply our patch [27] to get standard conformant output
• Pipe Viewer [28] can be used to display the progress of longer operations
• BioPython [18] and NumPy [19] are required for running 5prime_stats.py
• mawk [29], awk is often used in the pipeline, and mawk is usually an order of magnitude faster
• vcflib [31] has a nice interface for working with vcf files (but new bcftools are good as well)

Add installed tools to your PATH¶

To easily manage locations of the tools that you’re using with the Scrimer pipeline, create a text file containing paths to directories, where binaries of your tools are located. The format is one path per line, for example:

/opt/bedtools/bin
/opt/samtools-0.1.18
/opt/lastz/bin
/opt/tabix
/opt/gmap/bin
/data/samba/liborm/sw_testbed/smalt-0.7.4
/data/samba/liborm/sw_testbed/FastQC
/data/samba/liborm/sw_testbed/kmer/sim4db

Put this file to your virtual environment directory, e.g. ~/scrimer-env/paths. You can run the following snippet when starting your work session:

export PATH=$( cat ~/scrimer-env/paths | tr "\n" ":" ):$PATH

References¶

Virtual machine installation¶

Installation steps (it should take less than 10 minutes):

• Install VirtualBox (https://www.virtualbox.org/wiki/Downloads). It works on Windows, Linux and Mac.
• Download the virtual machine image from http://goo.gl/Xf2cVU. You’ll get a single file with .ova extension on your hard drive.
• You can either double click the .ova file, or run VirtualBox, and choose File > Import Appliance. Follow the instructions after the import is started.

After a successful instalation you should see something like this (only the machine list will contain jsut one machine). Check whether you can start the virtual machine: click Start in the main VirtualBox window:

After a while you should see something like this:

You don’t need to type anything into that window, just checking that it looks like the screenshot is enough.

Machine configuration details:

• Administrative user: root, password: debian
• Normal user: user, password: user
• ssh on port 2222

Connecting to the virtual machine¶

Connect to control the machine¶

To control the machine, you need to connect to the ssh service.

In Windows this is done with PuTTY.

• start PuTTY
• fill Host Name: localhost
• fill Port: 2222
• click Open or press <Enter>

In the black wnidow that appears, type your credentials:

• login as: user
• user@localhost's password: user

In Mac OS X or Linux, this is done with ssh tool:

ssh -p 2222 user@localhost

Connect to copy files¶

In Windows, WinSCP is used to copy files to Linux machines. You use the same information as for PuTTY to log in.

In Mac OS X or Linux, the most simple command to copy a file into a home directory of user on a running virtual machine is:

scp -P 2222 myfile user@localhost:~

Directory layout¶

Here we present the layout that we use to organize all the data needed to run the pipeline. The inputs together with intermediate (and final) results total to hundreds of files. Having those files organized can help prevent mistakes.

project.sh¶

Create a file called project.sh in your project directory. It will consist of KEY=VALUE lines that will define your project specific settings, and each time you want to use Scrimer you’ll start by:

cd my/project/directory
. project.sh

Example project.sh :

# number of cores you want to use for parallel calculations
CPUS=8

# location of genome data in your system
# you need write access to add a new reference genome to that location
GENOMES=/data/genomes

# reference genome used
GENOME=taeGut1
GENOMEDIR=$GENOMES/$GENOME
GENOMEFA=$GENOMEDIR/$GENOME.fa

# genome in blat format
GENOME2BIT=$GENOMEDIR/$GENOME.2bit

# gmap index location
GMAP_IDX_DIR=$GENOMEDIR
GMAP_IDX=gmap_${GENOME}

# smalt index
SMALT_IDX=$GENOMEDIR/smalt/${GENOME}k13s4

# primers used to synthetize cDNA
# (sequences were found in .pdf report from the company that did the normalization)
PRIMER1=AAGCAGTGGTATCAACGCAGAGTTTTTGTTTTTTTCTTTTTTTTTTNN
PRIMER2=AAGCAGTGGTATCAACGCAGAGTACGCGGG
PRIMER3=AAGCAGTGGTATCAACGCAGAGT

Adding the tools to your PATH¶

For each scrimer session, all the executables that are used have to be in one of the directories mentioned in PATH. You can set up your PATH easily by using the file you created during installation:

export PATH=$( cat ~/scrimer-env/paths | tr "\n" ":" ):$PATH

Such line can be at the end of your project.sh file, so everythig is set up at once.

Alternatively you can copy all the tool executables into your virtual environment bin directory (~/scrimer-env/bin).

Download and prepare the reference genome¶

• A list of available genomes is at http://hgdownload.cse.ucsc.edu/downloads.html
• We download the full data set, but it’s possible to interrupt the download during xenoMrna (not needed, too big).
• md5sum is a basic utility that should be present in your system, otherwise check your packages (yum, apt-get, ...)

# location of genome data that can be shared among users
mkdir -p $GENOMEDIR
cd $GENOMEDIR

rsync -avzP rsync://hgdownload.cse.ucsc.edu/goldenPath/$GENOME/bigZips/ .

# check downloaded data integrity
md5sum -c md5sum.txt
cat *.md5 | md5sum -c

Now unpack the genome. This process differs for different genomes - some are in single .fa, some are split by chromosomes. Some archives are tarbombs, so unpack to directory chromFa to avoid a possible mess:

mkdir chromFa
tar xvzf chromFa.tar.gz -C chromFa

Create concatenated genome, use Heng Li’s sort-alt to get the common ordering of chromosomes:

find chromFa -type f | sort-alt -N | xargs cat > $GENOME.fa

Download all needed annotations¶

Annotation data is best obtained in UCSC table browser in BED format and then sorted and indexed by BEDtools

For example: http://genome.ucsc.edu/cgi-bin/hgTables?db=taeGut1:

# directory where annotations are stored
ANNOT=annot
sortBed -i $ANNOT/ensGene.bed > $ANNOT/ensGene.sorted.bed
bgzip $ANNOT/ensGene.sorted.bed
tabix -p bed $ANNOT/ensGene.sorted.bed.gz

FIXME: rozepsat Or using compressed files:

zcat -d $ANNOT/refSeqGenes.bed.gz | sortBed | bgzip > $ANNOT/refSeqGenes.sorted.bed.gz
zcat -d $ANNOT/ensGenes.bed.gz | sortBed | bgzip > $ANNOT/ensGenes.sorted.bed.gz
tabix -p bed $ANNOT/ensGenes.sorted.bed.gz
tabix -p bed $ANNOT/refSeqGenes.sorted.bed.gz

Build indexes for all programs used in the pipeline¶

Some programs need a preprocessed form of the genome, to speed up their operation.

# index chromosome positions in the genome file for samtools
samtools faidx $GENOMEFA

# build gmap index for zebra finch
# with some newer versions it is necessary to use -B <path/to/bindir>
# beware, this requires quite a lot of memory (gigabytes)
gmap_build -d $GMAP_IDX -D $GMAP_IDX_DIR $GENOMEFA

# smalt index
# needed only for speeding up sim4db
mkdir -p $GENOMEDIR/smalt
smalt index -s 4 $SMALT_IDX $GENOMEFA

# convert to blat format
faToTwoBit $GENOMEFA $GENOME2BIT

Quality check of the raw data¶

Results of the quality check of the raw data can be used as a reference point to check the improvements done by this step.

IN=00-raw
OUT=10-fastqc
mkdir $OUT
fastqc --outdir=$OUT --noextract --threads=$CPUS $IN/*.fastq

Split the files by multiplex tags¶

If you used multiplexing during library preparation, your reads have either already been split in your sequencing facility, or you have to split the reads according to the ‘barcodes’ stored in the sequences yourself.

• for 454, this is done by sfffile, you have to convert resulting sff files to fastq format
• for Illumina this can be done e.g. by eautils

IN=03-sff
OUT=04-sff-split
mkdir -p $OUT

parallel -j 1 sfffile -s RLMIDs -o $OUT/{/.} {} ::: $IN/*.sff

Remove cDNA synthesis primers with cutadapt¶

Another source of noise in the data is the primers that were used for reverse transcription of mRNA and for the following PCR amplification of cDNA. We remove them using cutadapt.

If you have more or less than three primers, you have to change the command by adding/removing --anywhere= parts.

# data from previous steps
IN=10-mid-split
OUT=12-cutadapt
mkdir $OUT

# cut out the evrogen sequences using GNU parallel and cutadapt
# cutadapt supports only 'N' wildcards, no ambiguity codes
parallel -j $CPUS "cutadapt --anywhere=$PRIMER1 --anywhere=$PRIMER2 --anywhere=$PRIMER3 \
  --error-rate=0.2 --overlap=15 --minimum-length=40 \
  --output=$OUT/{/.}.fastq --rest-file=$OUT/{/.}.rest {}" ::: $IN/*.fastq > $OUT/cutadapt.log

Check the results¶

It is necessary to check the results of adaptor cutting.

First you can check how many of the primers were missed by cutadapt. agrep uses a different matching algorithm than cutadapt, so some remaining hits are usually found. /dev/null is used as the second input to agrep so the filenames are output.

NERR=5
parallel agrep -c -$NERR "$PRIMER3" {} /dev/null ::: $OUT/*.fastq|grep -v /dev

The next thing to check is the logs produced by cutadapt. Results for our data - 454 Titanium data from Smart kit synthesized cDNA:

• ~70% trimmed reads
• ~10% trimmed basepairs
• ~10% too short reads

grep -A5 Processed $OUT/cutadapt.log | less

The mean length of the removed sequences should be close to length of the adapter (31 in this case):

less $OUT/cutadapt.log

# Lengths of removed sequences (5')
# length  count   expected
# 5       350     264.7
# 6       146     66.2
# ...
# 30      6414    0.0
# 31      63398   0.0
# 32      6656    0.0
# ...

The size of the .rest files is 1/500 of the .fastq (should be 1/250 for .fasta)

ls -l $OUT

The fastqc checks should be +- ok.

fastqc --outdir=13-fastqc --noextract --threads=8 $OUT/*.fastq

Visual debugging¶

If something in the previous checks looks weird, look directly at the data. Substitute filenames below with the names of your files.

Look where the primers are in the sequence. tre-agrep is used to color the output of agrep, because agrep throughput is ~ 42 MB/s while tre-agrep throughput is ~ 2 MB/s.

FQFILE=$IN/G3UKN3Q01.fasta
NERR=5
agrep -n -$NERR "$PRIMER3" $FQFILE |tre-agrep -$NERR "$PRIMER3" --color|less -S -R

To find out how many differences should be allowed in the pattern matching, we try to find a value of NERR where the primer sequence starts to match randomly inside the reads, and not only in at the beginning. Notice the ^ in the first command, marking the start of the read.

agrep -c -$NERR "^$PRIMER3" $FQFILE && agrep -c -$NERR "$PRIMER3" $FQFILE

# numbers for tag-cleaned G59B..
# 4 errors: 11971 12767
# 5 errors: 16366 17566
# 6 errors: 17146 23858
# 7 errors: 18041 67844

In our sample results, numbers start to diverge for NERR > 5, so 5 is a good choice.

Read count statistics¶

For a single file:

# read count statistics
# @ can be in the beginning of quality string, so filter the rows in order

# count of sequences
awk '((NR%4)  == 1)' $FQFILE | wc -l
# or more effective
echo $(( $(wc -l $FQFILE) / 4 ))

# count of sequenced bases
awk '((NR%4)  == 2)' $FQFILE | wc -m

For all files in OUT:

# parallel, IO bound task, so run one process a time
OUT=12-cutadapt
echo "read_count base_count filename"
parallel -j 1 'echo $( gawk "((NR%4)  == 1)" {} | wc -l ) $( gawk "((NR%4)  == 2)" {} | wc -m ) {}' ::: $OUT/*.fastq

Use newbler to assemble the reads¶

Here newbler is used to assemble the contigs. For 3 GB of read data the assembly took 25 CPU hours and 15 GB RAM.

IN=12-cutadapt
OUT=22-newbler
CPUS=19

# run full assembly
runAssembly -o $OUT -cdna -cpu $CPUS -m $IN/*.fastq

Remove contigs that are similar to each other¶

The aim is to get one transcript per locus, preferably the longest one. Otherwise the read mapping process would be faced with many ambiguous locations. We achieve this by:

• doing all-to-all comparison within the isotigs
• grouping isotigs that are similar up to given thresholds of coverage and identity, (disjoint sets forest graph algorithm is used)
• and selecting only the longest contig for each group

IN=22-newbler
OUT=$IN
SEQFILE=$IN/454Isotigs.fna
MINCOVERAGE=90
MINIDENTITY=95

# taken from lastz human-chimp example, should be report only highly similar hits
# filter self matches with awk
lastz $SEQFILE[multiple] $SEQFILE \
      --step=10 --seed=match12 --notransition --exact=20 --noytrim \
      --match=1,5 --ambiguous=n \
      --coverage=$MINCOVERAGE --identity=$MINIDENTITY \
      --format=general:name1,size1,start1,name2,size2,start2,strand2,identity,coverage \
      | awk '($1 != $4)' > $SEQFILE.lastz-self

# takes 8 minutes, finds 190K pairs

# find the redundant sequences
tail -n +2 $SEQFILE.lastz-self | find_redundant_sequences.py > $SEQFILE.redundant

# add the short sequences to discard list
grep '>' $SEQFILE | awk '{ sub(/length=/,"",$3); sub(/^>/, "", $1); if($3 - 0 < 300) print $1;}' >> $SEQFILE.redundant

# get rid of the redundant ones
seq_filter_by_id.py $SEQFILE.redundant 1 $SEQFILE fasta - $SEQFILE.filtered

Check the results¶

# check the input contig length distribution
grep '>' $SEQFILE | awk '{ sub(/length=/,"",$3); print $3}' | histogram.py -b 30

# view how many contigs we got after filtering
grep -c '>' $SEQFILE.filtered

Assembling Illumina data¶

Trinity gives fairly good results for transcriptome assembly from Illumina data. Simple Trinity usage looks like:

# as mentioned on the homepage
Trinity --seqType fq --JM 50G --left reads_1.fq  --right reads_2.fq --CPU 6

Some more tips on assembling ‘perfect’ transcripts by Don Gilbert.

Mapping options¶

INFILE=20-jp-contigs/lu_master500_v2.fna.filtered
OUT=30-tg-gmap
mkdir $OUT

OUTFILE=$OUT/$( basename ${INFILE%%.*} ).gmap.gff3

gmap -D $GMAP_IDX_DIR -d $GMAP_IDX -f gff3_gene -B 3 -x 30 -t $CPUS\
    --cross-species $INFILE  > $OUTFILE

INFILE=20-newbler/454Isotigs.fna.filtered
OUT=31-tg-sim4db
mkdir $OUT

# these values are derived, it's not necessary to change them
FRAGS=$OUT/${INFILE##*/}.frags
SMALT_OUT=$FRAGS.cigar
SIM4_SCR=${FRAGS%.*}.sim4scr
OUT0=${FRAGS%.*}.tmp.gff3
OUTFILE=${FRAGS%.*}.gff3

# create fragments, using slightly modified fasta_fragments.py from lastz distribution
cat $INFILE | fasta_fragments.py --step=80 > $FRAGS

# map the fragments with smalt (takes few minutes), reporting all hits (-d -1) scoring over 60
smalt_x86_64 map -n 8 -f cigar -o $SMALT_OUT -d -1 -m 60 $SMALT_IDX $FRAGS

# construct the script for sim4db
cat $SMALT_OUT | cigar_to_sim4db_scr.py $GENOMEFA.fai | sort --key=5n,5 > $SIM4_SCR

sim4db -genomic $GENOMEFA -cdna $INFILE -script $SIM4_SCR -output $OUT0 -gff3 -interspecies -mincoverage 70 -minidentity 90 -minlength 60 -alignments -threads $CPUS

# fix chromosome names
sed s/^[0-9][0-9]*:chr/chr/ $OUT0 > $OUTFILE

Transfer genome annotations to our contigs¶

Annotate our sequences using data from similar sequences in the reference genome. Annotations are transferred in coordinates relative to each of the mapped contigs. The input annotation data have to be sorted and indexed with tabix. Multiple contig mappings and multiple reference annotations can be used.

# use multiple mappings like this:
# COORDS="30-tg-gmap/lu_master300_v2.gmap.gff3 31-tg-sim4db/lu_master500_v2.fna.filtered.gff3"
# use multiple annots like this:
# ANNOTS="$GENOMEDIR/annot/ensGene_s.bed.gz $GENOMEDIR/annot/refSeq_s.bed.gz"
COORDS=30-tg-gmap/454Isotigs.gmap.gff3
ANNOTS=$GENOMEDIR/ensGene.sorted.gz
OUT=32-liftover
mkdir -p $OUT

for C in $COORDS
do
    liftover.py "$C" $ANNOTS > $OUT/${C##*/}-lo.gff3
done

Create a ‘transcript scaffold’ using the annotations¶

Construct a ‘transcript scaffold’ (contigs joined in their order of appearance on reference genome chromosomes). This is mainly because of viewing convenience with IGV. ‘N’ gaps should be larger than the longest read size to avoid mapping of the reads across gaps:

# filtered contigs
INFILE=20-newbler/454Isotigs.fna.filtered
# transferred annotations from previous step
ANNOTS=32-liftover/*-lo.gff3
# output directory
OUT=33-scaffold
# name of the output 'genome'
GNAME=sc-demo

mkdir $OUT
OUTGFF=$OUT/$GNAME.gff3

scaffold.py $INFILE $ANNOTS $OUT/$GNAME.fasta $OUTGFF

# index the new genome
samtools faidx $OUT/$GNAME.fasta

# sort, compress and index the merged annotations
# so they can be used further down in the pipeline
OUTFILE=${OUTGFF%.*}.sorted.gff3

sortBed -i $OUTGFF > $OUTFILE
bgzip $OUTFILE
tabix -p gff $OUTFILE.gz

The transcript scaffold with the sorted .sorted.gff3 is the first thing worth loading to IGV.

Map all the reads using smalt¶

Set up variables:

# data from previous steps
SCAFFOLD=33-scaffold/sc-demo.fasta
INFILES=10-cutadapt/*.fastq
OUT=40-map-smalt

SMALT_IDX=${SCAFFOLD%/*}/smalt/${SCAFFOLD##*/}-k13s4

Create index for the scaffold and map the reads. Mapping 3 GB of reads (fastq format) takes ~5 hours in 8 threads on Intel Xeon E5620, 0.5 GB memory per each mapping. This step would benefit from parallelization even on multiple machines (not implemented here).

# create smalt index
mkdir -p ${SMALT_IDX%/*}
smalt index -s 4 $SMALT_IDX $SCAFFOLD

# map each file, smalt is multithreaded so feed the files sequentially
mkdir -p $OUT
parallel -j 1 "smalt map -n $CPUS -p -f sam -o $OUT/{/.}.sam $SMALT_IDX {}" ::: $INFILES

# Illumina reads can be maped e.g. with bwa
bwa index $SCAFFOLD
parallel -j 1 "bwa mem -t $CPUS $SCAFFOLD {} > $OUT/{/.}.sam" ::: $INFILES

samtools faidx $SCAFFOLD

Create readgroups.txt¶

According to your sample wet lab details, create a readgroups.txt file. Because samtools merge -r attaches read group to each alignment (line) in the input according to the original filename, the format is ($BASENAME is the fastq file name without suffix, $SAMPLE is your biological sample, ${BASENAME%%.*} is the dna library name, all tab separated):

@RG ID:$BASENAME    SM:$SAMPLE      LB:${BASENAME%%.*}      PL:LS454 DS:$SPECIES

The library name (LB) is important because of rmdup, description (DS) is here used to identify the species.

The following code generates most of the readgroups.txt file, you have to reorder lines and fill the places marked with ‘??’:

OUT=40-map-smalt
DIR=10-cutadapt

find $DIR -name '*.fastq' | xargs -n1 basename | sed s/.fastq// | awk '{OFS="\t";lb=$0;sub(/\..*$/,"",lb);print "@RG", "ID:" $0, "SM:??", "LB:" lb, "PL:LS454", "DS:??";}' > $OUT/readgroups.txt

# edit the file (ctrl-o enter ctrl-x saves and exits the editor)
nano $OUT/readgroups.txt

# sort the readgroups according to species
<$OUT/readgroups.txt sort -k6,6 > $OUT/readgroups-s.txt

Prepare the sam files¶

Extract the sequence headers from the first .sam file (other files should have identical headers):

SAMFILE=$( echo $OUT/*.sam | awk '{print $1;}' )
samtools view -S -t $SCAFFOLD.fai -H $SAMFILE > $OUT/sequences.txt
cat $OUT/sequences.txt $OUT/readgroups-s.txt > $OUT/sam-header.txt

samtools merge requires sorted alignments, sort them in parallel. This creates .bam files in the output directory:

parallel -j $CPUS "samtools view -but $SCAFFOLD.fai {} | samtools sort - {.}" ::: $OUT/*.sam

Merge¶

Merge all the alignments. Do not remove duplicates because the duplicate detection algorithm is based on the read properties of genomic DNA ([1], [2]).

samtools merge -ru -h $OUT/sam-header.txt - $OUT/*.bam | samtools sort - $OUT/alldup
samtools index $OUT/alldup.bam

parallel -j $CPUS 'echo $( cut -f2 {}|grep -c "^4$" ) {}' ::: $OUT/*.sam

samtools idxstats $OUT/alldup.bam | awk '{map += $3; unmap += $4;} END {print  unmap/map;}'

igvtools count -z 5 -w 25 -e 250 $OUT/alldup.bam  $OUT/alldup.bam.tdf ${CONTIGS%.*}.genome

Call variants with samtools¶

Set up input and output for the current step:

SCAFFOLD=33-scaffold/sc-demo.fasta
ALIGNS=40-map-smalt/alldup.bam

# outputs
OUT=50-variants
VARIANTS=$OUT/demo-variants
mkdir -p $OUT

Run the variant calling in parallel. Takes 3 hours for 15 samples on a single Intel Xeon E5620:

vcfutils.pl splitchr $SCAFFOLD.fai | parallel -j $CPUS "samtools mpileup -DSu -L 10000 -f $SCAFFOLD -r {} $ALIGNS | bcftools view -bvcg - > $OUT/part-{}.bcf"

# samtools > 0.1.19
vcfutils.pl splitchr $SCAFFOLD.fai | parallel -j $CPUS "samtools mpileup -u -t DP,SP -L 10000 -f $SCAFFOLD -r {} $ALIGNS | bcftools call -O b -vm - > $OUT/part-{}.bcf"

Merge the intermediate results. vcfutils.pl is used to generate the correct ordering of parts, as in the .fai:

vcfutils.pl splitchr $SCAFFOLD.fai | xargs -i echo $OUT/part-{}.bcf | xargs -x bcftools cat > $VARIANTS-raw.bcf
bcftools index $VARIANTS-raw.bcf

# samtools > 0.1.19
# workaround for https://github.com/samtools/bcftools/issues/134
vcfutils.pl splitchr $SCAFFOLD.fai | parallel "bcftools view -O z $OUT/part-{}.bcf > $OUT/{.}.vcf.gz && bcftools index $OUT/{.}.vcf.gz"
bcftools concat -O b $OUT/*.vcf.gz > $VARIANTS-raw.bcf
bcftools index $VARIANTS-raw.bcf

Remove the intermediate results, if the merge was ok:

vcfutils.pl splitchr $SCAFFOLD.fai | xargs -i echo $OUT/part-{}.bcf | xargs rm

Call variants with FreeBayes¶

This is an alternative to the previous section. FreeBayes uses local realignment around INDELs, so the variant calling for 454 data should be better.

SCAFFOLD=33-scaffold/sc-demo.fasta
ALIGNS=40-map-smalt/alldup.bam
OUT=51-var-freebayes

GNAME=$( echo ${SCAFFOLD##*/} | cut -d. -f1 )
mkdir -p $OUT
vcfutils.pl splitchr $SCAFFOLD.fai | parallel -j $CPUS "freebayes -f $SCAFFOLD -r {} $ALIGNS > $OUT/${GNAME}-{}.vcf"

# join the results
OFILE=$OUT/variants-raw.vcf

# headers
FILE=$( find $OUT -name ${GNAME}-*.vcf | sort | head -1 )
<$FILE egrep '^#' > $OFILE
# the rest in .fai order
vcfutils.pl splitchr $SCAFFOLD.fai | parallel -j 1 "egrep -v '^#' $OUT/${GNAME}-{}.vcf >> $OFILE"

# filter the variants on quality (ignore the warning messages)
<$OFILE bcftools view -i 'QUAL > 20' -O z > $OUT/variants-qual.vcf.gz
tabix -p vcf $OUT/variants-qual.vcf.gz

Filter variants¶

We’re interested in two kinds of variant qualities

• all possible variants so they can be avoided in primer design
• high confidence variants that can be used to answer our questions

Filtering strategy:

• use predefined samtools filtering
  • indels caused by 454 homopolymer problems generally have low quality scores, so they should be filtered at this stage
• remove uninteresting information (for convenient viewing in IGV)
  • overall low coverage sites (less than 3 reads per sample - averaged, to avoid discarding some otherwise interesting information because of one bad sample)
• select the interesting variants, leave the rest in the file flagged as ‘uninteresting’
  • only SNPs
  • at least 3 reads per sample
  • no shared variants between the two species

bcftools view $VARIANTS-raw.bcf | vcfutils.pl varFilter -1 0 | bgzip > $VARIANTS-filtered.vcf.gz
tabix -p vcf $VARIANTS-filtered.vcf.gz

# filter on average read depth and site quality
VCFINPUT=$VARIANTS-filtered.vcf.gz
VCFOUTPUT=$VARIANTS-filt2.vcf.gz
pv -p $VCFINPUT | bgzip -d | vcf_filter.py --no-filtered - avg-dps --avg-depth-per-sample 5 sq --site-quality 30 | bgzip > $VCFOUTPUT

# samtools > 0.1.19 produce conflicting info tags, get rid of it if the above filtering fails
pv -p $VCFINPUT | bgzip -d | sed 's/,Version="3"//' | vcf_filter.py --no-filtered - avg-dps --avg-depth-per-sample 5 sq --site-quality 30 | bgzip > $VCFOUTPUT

# freebayes output
zcat $OUT/variants-qual.vcf.gz| vcfutils.pl varFilter -1 0 | vcf_filter.py --no-filtered - avg-dps --avg-depth-per-sample 5 sq --site-quality 30 | bgzip > $OUT/demo-filt1.vcf.gz

# samtools files
VCFINPUT=$VARIANTS-filt2.vcf.gz
VCFOUTPUT=$VARIANTS-selected.vcf.gz

# freebayes files
VCFINPUT=$OUT/demo-filt1.vcf.gz
VCFOUTPUT=$OUT/demo-selected.vcf.gz

<$VCFINPUT pv -p | zcat | vcf_filter.py - dps --depth-per-sample 3 snp-only contrast-samples --sample-names lu02 lu14 lu15 | bgzip > $VCFOUTPUT
tabix -p vcf $VCFOUTPUT

Check the results¶

Extract calculated variant qualities, so the distribution can be checked (-> common power law distribution, additional peak at 999):

zcat $VCFINPUT | grep -v '^#' | cut -f6 > $VCFINPUT.qual
# and visualize externally

# or directly in terminal, using bit.ly data_hacks tools
zcat $VCFINPUT | grep -v '^#' | cut -f6 | histogram.py -b 30

Count selected variants:

zcat -d $VCFOUTPUT | grep -c PASS

Count variants on chromosome Z:

zcat -d $VCFOUTPUT | grep PASS | grep -c ^chrZ

Design primers with Primer3¶

Set inputs and outputs for this step:

VARIANTS=50-variants/demo-variants-selected.vcf.gz
SCAFFOLD=33-scaffold/sc-demo.fasta
ANNOTS=33-scaffold/sc-demo.sorted.gff3.gz

GFFILE=demo-primers.gff3
OUT=60-gff-primers

GFF=$OUT/$GFFILE
PRIMERS=${GFF%.*}.sorted.gff3.gz
mkdir -p $OUT

# for all selected variants design pcr and genotyping primers
# takes about a minute for 1000 selected variants, 5 MB gzipped vcf, 26 MB uncompressed genome, 5 MB gzipped gff
# default values are set for SNaPshot
export PRIMER3_CONFIG=/opt/primer3/primer3_config/
design_primers.py $SCAFFOLD $ANNOTS $VARIANTS > $GFF

# use --primer-pref to set preferred length of genotyping primer
# this is useful for other genotyping methods, like MALDI-TOF
design_primers.py --primer-pref 15 --primer-max 25 $SCAFFOLD $ANNOTS $VARIANTS > $GFF

Sort and index the annotation before using it in IGV. For a small set of primers it is not necessary to compress and index the file, IGV can handle raw files as well.

sortBed -i $GFF | bgzip > $PRIMERS
tabix -f -p gff $PRIMERS

Create a region list for IGV to quickly inspect all the primers.

<$GFF awk 'BEGIN{OFS="\t";} /pcr-product/ {match($9, "ID=[^;]+"); print $1, $4, $5, substr($9, RSTART+3, RLENGTH);}' > ${GFF%.*}.bed

Convert scaffold to the blat format¶

SCAFFOLD2BIT=$OUT/${SCAFFOLD##*/}.2bit
faToTwoBit $SCAFFOLD $SCAFFOLD2BIT

Validate primers with blat/isPcr¶

Recommended parameters for PCR primers in blat [1]: -tileSize=11, -stepSize=5

Get the primer sequences, in formats for isPcr and blat:

PRIMERS_BASE=${PRIMERS%%.*}
extract_primers.py $PRIMERS isPcr > $PRIMERS_BASE.isPcr
extract_primers.py $PRIMERS > $PRIMERS_BASE.fa

Check against transcriptome data and the reference genome:

# select one of those a time:
TARGET=$SCAFFOLD2BIT
TARGET=$GENOME2BIT

TARGET_TAG=$( basename ${TARGET%%.*} )
isPcr -out=psl $TARGET $PRIMERS_BASE.isPcr $PRIMERS_BASE.isPcr.$TARGET_TAG.psl
blat -minScore=15 -tileSize=11 -stepSize=5 -maxIntron=0 $TARGET $PRIMERS_BASE.fa $PRIMERS_BASE.$TARGET_TAG.psl

Check the results¶

Count and check all places where primer3 reported problems:

<$GFF grep gt-primer | grep -c 'PROBLEMS='
<$GFF grep gt-primer | grep 'PROBLEMS=' | less -S

# count unique variants with available primer set
<$GFF grep gt-primer|grep -v PROBLEM|egrep -o 'ID=[^;]+'|cut -c-13|sort -u|wc -l

Use agrep to find similar sequences in the transcript scaffold, to check if the sensitivity settings of blat are OK. Line wrapping in fasta can lead to false negatives, but at least some primers should yield hits:

# agrep is quite enough for simple checks on assemblies of this size (30 MB)
SEQ=GCACATTTCATGGTCTCCAA
agrep $SEQ $SCAFFOLD|grep $SEQ

Import your primers to any spreadsheet program with some selected information on each primer. Use copy and paste, the file format is tab separated values. When there is more than one genotyping primer for one PCR product, the information on the PCR product is repeated.

extract_primers.py $PRIMERS table > $PRIMERS_BASE.tsv

How to get to this view¶

• run IGV (version 2.2 is used here)
• Genomes > Load genome from file, choose your scaffold
• load all the tracks by File > Load from File
• rearrange tracks according to your preference by dragging the label with the mouse
• choose Regions > Import Regions, pick the .bed file created in Design primers, choose Regions > Region Navigator

Components bound to the pipeline¶

Tools for FASTA/FASTQ manipulation¶

General purpose tools¶