Skip to content

Workflow Examples

Jump to bottom Edit New page
Isaac Turner edited this page Feb 19, 2015 · 44 revisions

Table of Contents

## Graph Construction

Construct graphs for three samples (using 70GB ram):

ctx31 build -m 70G -k 31 --sample NA12878 --seq input.fq.gz NA12878.ctx
ctx31 build -m 70G -k 31 --sample Mickey --seq2 reads.1.fq.gz reads.2.fq.gz Mickey.ctx
ctx31 build -m 70G -k 31 --sample Minnie --seq data.bam Minnie.ctx

Construct graph for the reference (hg19)

ctx31 build -m 70G -k 31 --sample hg19 --seq hg19.fa.gz hg19.ctx
## Merging Graphs

Merge multiple graphs together. Graphs are kept as separate colours unless the --flatten option is used. To combine three graphs into refAndSamples.ctx:

ctx31 join --out refAndSamples.ctx NA12878.ctx Mickey.ctx Minnie.ctx
## Error Cleaning There are two options here depending on sequence coverage.

A) 'Clean' graphs to remove sequencing error (per sample, for high coverage samples):

ctx31 clean NA12878.clean.ctx NA12878.ctx
ctx31 clean Mickey.clean.ctx Mickey.ctx
ctx31 clean Minnie.clean.ctx Minnie.ctx

...then merge graphs into file refAndSamples.ctx. Uses 80GB ram:

ctx31 join -m 80G refAndSamples.clean.ctx hg19.ctx NA12878.clean.ctx Mickey.clean.ctx Minnie.clean.ctx

B) Alternatively merge uncleaned samples then clean on the population (multiple low depth samples). Uses 100GB ram:

ctx31 clean -m 100G --out samples.clean.ctx NA12878.ctx Mickey.ctx Minnie.ctx
ctx31 join -m 80G refAndSamples.ctx hg19.ctx samples.clean.ctx

Now we have cortex graphs of the reference and our samples with sequencing error removed.

Plot Coverage Distribution

The command ctx31 clean has the option of writing out histograms of supernode lengths and coverages with the --covg-before, --covg-after, --len-before and --len-after arguments. You can plot the results of these using a pair of scripts bundled with cortex in the scripts/ directory (requires R and R packages ggplot2 and gridExtra):

R --vanilla --file=scripts/plot-covg-hist.R --args covg.before_cleaning.csv covg.before_cleaning.pdf
R --vanilla --file=scripts/plot-covg-hist.R --args covg.after_cleaning.csv covg.after_cleaning.pdf
R --vanilla --file=scripts/plot-length-hist.R --args length.before_cleaning.csv length.before_cleaning.pdf
R --vanilla --file=scripts/plot-length-hist.R --args length.after_cleaning.csv length.after_cleaning.pdf
## Read Threading

Read threading aligns reads against the graph and generates graph `links' which store long range connectivity information. This is per-sample information (each sample has its own set of links).

Before threading reads through the graph, you must run inferedges:

ctx31 inferedges -m 50G in.ctx

This edits in.ctx to add missing edges between kmers. Later steps assume all edges between kmers already exist.

Thread reads through the graph to assist in assembly, generating file refAndSamples.ctp.gz:

ctx31 thread -m 80G --seq2 NA12878.1.fq NA12878.2.fq --out NA12878.ctp.gz refAndSamples.ctx:1
ctx31 thread -m 80G --seqi Mickey.bam --out Mickey.ctp.gz refAndSamples.ctx:2
ctx31 thread -m 80G --seqi Minnie.sam --out Minnie.ctp.gz refAndSamples.ctx:3

Now the link files NA12878.ctp.gz, Mickey.ctp.gz and Minnie.ctp.gz are ready to be used alongside their respective graphs, e.g.:

ctx31 contigs -m 100G -p 1:NA12878.ctp.gz -p 2:Mickey.ctp.gz -p 3:Minnie.ctp.gz refAndSamples.ctx

Note: the reference is colour 0 in the above example.

## Link Cleaning

Once a set of links have been generated, we need to remove links due to sequencing error.

# In this example k=31
# First generate two CSV files using up to 1000 kmers: links.raw.effcovg.csv links.raw.list.csv
./scripts/cortex_links.pl list --limit 1000 <(gzip -fcd links.raw.ctp.gz) links.raw.effcovg.csv links.raw.list.csv
# Use links.raw.list.csv to pick a cleaning threshold (k=31)
R --slave --vanilla --quiet -f ./scripts/R/make_link_cutoffs.R --args 31 links.raw.list.csv > links.cleaning.txt
# Fetch threshold from the output
thresh=$(tail -1 links.cleaning.txt)
# Clean the links and write them to links.clean.ctp.gz
./scripts/cortex_links.pl clean <(gzip -fcd links.raw.ctp.gz) $thresh | gzip -c) > links.clean.ctp.gz

Cleaning links also removes redundant links. Combined with removing rare links, this reduces the memory required to store graph connectivity.

## Merge Link Files

Similar to the join command, pjoin merges link files.

ctx31 pjoin -o refAndNA12878.ctp.gz ref.ctp.gz NA12878.ctp.gz
## Calling Variants 
ctx31 bubbles -m 100G --out denovo.bubbles.gz --haploid 0 -p refAndSamples.ctp.gz refAndSamples.ctx

Uses the graph and links (refAndSamples.ctx, refAndSamples.ctp.gz) to call variants that are saved to the file denovo.bubbles.gz. We also specify that the first colour is a haploid sample (with --haploid 0). This is useful in removing bubbles caused by repeats rather than real variants.

Next we extract the 5' flanks from the bubble file:

./scripts/cortex_print_flanks.sh denovo.bubbles.gz > denovo.5pflank.fa

Creates the file denovo.5pflank.fa. We then map the putative variant flanks to our reference (in this case hg19). We can do this with stampy:

stampy.py -G hg19 hg19.fa  # create hg19.stidx
stampy.py -g hg19 -H hg19  # create hg19.sthash
stampy.py -g hg19 -h hg19 --inputformat=fasta -M denovo.5pflank.fa > denovo.5pflank.sam

or BWA:

bwa index hg19.fa
bwa mem hg19.fa denovo.5pflank.fa > denovo.5pflank.sam

Then we align the alleles to the reference and generate a VCF

ctx31 calls2vcf -F denovo.5pflank.sam -o denovo.vcf bubbles.txt.gz hg19.fa

This generates a VCF of sites. Genotyping variants is not yet implemented. See the VCF post processing section to sort, remove duplicates, index and clean up your VCF.

## Compare to a Reference

If you have an assembled genome, you can compare your sequenced samples to it to identify large variants and re-arrangements:

ctx31 breakpoints -m 100G --seq hg19.fa --out breakpoints.txt.gz -p PlutoAndGoofy.ctp.gz PlutoAndGoofy.ctx

The breakpoints command identifies where your samples diverge from the reference and rejoin it elsewhere. This could be small regions of clustered SNPs, inversions or translocations.

We can generate a VCF of small events using calls2vcf:

ctx31 calls2vcf -o breakpoints.vcf breakpoints.txt.gz hg19.fa

See the VCF post processing section to sort, remove duplicates, index and clean up your VCF.

Results can be plotted on a Circos plot using a perl script in scripts/:

gzip -dc breakpoints.txt.gz | scripts/make-circos.pl out-dir -
cd out-dir
circos

Plot appears in out-dir/circos.png.

## Post processing a VCF

Once you have called variants with the Bubble Caller or the Breakpoint Caller, you are left with an unsorted VCF calls.raw.vcf. Our proposed post-processing sorts, removes duplicates, renames variants, compresses and indexes the VCF:

./scripts/bash/vcf-sort calls.raw.vcf > calls.sorted.vcf
bcftools norm --remove-duplicates --fasta-ref ref.fa --multiallelics +both $< | \
  ./scripts/bash/vcf-rename > calls.norm.vcf
bgzip calls.norm.vcf
bcftools index calls.norm.vcf.gz
## Assemble Contigs

Pull out contigs from an individual in a population. Fills in low coverage gaps in sequence from the population if unconfounded.

ctx31 contigs -p refAndSamples.ctp.gz --colour 1 --out refAndSamples.contigs.fa refAndSamples.ctx

The contigs command uses random kmers to seed contigs, therefore it may print duplicate contigs and some contigs that are substrings of others. To remove duplicate contigs you can use the rmsubstr command:

ctx31 rmsubstr -m 10G --out refAndSamples.contigs.rmdup.fa refAndSamples.contigs.fa

More on the contig assembly page.

## Filter A Subgraph

Pull out a region of interest using a set of reads or kmers. The subgraph command searches out from the seed kmers by dist kmers (dist >= 0).

ctx31 subgraph --seq contig.fa contig.ctx 100 graph.ctx

Input is graph.ctx, output is saved to contig.ctx. You can specify particular colours on the input to keep in the output. Add the --invert option to invert the kmers selected - only take those further than dist kmers.

## Filter Reads Against A Graph

Filter reads by whether or not they share a kmer with a graph:

ctx31 reads [--fasta|--fastq] --seq in1.fa out.se --seq2 in.1.fq in.2.fq out.pair in.ctx

Will save reads from in1.fa to out.se.fa.gz and from in.1.fq,in.2.fq to out.pair.1.fa.gz,out.pair.2.fa.gz. Output options are --fasta or --fastq. Paired end reads (--seq2) are considered touching the graph if either read shares one or more kmers with the graph. Add the --invert option to select reads that DO NOT touch the graph.

## Plot Graph

Requires GraphViz which provides the dot command. Flattens any colours in the graph:

./scripts/cortex_to_graphviz.pl graph.ctx:0,3 | dot -Tpdf > kmers.pdf

Or to merge adjacent kmers add the --simplify option:

./scripts/cortex_to_graphviz.pl --simplify graph.ctx:0,3 | dot -Tpdf > contigs.pdf

Or you can just use ./scripts/cortex_to_graphviz.pl to generate dot files.

You can generate a plot of a de Bruijn graph straight from sequence files using scripts/seq2pdf.sh:

./scripts/seq2pdf.sh --simplify 3 <(echo ACAACACGT) <(echo CCACACAA) > out.pdf

Usage is: ./scripts/seq2pdf.sh [--simplify] <kmer> <file1> [file2] > out.pdf. Files can be of types FASTA, FASTQ, SAM, BAM, TXT or even gzipped. Requires graphviz.

Add a custom footer
Add a custom sidebar
Clone this wiki locally