NAME

bcbio_nextgen - bcbio_nextgen Documentation [image: bcbio banner] [image]

A python toolkit providing best-practice pipelines for fully automated high throughput sequencing analysis. You write a high level configuration file specifying your inputs and analysis parameters. This input drives a parallel pipeline that handles distributed execution, idempotent processing restarts and safe transactional steps. The goal is to provide a shared community resource that handles the data processing component of sequencing analysis, providing researchers with more time to focus on the downstream biology.

Contents:

bcbio-nextgen provides best-practice pipelines for automated analysis of high throughput sequencing data with the goal of being:

USERS

A sample of institutions using bcbio-nextgen for solving biological problems. Please submit your story if you're using the pipeline in your own research.

AUTOMATED

We provide an automated script that installs third party analysis tools, required genome data and python library dependencies for running human variant and RNA-seq analysis, bundled into an isolated directory or virtual environment:

wget https://raw.github.com/bcbio/bcbio-nextgen/master/scripts/bcbio_nextgen_install.py
python bcbio_nextgen_install.py /usr/local/share/bcbio --tooldir=/usr/local \


  --genomes GRCh37 --aligners bwa --aligners bowtie2

bcbio should install cleanly on Linux systems. For Mac OSX, we suggest trying bcbio-vm which runs bcbio on docs-cloud or isolates all the third party tools inside a Docker container. bcbio-vm is still a work in progress but not all of the dependencies bcbio uses install cleanly on OSX.

With the command line above, indexes and associated data files go in /usr/local/share/bcbio-nextgen and tools are in /usr/local. If you don't have write permissions to install into the /usr/local directories you can install in a user directory like ~/local or use sudo chmod to give your standard user permissions. Please don't run the installer with sudo or as the root user. Do not use directories with : in the name, it is not POSIX compliant and will cause installation failures.

The installation is highly customizable, and you can install additional software and data later using bcbio_nextgen.py upgrade. Run python bcbio_nextgen_install.py with no arguments to see options for configuring the installation process. Some useful arguments are:

export PATH=/path_to_bcbio/bin:$PATH

The machine will need to have some basic requirements for installing and running bcbio:

Optional tool specific requirements:

The bcbio-nextgen Dockerfile contains the packages needed to install on bare Ubuntu systems.

The automated installer creates a fully integrated environment that allows simultaneous updates of the framework, third party tools and biological data. This offers the advantage over manual installation of being able to manage and evolve a consistent analysis environment as algorithms continue to evolve and improve. Installing this way is as isolated and self-contained as possible without virtual machines or lightweight system containers like Docker. The Upgrade section has additional documentation on including additional genome data for supported bcbio genomes. For genome builds not included in the defaults, see the documentation on config-custom. Following installation, you should edit the pre-created system configuration file in /usr/local/share/bcbio-nextgen/galaxy/bcbio_system.yaml to match your local system or cluster configuration (see tuning-cores).

UPGRADE

We use the same automated installation process for performing upgrades of tools, software and data in place. Since there are multiple targets and we want to avoid upgrading anything unexpectedly, we have specific arguments for each. Generally, you'd want to upgrade the code, tools and data together with:

bcbio_nextgen.py upgrade -u stable --tools --data

Tune the upgrade with these options:

CUSTOMIZING DATA INSTALLATION

bcbio installs associated data files for sequence processing, and you're able to customize this to install larger files or change the defaults. Use the --datatarget flag (potentially multiple times) to customize or add new targets.

By default, bcbio will install data files for variation, rnaseq and smallrna but you can sub-select a single one of these if you don't require other analyses. The available targets are:

For somatic analyses, bcbio includes COSMIC v68 for hg19 and GRCh37 only. Due to license restrictions, we cannot include updated versions of this dataset and hg38 support with the installer. To prepare these datasets yourself you can use a utility script shipped with cloudbiolinux that downloads, sorts and merges the VCFs, then copies into your bcbio installation:

export COSMIC_USER="your@registered.email.edu"
export COSMIC_PASS="cosmic_password"
bcbio_python prepare_cosmic.py 83 /path/to/bcbio

EXTRA SOFTWARE

We're not able to automatically install some useful tools due to licensing restrictions, so we provide a mechanism to manually download and add these to bcbio-nextgen during an upgrade with the --toolplus command line.

gatk-register /path/to/GenomeAnalysisTK.tar.bz2

bcbio_nextgen.py upgrade --tools --toolplus gatk=/path/to/gatk/GenomeAnalysisTK.jar

bcbio_nextgen.py upgrade --tools --toolplus mutect=/path/to/mutect/mutect-1.1.7.jar

SYSTEM REQUIREMENTS

bcbio-nextgen provides a wrapper around external tools and data, so the actual tools used drive the system requirements. For small projects, it should install on workstations or laptops with a couple Gb of memory, and then scale as needed on clusters or multicore machines.

Disk space requirements for the tools, including all system packages are under 4Gb. Biological data requirements will depend on the genomes and aligner indices used, but a suggested install with GRCh37 and bowtie/bwa2 indexes uses approximately 35Gb of storage during preparation and ~25Gb after:

$ du -shc genomes/Hsapiens/GRCh37/*
3.8G  bowtie2
5.1G  bwa
3.0G  rnaseq-2014-05-02
3.0G  seq
340M  snpeff
4.2G  variation
4.4G  vep
23.5G total

TROUBLESHOOTING

$ git config --global url.https://github.com/.insteadOf git://github.com/

Unsupported major.minor version 51.0

$ unset PYTHONHOME
$ unset PYTHONPATH
$ export PYTHONNOUSERSITE=1

MANUAL PROCESS

The manual process does not allow the in-place updates and management of third party tools that the automated installer makes possible. It's a more error-prone and labor intensive process. If you find you can't use the installer we'd love to hear why to make it more amenable to your system. If you'd like to develop against a bcbio installation, see the documentation on setting up a code-devel-infrastructure.

fab -f cloudbiolinux/fabfile.py -H localhost install_biolinux:flavor=ngs_pipeline_minimal

`packages-conda.yaml`_

fab -f data_fabfile.py -H localhost -c your_fabricrc.txt install_data_s3:your_biodata.yaml

├── galaxy
│   ├── tool-data
│   │   ├── alignseq.loc
│   │   ├── bowtie_indices.loc
│   │   ├── bwa_index.loc
│   │   ├── sam_fa_indices.loc
│   │   └── twobit.loc
│   └── tool_data_table_conf.xml
├── genomes
│   ├── Hsapiens
│   │   ├── GRCh37
│   │   └── hg19
│   └── phiX174
│       └── phix
└── liftOver

`-- phix


    |-- bowtie


    |   |-- phix.1.ebwt


    |   |-- phix.2.ebwt


    |   |-- phix.3.ebwt


    |   |-- phix.4.ebwt


    |   |-- phix.rev.1.ebwt


    |   `-- phix.rev.2.ebwt


    |-- bowtie2


    |   |-- phix.1.bt2


    |   |-- phix.2.bt2


    |   |-- phix.3.bt2


    |   |-- phix.4.bt2


    |   |-- phix.rev.1.bt2


    |   `-- phix.rev.2.bt2


    |-- bwa


    |   |-- phix.fa.amb


    |   |-- phix.fa.ann


    |   |-- phix.fa.bwt


    |   |-- phix.fa.pac


    |   |-- phix.fa.rbwt


    |   |-- phix.fa.rpac


    |   |-- phix.fa.rsa


    |   `-- phix.fa.sa


    |-- novoalign


    |   `-- phix


    |-- seq


    |   |-- phix.dict


    |   |-- phix.fa


    |   `-- phix.fa.fai


    `-- ucsc


        `-- phix.2bit

GERMLINE VARIANT CALLING

bcbio implements configurable SNP, indel and structural variant calling for germline populations. We include whole genome and exome evaluations against reference calls from the Genome in a Bottle consortium and Illumina Platinum Genomes project, enabling continuous assessment of new alignment and variant calling algorithms. We regularly report on these comparisons and continue to improve approaches as the community makes new tools available. Here is some of the research that contributes to the current implementation:

bcbio automates post-variant calling annotation to make the outputs easier to feed directly into your biological analysis. We annotate variant effects using snpEff or Variant Effect Predictor (VEP), and prepare a GEMINI database that associates variants with multiple external annotations in a SQL-based query interface. GEMINI databases have the most associated external information for human samples (GRCh37/hg19 and hg38) but are available for any organism with the database populated using the VCF INFO column and predicted effects.

- description: Sample1


  metadata:


    batch: Batch1
- description: Sample2


  metadata:


    batch: Batch1

- description: Sample1


  algorithm:


    variantcaller: gatk-haplotype


    jointcaller: gatk-haplotype-joint


  metadata:


    batch: Batch1
- description: Sample2


  algorithm:


    variantcaller: gatk-haplotype


    jointcaller: gatk-haplotype-joint


  metadata:


    batch: Batch1

CANCER VARIANT CALLING

bcbio supports somatic cancer calling with tumor and optionally matched normal pairs using multiple SNP, indel and structural variant callers. A full evaluation of cancer calling validates callers against synthetic dataset 3 from the ICGC-TCGA DREAM challenge. bcbio uses a majority voting ensemble approach to combining calls from multiple SNP and indel callers, and also flattens structural variant calls into a combined representation.

The example configuration for the example-cancer validation is a good starting point for setting up a tumor/normal run on your own dataset. The configuration works similarly to population based calling. Supply a consistent batch for tumor/normal pairs and mark them with the phenotype:

- description: your-tumor


  algorithm:


    variantcaller: [vardict, strelka2, mutect2]


  metadata:


    batch: batch1


    phenotype: tumor
- description: your-normal


  algorithm:


    variantcaller: [vardict, strelka2, mutect2]


  metadata:


    batch: batch1


    phenotype: normal

Other config-cancer configuration options allow tweaking of the processing parameters. For pairs you want to analyze together, specify a consistent set of variantcaller options for both samples.

Cancer calling handles both tumor-normal paired calls and tumor-only calling. To specify a tumor-only sample, provide a single sample labeled with phenotype: tumor. Otherwise the configuration and setup is the same as with paired analyses. For tumor-only samples, bcbio will try to remove likely germline variants present in the public databases like 1000 genomes and ExAC, and not in COSMIC. This runs as long as you have a local GEMINI data installation (--datatarget gemini) and marks likely germline variants with a LowPriority filter. This post has more details on the approach and validation.

The standard variant outputs (sample-caller.vcf.gz) for tumor calling emphasize somatic differences, those likely variants unique to the cancer. If you have a tumor-only sample and GEMINI data installed, it will also output sample-caller-germline.vcf.gz, which tries to identify germline background mutations based on presence in public databases. If you have tumor/normal data and would like to also call likely germline mutations see the documentation on specifying a germline caller: Somatic with germline variants.

We're actively working on improving calling to better account for the heterogeneity and structural variability that define cancer genomes.

SOMATIC WITH GERMLINE VARIANTS

For tumor/normal somatic samples, bcbio can call both somatic (tumor-specific) and germline (pre-existing) variants. The typical outputs of Cancer variant calling are likely somatic variants acquired by the cancer, but pre-existing germline risk variants are often also diagnostic.

For tumor-only cases we suggest running standard Cancer variant calling. Tumor-only inputs mix somatic and germline variants, making it difficult to separate events. For small variants (SNPs and indels) bcbio will attempt to distinguish somatic and germline mutations using the presence of variants in population databases.

To option somatic and germline calls for your tumor/normal inputs, specify which callers to use for each step in the variant-config configuration:

description: your-normal
variantcaller:


   somatic: vardict


   germline: freebayes

bcbio does a single alignment for the normal sample, then splits at the variant calling steps using this normal sample to do germline calling. In this example, the output files are:

Germline calling supports multiple callers, and other configuration options like ensemble and structural variant calling inherit from the remainder configuration. For example, to use 3 callers for somatic and germline calling, create ensemble calls for both and include germline and somatic events from two structural variant callers:

variantcaller:


   somatic: [vardict, strelka2, mutect2]


   germline: [freebayes, gatk-haplotype, strelka2]
ensemble:


   numpass: 2
svcaller: [manta, cnvkit]

In addition to the somatic and germline outputs attached to the tumor and normal sample outputs as described above, you'll get:

STRUCTURAL VARIANT CALLING

bcbio can detect larger structural variants like deletions, insertions, inversions and copy number changes for both germline population and cancer variant calling, based on validation against existing truth sets:

To enable structural variant calling, specify svcaller options in the algorithm section of your configuration:

- description: Sample


  algorithm:


    svcaller: [lumpy, manta, cnvkit]

The best supported callers are Lumpy and Manta, for paired end and split read calling, CNVkit for read-depth based CNV calling, and WHAM for association testing. We also support DELLY, another excellent paired end and split read caller, although it is slow on large whole genome datasets.

In addition to results from individual callers, bcbio can create a summarized ensemble callset using MetaSV. We're actively working on improved structural variant reporting to highlight potential variants of interest.

RNA-SEQ

bcbio can also be used to analyze RNA-seq data. It includes steps for quality control, adapter trimming, alignment, variant calling, transcriptome reconstruction and post-alignment quantitation at the level of the gene and isoform.

We recommend using the STAR aligner for all genomes where there are no alt alleles. For genomes such as hg38 that have alt alleles, hisat2 should be used as it handles the alts correctly and STAR does not yet. Use Tophat2 only if you do not have enough RAM available to run STAR (about 30 GB).

Our current recommendation is to run adapter trimming only if using the Tophat2 aligner. Adapter trimming is very slow, and aligners that soft clip the ends of reads such as STAR and hisat2, or algorithms using pseudoalignments like Sailfish handle contaminant sequences at the ends properly. This makes trimming unnecessary. Tophat2 does not perform soft clipping so if using Tophat2, trimming must still be done.

Salmon, which is an extremely fast alignment-free method of quantitation, is run for all experiments. Salmon can accurately quantitate the expression of genes, even ones which are hard to quantitate with other methods (see this paper for example for Sailfish, which performs similarly to Salmon.). Salmon can also quantitate at the transcript level which can help gene-level analyses (see this paper for example). We recommend using the Salmon quantitation rather than the counts from featureCounts to perform downstream quantification.

Although we do not recommend using the featureCounts based counts, the alignments are still useful because they give you many more quality metrics than the quasi-alignments from Salmon.

After a bcbio RNA-seq run there will be in the upload directory a directory for each sample which contains a BAM file of the aligned and unaligned reads, a sailfish directory with the output of Salmon, including TPM values, and a qc directory with plots from FastQC and qualimap.

In addition to directories for each sample, in the upload directory there is a project directory which contains a YAML file describing some summary statistics for each sample and some provenance data about the bcbio run. In that directory is also a combined.counts file with the featureCounts derived counts per cell.

FAST RNA-SEQ

This mode of bcbio-nextgen quantitates transcript expression using Salmon and does nothing else. It is an order of magnitude faster or more than running the full RNA-seq analysis. The cost of the increased speed is that you will have much less information about your samples at the end of the run, which can make troubleshooting trickier. Invoke with analysis: fastrna-seq.

SINGLE-CELL RNA-SEQ

bcbio-nextgen supports universal molecular identifiers (UMI) based single-cell RNA-seq analyses. If your single-cell prep does not use universal molecular identifiers (UMI), you can most likely just run the standard RNA-seq pipeline and use the results from that. The UMI are used to discard reads which are possibly PCR duplicates and is very helpful for removing some of the PCR duplicate noise that can dominate single-cell experiments.

Unlike the standard RNA-seq pipeline, the single-cell pipeline expects the FASTQ input files to not be separated by cellular barcode, so each file is a mix of cells identified by a cellular barcode (CB), and unique reads from a transcript are identified with a UMI. bcbio-nextgen inspects each read, identifies the cellular barcode and UMI and puts them in the read name. Then the reads are aligned to the transcriptome with RapMap and the number of reads aligning to each transcript is counted for each cellular barcode. The output is a table of counts with transcripts as the rows and columns as the cellular barcodes for each input FASTQ file.

Optionally the reads can be quantitated with kallisto to output transcript compatibility counts rather than counts per gene (TCC paper).

To extract the UMI and cellular barcodes from the read, bcbio-nextgen needs to know where the UMI and the cellular barcode are expected to be in the read. Currently there is support for two schemes, the inDrop system from the Harvard single-cell core facility and CEL-seq. If bcbio-nextgen does not support your UMI and barcoding scheme, please open up an issue and we will help implement support for it.

Most of the heavy lifting for this part of bcbio-nextgen is implemented in the umis repository.

SMALLRNA-SEQ

bcbio-nextgen also implements a configurable best-practices pipeline for smallRNA-seq quality controls, adapter trimming, miRNA/isomiR quantification and other small RNA detection.

The pipeline generates a _RMD_ template file inside report folder that can be rendered with knitr. An example of the report is at here. Count table (counts_mirna.tst) from mirbase miRNAs will be inside mirbase or final project folder. Input files for isomiRs package for isomiRs analysis will be inside each sample in mirbase folder.. If mirdeep2 can run, count table (counts_mirna_novel.tsv) for novel miRNAs will be inside mirdeep2 or final project folder. tdrmapper results will be inside each sample inside tdrmapper or final project folder.

CHIP-SEQ

The bcbio-nextgen implementation of ChIP-seq aligns, removes multimapping reads, calls peaks with a paired input file using MACS2 and outputs a set of greylist regions for filtering possible false peaks in regions of high depth in the input file.

STANDARD

This pipeline implements alignment and qc tools. Furthermore, it will run qsignature to detect possible duplicated samples, or mislabeling. It uses SNPs signature to create a distance matrix that helps easily to create groups. The project yaml file will show the number of total samples analyzed, the number of very similar samples, and samples that could be duplicated.

fc_date: '070113'
fc_name: control_experiment
upload:


  dir: final
details:


  - files: [/full/path/to/control_1.fastq, /full/path/to/control_2.fastq]


    description: 'Control_rep1'


    genome_build: GRCh37


    analysis: RNA-seq


    algorithm:


         aligner: tophat2


         quality_format: Standard


         trim_reads: read_through


         adapters: [truseq, polya]


         strandedness: unstranded

fc_date: '070113'
fc_name: mouse_analysis
upload:


  dir: final
details:


  - files: [/full/path/to/control_rep1.fastq]


    description: 'Control_rep1'


    genome_build: mm10


    analysis: RNA-seq


    algorithm:


         aligner: tophat2


         quality_format: Standard


         trim_reads: read_through


         adapters: [nextera, polya]


  - files: [/full/path/to/control_rep2.fastq]


    description: 'Control_rep2'


    genome_build: mm10


    analysis: RNA-seq


    algorithm:


         aligner: tophat2


         quality_format: Standard


         trim_reads: read_through


         adapters: [nextera, polya]


  - files: [/full/path/to/case_rep1.fastq]


    description: 'Case_rep1'


    genome_build: mm10


    analysis: RNA-seq


    algorithm:


         aligner: tophat2


         quality_format: Standard


         trim_reads: read_through


         adapters: [nextera, polya]


  - files: [/full/path/to/case_rep2.fastq]


    description: 'Case_rep2'


    genome_build: mm10


    analysis: RNA-seq


    algorithm:


         aligner: tophat2


         quality_format: Standard


         trim_reads: read_through


         adapters: [nextera, polya]

bcbio_nextgen.py -w template illumina-mouse-rnaseq.yaml mouse_analysis
/full/path/to/control_rep1.fastq /full/path/to/control_rep2.fastq
/full/path/to/case_rep1.fastq /full/path/to/case_rep2.fastq

OVERVIEW

bcbio_nextgen.py -w template gatk-variant project1 sample1.bam sample2_1.fq sample2_2.fq

bcbio_nextgen.py bcbio_sample.yaml -n 8

PROJECT DIRECTORY

bcbio encourages a project structure like:

my-project/
├── config
├── final
└── work

with the input configuration in the config directory, the outputs of the pipeline in the final directory, and the actual processing done in the work directory. Run the bcbio_nextgen.py script from inside the work directory to keep all intermediates there. The final directory, relative to the parent directory of the work directory, is the default location specified in the example configuration files and gets created during processing. The final directory has all of the finished outputs and you can remove the work intermediates to cleanup disk space after confirming the results. All of these locations are configurable and this project structure is only a recommendation.

LOGGING

There are 3 logging files in the log directory within your working folder:

EXAMPLE PIPELINES

We supply example input configuration files for validation and to help in understanding the pipeline.

mkdir -p NA12878-trio-eval/config NA12878-trio-eval/input NA12878-trio-eval/work
cd NA12878-trio-eval/config
wget https://raw.github.com/bcbio/bcbio-nextgen/master/config/examples/NA12878-trio-wgs-validate.yaml
cd ../input
wget https://raw.github.com/bcbio/bcbio-nextgen/master/config/examples/NA12878-trio-wgs-validate-getdata.sh
bash NA12878-trio-wgs-validate-getdata.sh
cd ../work
bcbio_nextgen.py ../config/NA12878-trio-wgs-validate.yaml -n 16

mkdir -p NA12878-trio-eval/work_joint
cd NA12878-trio-eval/config
wget https://raw.github.com/bcbio/bcbio-nextgen/master/config/examples/NA12878-trio-wgs-joint.yaml
cd ../work_joint
bcbio_nextgen.py ../config/NA12878-trio-wgs-joint.yaml -n 16

mkdir -p NA12878-exome-eval
cd NA12878-exome-eval
wget https://raw.github.com/bcbio/bcbio-nextgen/master/config/examples/NA12878-exome-methodcmp-getdata.sh
bash NA12878-exome-methodcmp-getdata.sh

cd work
bcbio_nextgen.py ../config/NA12878-exome-methodcmp.yaml -n 8

mkdir -p cancer-dream-syn3/config cancer-dream-syn3/input cancer-dream-syn3/work
cd cancer-dream-syn3/config
wget https://raw.githubusercontent.com/bcbio/bcbio-nextgen/master/config/examples/cancer-dream-syn3.yaml
cd ../input
wget https://raw.githubusercontent.com/bcbio/bcbio-nextgen/master/config/examples/cancer-dream-syn3-getdata.sh
bash cancer-dream-syn3-getdata.sh

cd ../work
bcbio_nextgen.py ../config/cancer-dream-syn3.yaml -n 8

wget https://raw.githubusercontent.com/bcbio/bcbio-nextgen/master/config/examples/cancer-giab-na12878-na24385-getdata.sh
bash cancer-giab-na12878-na24385-getdata.sh

cd work
bcbio_nextgen.py ../config/cancer-giab-na12878-na24385.yaml -n 16

mkdir -p NA12878-sv-eval
cd NA12878-sv-eval
wget https://raw.github.com/bcbio/bcbio-nextgen/master/config/examples/NA12878-sv-getdata.sh
bash NA12878-sv-getdata.sh
cd work
bcbio_nextgen.py ../config/NA12878-sv.yaml -n 16

wget https://raw.github.com/bcbio/bcbio-nextgen/master/config/examples/rnaseq-seqc-getdata.sh
bash rnaseq-seqc-getdata.sh

cd seqc/work
bcbio_nextgen.py ../config/seqc.yaml -n 8

mkdir -p NA12878-hg38-val
cd NA12878-hg38-val
wget https://raw.github.com/bcbio/bcbio-nextgen/master/config/examples/NA12878-hg38-validate-getdata.sh
bash NA12878-hg38-validate-getdata.sh
mkdir -p work
cd work
bcbio_nextgen.py ../config/NA12878-hg38-validate.yaml -n 16

├── config
│   └── NA12878-illumina.yaml
├── input
└── work

cd input
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR091/ERR091571/ERR091571_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR091/ERR091571/ERR091571_2.fastq.gz

cd config
wget https://raw.github.com/bcbio/bcbio-nextgen/master/config/examples/NA12878-illumina.yaml

cd work
bcbio_nextgen.py ../config/NA12878-illumina.yaml -n 16

TEST SUITE

The test suite exercises the scripts driving the analysis, so are a good starting point to ensure correct installation. Tests use the pytest framework. The tests are available in the bcbio source code:

$ git clone https://github.com/bcbio/bcbio-nextgen.git

There is a small wrapper script that finds the py.test and other dependencies pre-installed with bcbio you can use to run tests:

$ cd tests
$ ./run_tests.sh

You can use this to run specific test targets:

$ ./run_tests.sh cancer
$ ./run_tests.sh rnaseq
$ ./run_tests.sh devel
$ ./run_tests.sh docker

Optionally, you can run pytest directly from the bcbio install to tweak more options. It will be in /path/to/bcbio/anaconda/bin/py.test. Pass -s to py.test to see the stdout log, and -v to make py.test output more verbose. The tests are marked with labels which you can use to run a specific subset of the tests using the -m argument:

$ py.test -m rnaseq

To run unit tests:

$ py.test tests/unit

To run integration pipeline tests:

$ py.test tests/integration

To run tests which use bcbio_vm:

$ py.test tests/bcbio_vm

To see the test coverage, add the --cov=bcbio argument to py.test.

By default the test suite will use your installed system configuration for running tests, substituting the test genome information instead of using full genomes. If you need a specific testing environment, copy tests/data/automated/post_process-sample.yaml to tests/data/automated/post_process.yaml to provide a test-only configuration.

Two configuration files, in easy to write YAML format, specify details about your system and samples to run:

Commented system and sample example files are available in the config directory. The example-pipelines section contains additional examples of ready to run sample files.

AUTOMATED SAMPLE CONFIGURATION

bcbio-nextgen provides a utility to create configuration files for multiple sample inputs using a base template. Start with one of the best-practice templates, or define your own, then apply to multiple samples using the template workflow command:

bcbio_nextgen.py -w template freebayes-variant project1.csv sample1.bam sample2_1.fq sample2_2.fq

samplename,description,batch,phenotype,sex,variant_regions
sample1,ERR256785,batch1,normal,female,/path/to/regions.bed
sample2,ERR256786,batch1,tumor,,/path/to/regions.bed

samplename,description,phenotype,batch,sex,ethnicity,maternal_id,paternal_id,family_id
NA12878.bam,NA12878,-9,CEPH,female,-9,NA12892,NA12891,NA12878FAM

samplename,description,phenotype,batch
normal.bam,two_normal,normal,Batch1;Batch2

samplename,description,phenotype,variantcaller
tumor.bam,sample1,tumor,germline:freebayes;gatk-haplotype::somatic:vardict;freebayes

To make it easier to define your own project specific template, an optional first step is to download and edit a local template. First retrieve a standard template:

bcbio_nextgen.py -w template freebayes-variant project1

This pulls the current GATK best practice variant calling template into your project directory in project1/config/project1-template.yaml. Manually edit this file to define your options, then run the full template creation for your samples, pointing to this custom configuration file:

bcbio_nextgen.py -w template project1/config/project1-template.yaml project1.csv folder/*

If your sample folder contains additional BAM or FASTQ files you do not wish to include in the sample YAML configuration, you can restrict the output to only include samples in the metadata CSV with --only-metadata. The output will print warnings about samples not present in the metadata file, then leave these out of the final output YAML:

bcbio_nextgen.py -w template --only-metadata project1/config/project1-template.yaml project1.csv folder/*

MULTIPLE FILES PER SAMPLE

In case you have multiple FASTQ or BAM files for each sample you can use bcbio_prepare_samples.py. The main parameters are:

samplename,description,batch,phenotype,sex,variant_regions
file1.fastq,sample1,batch1,normal,female,/path/to/regions.bed
file2.fastq,sample1,batch1,normal,female,/path/to/regions.bed
file1.fastq,sample2,batch1,tumor,,/path/to/regions.bed

An example of usage is:

bcbio_prepare_samples.py --out merged --csv project1.csv

The script will create the sample1.fastq,sample2.fastq in the merged folder, and a new CSV file in the same folder than the input CSV :project1-merged.csv. Later, it can be used for bcbio:

bcbio_nextgen.py -w template project1/config/project1-template.yaml project1-merged.csv merged/*fastq

The new CSV file will look like:

samplename,description,batch,phenotype,sex,variant_regions
sample1.fastq,sample1,batch1,normal,female,/path/to/regions.bed
sample2.fastq,sample2,batch1,tumor,,/path/to/regions.bed

It supports parallelization the same way bcbio_nextgen.py does:

python $BCBIO_PATH/scripts/utils/bcbio_prepare_samples.py --out merged --csv project1.csv -t ipython -q queue_name -s lsf -n 1

See more examples at parallelize pipeline.

In case of paired reads, the CSV file should contain all files:

samplename,description,batch,phenotype,sex,variant_regions
file1_R1.fastq,sample1,batch1,normal,female,/path/to/regions.bed
file2_R1.fastq,sample1,batch1,normal,female,/path/to/regions.bed
file1_R2.fastq,sample1,batch1,normal,femela,/path/to/regions.bed
file2_R2.fastq,sample1,batch1,normal,female,/path/to/regions.bed

The script will try to guess the paired files the same way that bcbio_nextgen.py -w template does. It would detect paired files if the difference among two files is only _R1/_R2 or -1/-2 or _1/_2 or .1/.2

The output CSV will look like and is compatible with bcbio:

samplename,description,batch,phenotype,sex,variant_regions
sample1,sample1,batch1,normal,female,/path/to/regions.bed

SAMPLE INFORMATION

The sample configuration file defines details of each sample to process:

details:


  - analysis: variant2


    lane: 1


    description: Example1


    files: [in_pair_1.fq, in_pair_2.fq]


    genome_build: hg19


    algorithm:


      platform: illumina


    metadata:


      batch: Batch1


      sex: female


      platform_unit: flowcell-barcode.lane


      library: library_type

UPLOAD

The upload section of the sample configuration file describes where to put the final output files of the pipeline. At its simplest, you can configure bcbio-nextgen to upload results to a local directory, for example a folder shared amongst collaborators or a Dropbox account. You can also configure it to upload results automatically to a Galaxy instance, to Amazon S3 or to iRODS. Here is the simplest configuration, uploading to a local directory:

upload:


  dir: /local/filesystem/directory

General parameters, always required:

Galaxy parameters:

Here is an example configuration for uploading to a Galaxy instance. This assumes you have a shared mounted filesystem that your Galaxy instance can also access:

upload:


  method: galaxy


  dir: /path/to/shared/galaxy/filesystem/folder


  galaxy_url: http://url-to-galaxy-instance


  galaxy_api_key: YOURAPIKEY


  galaxy_library: data_library_to_upload_to

Your Galaxy universe_wsgi.ini configuration needs to have allow_library_path_paste = True set to enable uploads.

S3 parameters:

For S3 access credentials, set the standard environmental variables, AWS_ACCESS_KEY_ID, AWS_SECRET_ACCESS_KEY, and AWS_DEFAULT_REGION or use IAM access roles with an instance profile on EC2 to give your instances permission to create temporary S3 access.

iRODS parameters:

example configuration

Uploads to iRODS depend on a valid installation of the iCommands CLI, and a preconfigured connection through the iinit command.

GLOBALS

You can define files used multiple times in the algorithm section of your configuration in a top level globals dictionary. This saves copying and pasting across the configuration and makes it easier to manually adjust the configuration if inputs change:

globals:


  my_custom_locations: /path/to/file.bed
details:


  - description: sample1


    algorithm:


      variant_regions: my_custom_locations


  - description: sample2


    algorithm:


      variant_regions: my_custom_locations

ALGORITHM PARAMETERS

The YAML configuration file provides a number of hooks to customize analysis in the sample configuration file. Place these under the algorithm keyword.

ploidy:


  default: 2


  mitochondrial: 1


  female: 2


  male: 1

validate: giab-NA12878/truth_small_variants.vcf.gz
validate_regions: giab-NA12878/truth_regions.bed
svvalidate:


  DEL: giab-NA12878/truth_DEL.bed

bcbio_fastq_umi_prep.py autopair -c <cores_to_use> <list> <of> <fastq> <files>

bcbio_fastq_umi_prep.py autopair -c <cores_to_use> --tag1 5 --tag2 5 <list> <of> <fastq> <files>

resources:


  fgbio:


    options: [--min-reads, 2]

bcbiornaseq:


  organism: homo sapiens


  interesting_groups: [treatment, genotype, etc, etc]

resources:


  macs2:


    options: ["--broad"]

resources:


  preseq:


    extrap_fraction: 5


    steps: 500


    seg_len: 5000

resources:


  preseq:


    extrap: 10000000


    step: 30000


    options: ["-D"]

resources:


  multiqc:


    options: ["--cl_config", "'read_count_multiplier: 1'"]

resources:


  multiqc:


    options: ["--cl_config", "'{read_count_multiplier: 0.001, read_count_prefix: K}'"]

variantcaller: [mutect2, varscan, freebayes, vardict]
ensemble:


  numpass: 2


  use_filtered: false

variantcaller: [gatk, freebayes, samtools, gatk-haplotype, varscan]
ensemble:


  format-filters: [DP < 4]


  classifier-params:


    type: svm


  classifiers:


    balance: [AD, FS, Entropy]


    calling: [ReadPosEndDist, PL, PLratio, Entropy, NBQ]


  trusted-pct: 0.65

RESOURCES

The resources section allows customization of locations of programs and memory and compute resources to devote to them:

resources:


  bwa:


    cores: 12


    cmd: /an/alternative/path/to/bwa


  samtools:


    cores: 16


    memory: 2G


  gatk:


    jvm_opts: ["-Xms2g", "-Xmx4g"]

resources:


  sentieon:


    keyfile:


      SENTIEON_LICENSE_SERVER: 100.100.100.100:8888


      SENTIEON_AUTH_MECH: XXX


      SENTIEON_AUTH_DATA: signature

resources:


  tmp:


    dir: /scratch

resources:


  tmp:


    dir: $YOUR_SCRATCH_LOCATION

- description: Example


  analysis: variant2


  resources:


    novoalign:


      options: ["-o", "FullNW", "--rOQ"]


    tmp:


      dir: tmp/sampletmpdir

details:


  - description: Example
resources:


  default:


    cores: 16

log_dir: /path/to/logs

INPUT FILE PREPARATION

Input files for supplementing analysis, like variant_regions need to match the specified reference genome. A common cause of confusion is the two chromosome naming schemes for human genome build 37: UCSC-style in hg19 (chr1, chr2) and Ensembl/NCBI style in GRCh37 (1, 2). To help avoid some of this confusion, in build 38 we only support the commonly agreed on chr1, chr2 style.

It's important to ensure that the chromosome naming in your input files match those in the reference genome selected. bcbio will try to detect this and provide helpful errors if you miss it.

To convert chromosome names, you can use Devon Ryan's collection of chromosome mappings as an input to sed. For instance, to convert hg19 chr-style coordinates to GRCh37:

wget --no-check-certificate -qO- http://raw.githubusercontent.com/dpryan79/ChromosomeMappings/master/GRCh37_UCSC2ensembl.txt \


   | awk '{if($1!=$2) print "s/^"$1"/"$2"/g"}' > remap.sed
sed -f remap.sed original.bed > final.bed

GENOME CONFIGURATION FILES

Each genome build has an associated buildname-resources.yaml configuration file which contains organism specific naming and resource files. bcbio-nextgen expects a resource file present next to the genome FASTA file. Example genome configuration files are available, and automatically installed for natively supported genomes. Create these by hand to support additional organisms or builds.

The major sections of the file are:

By default, we place the buildname-resources.yaml files next to the genome FASTA files in the reference directory. For custom setups, you specify an alternative directory in the ref:config-resources section of your bcbio_system.yaml file:

resources:


  genome:


    dir: /path/to/resources/files

REFERENCE GENOME FILES

The pipeline requires access to reference genomes, including the raw FASTA sequence and pre-built indexes for aligners. The automated installer will install reference files and indexes for commonly used genomes (see the upgrade-install documentation for command line options).

For human genomes, we recommend using build 38 (hg38). This is fully supported and validated in bcbio, and corrects a lot of issues in the previous build 37. We use the 1000 genomes distribution which includes HLAs and decoy sequences. For human build 37, GRCh37 and hg19, we use the 1000 genome references provided in the GATK resource bundle. These differ in chromosome naming: hg19 uses chr1, chr2, chr3 style contigs while GRCh37 uses 1, 2, 3. They also differ slightly in content: GRCh37 has masked Pseudoautosomal regions on chromosome Y allowing alignment to these regions on chromosome X.

You can use pre-existing data and reference indexes by pointing bcbio-nextgen at these resources. We use the Galaxy .loc files approach to describing the location of the sequence and index data, as described in data-requirements. This does not require a Galaxy installation since the installer sets up a minimal set of .loc files. It finds these by identifying the root galaxy directory, in which it expects a tool-data sub-directory with the .loc files. It can do this in two ways:

To manually make genomes available to bcbio-nextgen, edit the individual .loc files with locations to your reference and index genomes. You need to edit sam_fa_indices.loc to point at the FASTA files and then any genome indexes corresponding to aligners you'd like to use (for example: bwa_index.loc for bwa and bowtie2_indices.loc for bowtie2). The database key names used (like GRCh37 and mm10) should match those used in the genome_build of your sample input configuration file.

ADDING CUSTOM GENOMES

bcbio_setup_genome.py will help you to install a custom genome and apply all changes needed to the configuration files. It needs the genome in FASTA format, and the annotation file in GTF or GFF3 format. It can create index for all aligners used by bcbio. Moreover, it will create the folder rnaseq to allow you run the RNAseq pipeline without further configuration.

bcbio_setup_genome.py -f genome.fa -g annotation.gtf -i bowtie2 star seq -n Celegans -b WBcel135

If you want to add smallRNA-seq data files, you will need to add the 3 letters code of mirbase for your genome (i.e hsa for human) and the GTF file for the annotation of smallRNA data. Here you can use the same file than the transcriptome if no other available.

bcbio_setup_genome.py -f genome.fa -g annotation.gtf -i bowtie2 star seq -n Celegans -b WBcel135 --species cel --srna_gtf another_annotation.gtf

To use that genome just need to configure your YAML files as:

genome_build: WBcel135

$ snpEff databases | grep Lactobacillus | grep pentosus
Lactobacillus_pentosus_dsm_20314                                Lactobacillus_pentosus_dsm_20314                                              ENSEMBL_BFMPP_32_179            http://downloads.sourceforge.net/project/snpeff/databases/v4_3/snpEff_v4_3_ENSEMBL_BFMPP_32_179.zip
Lactobacillus_pentosus_kca1                                     Lactobacillus_pentosus_kca1                                                   ENSEMBL_BFMPP_32_179            http://downloads.sourceforge.net/project/snpeff/databases/v4_3/snpEff_v4_3_ENSEMBL_BFMPP_32_179.zip

$ cd /path/to/bcbio/genomes/Lacto/Lactobacillus_pentosus
$ mkdir snpEff
$ cd snpEff
$ wget http://downloads.sourceforge.net/project/snpeff/databases/v4_3/snpEff_v4_3_ENSEMBL_BFMPP_32_179.zip
$ unzip snpEff_v4_3_ENSEMBL_BFMPP_32_179.zip
$ find . -name "Lactobacillus_pentosus_dsm_20314"


 ./home/pcingola/snpEff/data/Lactobacillus_pentosus_dsm_20314
$ mv ./home/pcingola/snpEff/data/Lactobacillus_pentosus_dsm_20314 .

aliases:


  snpeff: Lactobacillus_pentosus_dsm_20314

TUNING CORE AND MEMORY USAGE

bcbio has two ways to specify core usage, helping provide options for parallelizing different types of processes:

resources:


  default:


    memory: 2G


    cores: 16


    jvm_opts: ["-Xms750m", "-Xmx2000m"]

bcbio uses these settings, along with memory requests, to determine how to partition jobs. For example, if you had -n 32 and cores: 16 for a run on a single 32 core machine, this would run two simultaneous bwa mapping jobs using 16 cores each.

Memory specifications (both in memory and jvm_opts) are per-core. bcbio takes care of adjusting this memory to match the cores used. In the example above, if bcbio was running a 16 core java process, it would use 32Gb of memory for the JVM, adjusting Xmx and Xms to match cores used. Internally bcbio looks at the memory and CPU usage on a machine and matches your configuration options to the available system resources. It will scale down core requests if memory is limiting, avoiding over-scheduling resources during the run. You ideally want to set both memory and jvm_opts to match the average memory per core on the run machine and adjust upwards if this does not provide enough memory for some processes during the run.

For single machine runs with a small number of samples, you generally want to set cores close to or equal the number of total cores you're allocating to the job with -n. This will allow individual samples to process as fast as possible and take advantage of multicore software.

For distributed jobs, you want to set cores to match the available cores on a single node in your cluster, then use -n as a multiple of this to determine how many nodes to spin up. For example, cores: 16 and -n 64 would try to make four 16 core machines available for analysis.

MULTIPLE CORES

Running using multiple cores only requires setting the -n command line flag:

bcbio_nextgen.py bcbio_sample.yaml -t local -n 12

IPYTHON PARALLEL

IPython parallel provides a distributed framework for performing parallel computation in standard cluster environments. The bcbio-nextgen setup script installs both IPython and pyzmq, which provides Python bindings for the ZeroMQ messaging library. The only additional requirement is that the work directory where you run the analysis is accessible to all processing nodes. This is typically accomplished with a distributed file system like NFS, Gluster or Lustre.

Run an analysis using ipython for parallel execution:

bcbio_nextgen.py bcbio_sample.yaml -t ipython -n 12 -s lsf -q queue

The -s flag specifies a type of scheduler to use (lsf, sge, torque, slurm, pbspro).

The -q flag specifies the queue to submit jobs to.

The -n flag defines the total number of cores to use on the cluster during processing. The framework will select the appropriate number of cores and type of cluster (single core versus multi-core) to use based on the pipeline stage (see the internals-parallel section in the internals documentation for more details). For multiple core steps, the number of cores to use for programs like bwa, novoalign and gatk comes from the config-resources section of the configuration. Ensure the cores specification matches the physical cores available on machines in your cluster, and the pipeline will divide the total cores specified by -n into the appropriate number of multicore jobs to run.

The pipeline default parameters assume a system with minimal time to obtain processing cores and consistent file system accessibility. These defaults allow the system to fail fast in the case of cluster issues which need diagnosis. For running on shared systems with high resource usage and potential failures due to intermittent cluster issues, there are turning parameters that increase resiliency. The --timeout flag specifies the numbers of minutes to wait for a cluster to start up before timing out. This defaults to 15 minutes. The --retries flag specify the number of times to retry a job on failure. In systems with transient distributed file system hiccups like lock errors or disk availability, this will provide recoverability at the cost of resubmitting jobs that may have failed for reproducible reasons.

Finally, the -r resources flag specifies resource options to pass along to the underlying queue scheduler. This currently supports SGE's -l parameter, Torque's -l parameter and LSF and SLURM native flags. This allows specification or resources to the scheduler (see the qsub man page). You may specify multiple resources, so -r mem=4g -r ct=01:40:00 translates to -l mem=4g -l ct=01:40:00 when passed to qsub or -r "account=a2010002" -r "timelimit=04:00:00" when using SLURM, for instance. SLURM and Torque support specification of an account parameter with -r account=your_name, which IPython transfers into -A.

SGE supports special parameters passed using resources to help handle the heterogeneity of possible setups.

Specify an SGE parallel environment that supports using multiple cores on a single node with -r pename=your_pe. Since this setup is system specific it is hard to write general code for find a suitable environment. Specifically, when there are multiple usable parallel environments, it will select the first one which may not be correct. Manually specifying it with a pename= flag to resources will ensure correct selection of the right environment. If you're administering a grid engine cluster and not sure how to set this up you'd typically want a smp queue using allocation_rule: $pe_slots like in this example pename configuration or smp template.

SGE has other specific flags you may want to tune, depending on your setup. To specify an advanced reservation with the -ar flag, use -r ar=ar_id. To specify an alternative memory management model instead of mem_free use -r memtype=approach. It is further recommended to configure mem_free (or any other chosen memory management model) as a consumable, requestable resource in SGE to prevent overfilling hosts that do not have sufficient memory per slot. This can be done in two steps. First, launch qmon as an admin, select Complex Configuration in qmon, click on mem_free`, under the ``Consumable dialog select JOB (instead of YES or NO) and finally click Modify for the changes to take effect. Secondly, for each host in the queue, configure mem_free as a complex value. If a host called myngshost has 128GB of RAM, the corresponding command would be qconf -mattr exechost complex_values mem_free=128G myngshost

There are also special -r resources parameters to support pipeline configuration:

TROUBLESHOOTING

bcbio_nextgen.py your_input.yaml -n 1

python -c 'import socket; print socket.gethostbyname(socket.gethostname())'

host-ip hostname

MEMORY MANAGEMENT

The memory information specified in the system configuration config-resources enables scheduling of memory intensive processes. The values are specified on a memory-per-core basis and thus bcbio-nextgen handles memory scheduling by:

As a result of these calculations, the cores used during processing will not always correspond to the maximum cores provided in the input -n parameter. The goal is rather to intelligently maximize cores and memory while staying within system resources. Note that memory specifications are for a single core, and the pipeline takes care of adjusting this to actual cores used during processing.

DETERMINING AVAILABLE CORES AND MEMORY PER MACHINE

bcbio automatically tries to determine the total available memory and cores per machine for balancing resource usage. For multicore runs, it retrieves total memory from the current machine. For parallel runs, it spawns a job on the queue and extracts the system information from that machine. This expects a homogeneous set of machines within a cluster queue. You can see the determined cores and total memory in provenance/system-ipython-queue.yaml.

For heterogeneous clusters or other cases where bcbio does not correctly identify available system resources, you can manually set the machine cores and total memory in the resource section of either a sample YAML file (sample-resources) or bcbio_system.yaml:

resources:


  machine:


    memory: 48.0


    cores: 16

The memory usage is total available on the machine in Gb, so this specifies that individual machines have 48Gb of total memory and 16 cores.

TUNING SYSTEMS FOR SCALE

bcbio-nextgen scales out on clusters including hundreds of cores and is stress tested on systems with 1000 simultaneous processes. Scaling up often requires system specific tuning to handle simultaneous processes. This section collects useful tips and tricks for managing scaling issues.

execd_params                 S_DESCRIPTORS=20000

CURRENT STATUS

bcbio creates CWL for alignment, small variant calls (SNPs and indels), coverage assessment, HLA typing, quality control and structural variant calling. It generates a CWL v1.0.2 compatible workflow. The actual biological code execution during runs works with either a bcbio docker container or a local installation of bcbio.

The implementation includes bcbio's approaches to splitting and batching analyses. At the top level workflow, we parallelize by samples. Using sub-workflows, we split fastq inputs into sections for parallel alignment over multiple machines following by merging. We also use sub-workflows, along with CWL records, to batch multiple samples and run in parallel. This enables pooled and tumor/normal cancer calling with parallelization by chromosome regions based on coverage calculations.

bcbio supports these CWL-compatible tools:

We plan to continue to expand CWL support to include more components of bcbio, and also need to evaluate the workflow on larger, real life analyses. This includes supporting additional CWL runners. We're working on evaluating Galaxy/Planemo for integration with the Galaxy community.

INSTALLATION

bcbio-vm installs all dependencies required to generate CWL and run bcbio, along with supported CWL runners. There are two install choices, depending on your usage of bcbio: running CWL with a existing local bcbio install, or running with containers.

bcbio_nextgen.py upgrade --cwl

wget http://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh
export TARGETDIR=~/install/bcbio-vm/anaconda
export BINDIR=/usr/local
bash Miniconda2-latest-Linux-x86_64.sh -b -p $TARGETDIR
$TARGETDIR/bin/conda install --yes -c conda-forge -c bioconda bcbio-nextgen
$TARGETDIR/bin/conda install --yes -c conda-forge -c bioconda bcbio-nextgen-vm
ln -s $TARGETDIR/bin/bcbio_vm.py $BINDIR/bcbio_vm.py
ln -s $TARGETDIR/bin/conda $BINDIR/bcbiovm_conda
ln -s $TARGETDIR/bin/python $BINDIR/bcbiovm_python

export PATH=$TARGETDIR/bin:$PATH

GETTING STARTED

To make it easy to get started, we have pre-built CWL descriptions that use test data. These run in under 5 minutes on a local machine and don't require a bcbio installation if you have Docker available on your machine:

wget -O test_bcbio_cwl.tar.gz https://github.com/bcbio/test_bcbio_cwl/archive/master.tar.gz
tar -xzvpf test_bcbio_cwl.tar.gz
cd test_bcbio_cwl-master/somatic

bash run_cromwell.sh
bash run_bunny.sh
bash run_toil.sh

bcbio_vm.py cwlrun cromwell somatic-workflow
bcbio_vm.py cwlrun toil somatic-workflow
bcbio_vm.py cwlrun bunny somatic-workflow

GENERATING CWL FOR INPUT TO A TOOL

The first step in running your analysis project in bcbio is to generate CWL. If you're already familiar with bcbio, the process of preparing information about your sample inputs and analysis are almost identical:

local:


  ref: /path/to/bcbio/genomes/Hsapiens


  inputs:


    - /path/to/input/files
resources:


  default:


    cores: 16


    memory: 3500M


    jvm_opts: [-Xms1g, -Xmx3500m]

Generate CWL with:

bcbio_vm.py template --systemconfig bcbio_system.yaml template.yaml samples.csv [optional list of fastq or BAM inputs]
bcbio_vm.py cwl --systemconfig bcbio_system.yaml samples/config/samples.yaml

producing a sample-workflow output directory with the CWL.

On a first CWL generation run with a new genome, this process will run for a longer time as it needs to make your reference compatible with CWL. This includes creating single tar.gz files from some reference directories so they can get passed to CWL steps where they'll get unpacked. This process only happens a single time and keeps unpacked versions so your reference setup is compatible with both old bcbio IPython and new CWL runs.

You can now run this with any CWL compatible runner and the bcbio_vm.py cwlrun wrappers standardize running across multiple tools in different environments.

RUNNING WITH CROMWELL (LOCAL, HPC, CLOUD)

The Cromwell workflow management system runs bcbio either locally on a single machine or distributed on a cluster using a scheduler like SLURM, SGE or PBSPro.

To run a bcbio CWL workflow locally using Docker:

bcbio_vm.py cwlrun cromwell sample-workflow

If you want to run from a locally installed bcbio add --no-container to the commandline.

To run distributed on a SLURM cluster:

bcbio_vm.py cwlrun cromwell sample-workflow --no-container -q your_queue -s slurm -r timelimit=0-12:00

Tweak scheduler parameters using the same options as the older bcbio IPython approach.

To control the resources used Cromwell, set --joblimit to the allowed jobs allocated concurrently. This isn't total cores used, but rather the number of jobs either locally or remotely scheduled concurrently. Since CWL steps are heterogeneous and use only cores necessary for that job, the total cores used will max out at joblimit times maximum cores for an individual process. Setting this helps avoid over-committing jobs to a shared scheduler during highly parallel processes like variant calling.

Cromwell can also run directly on cloud resources: docs-cloud-gcp.

RUNNING WITH TOIL (LOCAL, HPC)

The Toil pipeline management system runs CWL workflows in parallel on a local machine, on a cluster or at AWS.

To run a bcbio CWL workflow locally with Toil using Docker:

bcbio_vm.py cwlrun toil sample-workflow

If you want to run from a locally installed bcbio add --no-container to the commandline.

To run distributed on a Slurm cluster:

bcbio_vm.py cwlrun toil sample-workflow -- --batchSystem slurm

RUNNING ON ARVADOS (HOSTED CLOUD)

bcbio generated CWL workflows run on Arvados and these instructions detail how to run on the Arvdos public instance. Arvados cwl-runner comes pre-installed with bcbio-vm. We have a publicly accessible project, called bcbio_resources that contains the latest Docker images, test data and genome references you can use for runs.

Retrieve API keys from the Arvados public instance. Login, then go to 'User Icon -> Personal Token'. Copy and paste the commands given there into your shell. You'll specifically need to set ARVADOS_API_HOST and ARVADOS_API_TOKEN.

To run an analysis:

arv-put --name testdata_genomes --project-uuid $PROJECT_ID testdata/genomes/hg19

arv-put --name testdata_inputs --project-uuid $PROJECT_ID testdata/100326_FC6107FAAXX testdata/automated testdata/reference_material

arvados:


  reference: qr1hi-4zz18-kuz1izsj3wkfisq


  input: [qr1hi-j7d0g-h691y6104tlg8b4]
resources:


  default: {cores: 4, memory: 2G, jvm_opts: [-Xms750m, -Xmx2500m]}

bcbio_vm.py template --systemconfig bcbio_system_arvados.yaml testcwl_template.yaml testcwl.csv

bcbio_vm.py cwl --systemconfig bcbio_system_arvados.yaml testcwl/config/testcwl.yaml

arv-copy $JOBS_ID --project-uuid $PROJECT_ID --src qr1hi --dst qr1hi
arv-copy $BCBIO_VC_ID --project-uuid $PROJECT_ID --src qr1hi --dst qr1hi

docker pull arvados/jobs:1.0.20180216164101
arv-keepdocker --project $PROJECT_ID -- arvados/jobs 1.0.20180216164101
docker pull quay.io/bcbio/bcbio-vc
arv-keepdocker --project $PROJECT_ID -- quay.io/bcbio/bcbio-vc latest

bcbio_vm.py cwlrun arvados arvados_testcwl-workflow -- --project-uuid $PROJECT_ID

RUNNING ON DNANEXUS (HOSTED CLOUD)

bcbio runs on the DNAnexus platform by converting bcbio generated CWL into DNAnexus workflows and apps using dx-cwl. This describes the process using the bcbio workflow app (bcbio-run-workflow) and bcbio workflow applet (bcbio_resources:/applets/bcbio-run-workflow) in the public bcbio_resources project, Both are regularly updated and maintained on the DNAnexus platform. Secondarily, we also show how to install and create workflows locally for additional control and debugging.

dx new project $PNAME

dx select $PNAME
dx upload -p --path /data/input *.bam

dnanexus:


  project: PNAME


  ref:


    project: bcbio_resources


    folder: /reference_genomes


  inputs:


    - /data/input


    - /data/input/regions
resources:


  default: {cores: 8, memory: 3000M, jvm_opts: [-Xms1g, -Xmx3000m]}

samplename,description,batch,phenotype
file1.bam,sample1,b1,tumor
file2.bam,sample2,b1,normal
file3.bam,sample3,b2,tumor
file4.bam,sample4,b2,normal

details:


 - algorithm:


     aligner: bwa


     variantcaller: gatk-haplotype


   analysis: variant2


   genome_build: hg38

TEMPLATE=germline
APP_VERSION=0.0.2
FOLDER=/bcbio/$PNAME
dx select "$PROJECT"
dx mkdir -p $FOLDER
for F in $TEMPLATE-template.yaml $PNAME.csv bcbio_system-dnanexus.yaml
do


        dx rm -a /$FOLDER/$F || true


        dx upload --path /$FOLDER/ $F
done
dx ls $FOLDER
dx rm -a -r /$FOLDER/dx-cwl-run || true
dx run bcbio-run-workflow/$APP_VERSION -iyaml_template=/$FOLDER/$TEMPLATE-template.yaml -isample_spec=/$FOLDER/$PNAME.csv -isystem_configuration=/$FOLDER/bcbio_system-dnanexus.yaml -ioutput_folder=/$FOLDER/dx-cwl-run

dx run bcbio_resources:/applets/bcbio-run-workflow -iyaml_template=/$FOLDER/$TEMPLATE-template.yaml -isample_spec=/$FOLDER/$PNAME.csv -isystem_configuration=/$FOLDER/bcbio_system-dnanexus.yaml -ioutput_folder=/$FOLDER/dx-cwl-run

The app will lookup all files, prepare a bcbio CWL workflow, convert into a DNAnexus workflow, and submit to the platform. The workflow runs as a standard DNAnexus workflow and you can monitor through the command line (with dx find executions --root job-YOURJOBID and dx watch) or the web interface (Monitor tab).

If you prefer not to use the DNAnexus app, you can also submit jobs locally by installing bcbio-vm on your local machine. This can also be useful to test generation of CWL and manually ensure identification of all your samples and associated files on the DNAnexus platform.

TEMPLATE=germline
rm -rf $PNAME $PNAME-workflow
bcbio_vm.py template --systemconfig bcbio_system-dnanexus.yaml $TEMPLATE-template.yaml $PNAME.csv
bcbio_vm.py cwl --systemconfig bcbio_system-dnanexus.yaml $PNAME/config/$PNAME.yaml

dx env

dx-cwl compile-workflow $PNAME-workflow/main-$PNAME.cwl \


   --project PROJECT_ID --token $DX_AUTH_TOKEN \


   --rootdir $FOLDER/dx-cwl-run

FOLDER=/bcbio/$PNAME
dx mkdir -p $DX_PROJECT_ID:$FOLDER/$PNAME-workflow
dx upload -p --path $DX_PROJECT_ID:$FOLDER/$PNAME-workflow $PNAME-workflow/main-$PNAME-samples.json
dx-cwl run-workflow $FOLDER/dx-cwl-run/main-$PNAME/main-$PNAME \


       $FOLDER/$PNAME-workflow/main-$PNAME-samples.json \


       --project PROJECT_ID --token $DX_AUTH_TOKEN \


       --rootdir $FOLDER/dx-cwl-run

DEVELOPMENT NOTES

bcbio generates a common workflow language description. Internally, bcbio represents the files and information related to processing as a comprehensive dictionary. This world object describes the state of a run and associated files, and new processing steps update or add information to it. The world object is roughly equivalent to CWL's JSON-based input object, but CWL enforces additional annotations to identify files and models new inputs/outputs at each step. The work in bcbio is to move from our laissez-faire approach to the more structured CWL model.

The generated CWL workflow is in run_info-cwl-workflow:

To help with defining the outputs at each step, there is a WorldWatcher object that can output changed files and world dictionary objects between steps in the pipeline when running a bcbio in the standard way. The variant pipeline has examples using it. This is useful when preparing the CWL definitions of inputs and outputs for new steps in the bcbio CWL step definitions.

TODO

bcbio has two approaches to running on cloud providers like Amazon Web Services (AWS), Google Cloud (GCP) and Microsoft Azure. For smaller projects we use a simplified ansible based approach which automates spinning up single multicore machines for running either traditional or docs-cwl bcbio runs.

For larger distributed projects, we're actively working on using docs-cwl support with runners like Cromwell that directly interface and run on cloud services. We'll document these approaches here as they're tested and available.

For getting started, the CWL docs-cwl-installation documentation describes how to install bcbio-vm, which provides a wrapper around bcbio that automates interaction with cloud providers and Docker. bcbio_vm.py also cleans up the command line usage to make it more intuitive and provides a superset of functionality available in bcbio_nextgen.py.

GOOGLE CLOUD (GCP)

Cromwell runs bcbio CWL pipelines on Google Cloud using the Google Pipelines API.

bcbio_conda install -c conda-forge -c bioconda google-cloud-sdk

gcloud auth login
gcloud projects create your-project
gcloud iam service-accounts create your-service-account
gcloud projects add-iam-policy-binding your-project --member \


  "serviceAccount:your-service-account@your-project.iam.gserviceaccount.com" --role "roles/owner"
gcloud iam service-accounts keys create ~/.config/glcoud/your-service-account.json \


  --iam-account your-service-account@your-project.iam.gserviceaccount.com
export GOOGLE_APPLICATION_CREDENTIALS=~/.config/gcloud/your-service-account.json

gcloud config set project your-project
gsutil mb gs://your-project

gsutil cp your_data.bam gs://your-project/inputs/

gs:


  ref: gs://bcbiodata/collections


  inputs:


    - gs://your-project/inputs
resources:


  default: {cores: 2, memory: 3G, jvm_opts: [-Xms750m, -Xmx3000m]}

bcbio_vm.py template --systemconfig bcbio_system-gcp.yaml ${TEMPLATE}-template.yaml $PNAME.csv
bcbio_vm.py cwl --systemconfig bcbio_system-gcp.yaml $PNAME/config/$PNAME.yaml

bcbio_vm.py cwlrun cromwell $PNAME --cloud-project your-project \


    --cloud-root gs://your-project/work_cromwell

AMAZON WEB SERVICES (OLD)

Amazon Web Services (AWS) provides a flexible cloud based environment for running analyses. Cloud approaches offer the ability to perform analyses at scale with no investment in local hardware. They also offer full programmatic control over the environment, allowing bcbio to automate the entire setup, run and teardown process.

bcbio-vm uses Elasticluster to build a cluster on AWS with an encrypted NFS mounted drive and an optional Lustre shared filesystem. We're phasing out this approach to cloud support in bcbio and will be actively moving to Common Workflow Language based approaches.

bcbio_vm.py template s3://your-project/your-analysis/template.yaml s3://your-project/your-analysis/name.csv

resources:


  gatk:


    jar: s3://bcbio-syn3-eval/jars/GenomeAnalysisTK.jar

export AWS_ACCESS_KEY_ID=your_access_key
export AWS_SECRET_ACCESS_KEY=your_secret_key

bcbio_vm.py aws iam --region=us-east-1

bcbio_vm.py aws vpc --region=us-east-1

bcbio_vm.py aws config edit

bcbio_vm.py aws info

bcbio_vm.py aws cluster start

bcbio_vm.py aws icel create

bcbio_vm.py aws icel mount

bcbio_vm.py aws cluster ssh

mkdir /encrypted/your-project
cd !$ && mkdir work && cd work

sudo mkdir /scratch/cancer-dream-syn3-exome
sudo chown ubuntu !$
cd !$ && mkdir work && cd work

bcbio_vm.py run -n 8 s3://your-project/your-analysis/name.yaml

bcbio_vm.py ipythonprep s3://your-project/your-analysis/name.yaml slurm cloud -n 60
sbatch bcbio_submit.sh

bcbio_vm.py graph bcbio-nextgen.log

``

`

bcbio_vm.py aws cluster stop

bcbio_vm.py aws icel stop

COMMON FILES

The output directory contains sample specific output files labeled by sample name and a more general project directory. The sample directories contain all of the sample specific output files, while the project directory contains global files like project summaries or batched population level variant calls. See the teaching documentation for a full variant calling example with additional details about configuration setting and resulting output files.

VARIANT CALLING

RNA-SEQ

SMALL RNA-SEQ

DOWNSTREAM ANALYSIS

This section collects useful scripts and tools to do downstream analysis of bcbio-nextgen outputs. If you have pointers to useful tools, please add them to the documentation.

This section provides useful concepts for getting started digging into the code and contributing new functionality. We welcome contributors and hope these notes help make it easier to get started.

DEVELOPMENT GOALS

During development we seek to maximize functionality and usefulness, while avoiding complexity. Since these goals are sometimes in conflict, it's useful to understand the design approaches:

OVERVIEW

The most useful modules inside bcbio, ordered by likely interest:

DEVELOPMENT INFRASTRUCTURE

bcbio-nextgen uses GitHub for code development, and we welcome pull requests. GitHub makes it easy to establish custom forks of the code and contribute those back. The Biopython documentation has great information on using git and GitHub for a community developed project. In short, make a fork of the bcbio code by clicking the Fork button in the upper right corner of the GitHub page, commit your changes to this custom fork and keep it up to date with the main bcbio repository as you develop. The github help pages have detailed information on keeping your fork updated with the main github repository (e.g. https://help.github.com/articles/syncing-a-fork/). After commiting changes, click New Pull Request from your fork when you'd like to submit your changes for integration in bcbio.

For developing and testing changes locally, you can install directly into a bcbio-nextgen installation. The automated bcbio-nextgen installer creates an isolated Python environment using Anaconda. This will be a subdirectory of your installation root, like /usr/local/share/bcbio-nextgen/anaconda. The installer also includes a bcbio_python executable target which is the python in this isolated anaconda directory. You generally will want to make changes to your local copy of the bcbio-nextgen code and then install these into the code directory. To install from your bcbio-nextgen source tree for testing do:

bcbio_python setup.py install

One tricky part that we don't yet know how to work around is that pip and standard setup.py install have different ideas about how to write Python eggs. setup.py install will create an isolated python egg directory like bcbio_nextgen-0.7.5a-py2.7.egg, while pip creates an egg pointing to a top level bcbio directory. Where this gets tricky is that the top level bcbio directory takes precedence. The best way to work around this problem is to manually remove the current pip installed bcbio-nextgen code (rm -rf /path/to/anaconda/lib/python2.7/site-packages/bcbio*) before managing it manually with bcbio_python setup.py install. We'd welcome tips about ways to force consistent installation across methods.

If you want to test with bcbio_nextgen code in a separate environment from your work directory, we recommend using the installer to install only the bcbio code into a separate directory:

python bcbio_nextgen_install.py /path/to/testbcbio --nodata --isolate

Then add this directory to your PATH before your bcbio installation with the tools: export PATH=/path/to/testbcbio/anaconda/bin:$PATH, or directly calling the testing bcbio /path/to/testbcbio/anaconda/bin/bcbio_nextgen.py.

BUILDING THE DOCUMENTATION LOCALLY

If you have added or modified this documentation, to build it locally and see how it looks like you can do so by running:

cd docs
make html

The documentation will be built under docs/_build/html, open index.html with your browser to load your local build.

If you want to use the same theme that Read The Docs uses, you can do so by installing sphinx_rtd_theme via pip. You will also need to add this in the docs/conf.py file to use the theme only locally:

html_theme = 'default'
on_rtd = os.environ.get('READTHEDOCS', False)
if not on_rtd:  # only import and set the theme if we're building docs locally


    import sphinx_rtd_theme


    html_theme = 'sphinx_rtd_theme'


    html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]

ADDING TOOLS

ADDING NEW ORGANISMS

While bcbio-nextgen and supporting tools receive the most testing and development on human or human-like diploid organisms, the algorithms are generic and we strive to support the wide diversity of organisms used in your research. We welcome contributors interested in setting up and maintaining support for their particular research organism, and this section defines the steps in integrating a new genome. We also welcome suggestions and implementations that improve this process.

Setup CloudBioLinux to automatically download and prepare the genome:

Add the organism to the supported installs within bcbio:

Test installation of genomes by pointing to your local cloudbiolinux edits during a data installation:

mkdir -p tmpbcbio-install
ln -s ~/bio/cloudbiolinux tmpbcbio-install
bcbio_nextgen.py upgrade --data --genomes DBKEY

Add configuration information to bcbio-nextgen by creating a config/genomes/DBKEY-resources.yaml file. Copy an existing minimal template like canFam3 and edit with pointers to snpEff and other genome resources. The VEP database directory has Ensembl names. SnpEff has a command to list available databases:

snpEff databases

Finally, send pull requests for CloudBioLinux and bcbio-nextgen and we'll happily integrate the new genome.

This will provide basic integration with bcbio and allow running a minimal pipeline with alignment and quality control. We also have utility scripts in CloudBioLinux to help with preparing dbSNP (utils/prepare_dbsnp.py) and RNA-seq (utils/prepare_tx_gff.py) resources for some genomes. For instance, to prepare RNA-seq transcripts for mm9:

bcbio_python prepare_tx_gff.py --genome-dir /path/to/bcbio/genomes Mmusculus mm9

We are still working on ways to best include these as part of the standard build and install since they either require additional tools to run locally, or require preparing copies in S3 buckets.

STANDARD FUNCTION ARGUMENTS

{'lane': '7_100326_FC6107FAAXX',


 'pl': 'illumina',


 'pu': '7_100326_FC6107FAAXX',


 'rg': '7',


 'sample': 'Test1'}

{'config': {'algorithm': {'aligner': 'bwa',


                          'callable_regions': 'analysis_blocks.bed',


                          'coverage_depth': 'low',


                          'coverage_interval': 'regional',


                          'mark_duplicates': 'samtools',


                          'nomap_split_size': 50,


                          'nomap_split_targets': 20,


                          'num_cores': 1,


                          'platform': 'illumina',


                          'quality_format': 'Standard',


                          'realign': 'gkno',


                          'recalibrate': 'gatk',


                          'save_diskspace': True,


                          'upload_fastq': False,


                          'validate': '../reference_material/7_100326_FC6107FAAXX-grade.vcf',


                          'variant_regions': '../data/automated/variant_regions-bam.bed',


                          'variantcaller': 'freebayes'},


            'resources': {'bcbio_variation': {'dir': '/usr/share/java/bcbio_variation'},


                          'bowtie': {'cores': None},


                          'bwa': {'cores': 4},


                          'cortex': {'dir': '~/install/CORTEX_release_v1.0.5.14'},


                          'cram': {'dir': '/usr/share/java/cram'},


                          'gatk': {'cores': 2,


                                   'dir': '/usr/share/java/gatk',


                                   'jvm_opts': ['-Xms750m', '-Xmx2000m'],


                                   'version': '2.4-9-g532efad'},


                          'gemini': {'cores': 4},


                          'novoalign': {'cores': 4,


                                        'memory': '4G',


                                        'options': ['-o', 'FullNW']},


                          'picard': {'cores': 1,


                                     'dir': '/usr/share/java/picard'},


                          'snpEff': {'dir': '/usr/share/java/snpeff',


                                     'jvm_opts': ['-Xms750m', '-Xmx3g']},


                          'stampy': {'dir': '~/install/stampy-1.0.18'},


                          'tophat': {'cores': None},


                          'varscan': {'dir': '/usr/share/java/varscan'},


                          'vcftools': {'dir': '~/install/vcftools_0.1.9'}}},
'genome_resources': {'aliases': {'ensembl': 'human',


                                  'human': True,


                                  'snpeff': 'hg19'},


                      'rnaseq': {'transcripts': '/path/to/rnaseq/ref-transcripts.gtf',


                                 'transcripts_mask': '/path/to/rnaseq/ref-transcripts-mask.gtf'},


                      'variation': {'dbsnp': '/path/to/variation/dbsnp_132.vcf',


                                    'train_1000g_omni': '/path/to/variation/1000G_omni2.5.vcf',


                                    'train_hapmap': '/path/to/hg19/variation/hapmap_3.3.vcf',


                                    'train_indels': '/path/to/variation/Mills_Devine_2hit.indels.vcf'},


                      'version': 1},


 'dirs': {'fastq': 'input fastq directory',


              'galaxy': 'directory with galaxy loc and other files',


              'work': 'base work directory'},


 'metadata': {'batch': 'TestBatch1'},


 'genome_build': 'hg19',


 'name': ('', 'Test1'),


 'sam_ref': '/path/to/hg19.fa'}

{'prep_recal': 'Test1/7_100326_FC6107FAAXX-sort.grp',


 'summary': {'metrics': [('Reference organism', 'hg19', ''),


                         ('Total', '39,172', '76bp paired'),


                         ('Aligned', '39,161', '(100.0\\%)'),


                         ('Pairs aligned', '39,150', '(99.9\\%)'),


                         ('Pair duplicates', '0', '(0.0\\%)'),


                         ('Insert size', '152.2', '+/- 31.4')],


             'pdf': '7_100326_FC6107FAAXX-sort-prep-summary.pdf',


             'project': 'project-summary.yaml'},


 'validate': {'concordant': 'Test1-ref-eval-concordance.vcf',


              'discordant': 'Test1-eval-ref-discordance-annotate.vcf',


              'grading': 'validate-grading.yaml',


              'summary': 'validate-summary.csv'},


 'variants': [{'population': {'db': 'gemini/TestBatch1-freebayes.db',


                              'vcf': None},


               'validate': None,


               'variantcaller': 'freebayes',


               'vrn_file': '7_100326_FC6107FAAXX-sort-variants-gatkann-filter-effects.vcf'}],


 'vrn_file': '7_100326_FC6107FAAXX-sort-variants-gatkann-filter-effects.vcf',


 'work_bam': '7_100326_FC6107FAAXX-sort-prep.bam'}

PARALLELIZATION FRAMEWORK

bcbio-nextgen supports parallel runs on local machines using multiple cores and distributed on a cluster using IPython using a general framework.

The first parallelization step starts up a set of resources for processing. On a cluster this spawns a IPython parallel controller and set of engines for processing. The prun (parallel run) start function is the entry point to spawning the cluster and the main argument is a parallel dictionary which contains arguments to the engine processing command. Here is an example input from an IPython parallel run:

{'cores': 12,


 'type': 'ipython'


 'progs': ['aligner', 'gatk'],


 'ensure_mem': {'star': 30, 'tophat': 8, 'tophat2': 8},


 'module': 'bcbio.distributed',


 'queue': 'batch',


 'scheduler': 'torque',


 'resources': [],


 'retries': 0,


 'tag': '',


 'timeout': 15}

The cores and type arguments must be present, identifying the total cores to use and type of processing, respectively. Following that are arguments to help identify the resources to use. progs specifies the programs used, here the aligner, which bcbio looks up from the input sample file, and gatk. ensure_mem is an optional argument that specifies minimum memory requirements to programs if used in the workflow. The remaining arguments are all specific to IPython to help it spin up engines on the appropriate computing cluster.

A shared component of all processing runs is the identification of used programs from the progs argument. The run creation process looks up required memory and CPU resources for each program from the config-resources section of your bcbio_system.yaml file. It combines these resources into required memory and cores using the logic described in the memory-management section of the parallel documentation. Passing these requirements to the cluster creation process ensures the available machines match program requirements.

bcbio-nextgen's pipeline.main code contains examples of starting and using set of available processing engines. This example starts up machines that use samtools, gatk and cufflinks then runs an RNA-seq expression analysis:

with prun.start(_wprogs(parallel, ["samtools", "gatk", "cufflinks"]),


                samples, config, dirs, "rnaseqcount") as run_parallel:


    samples = rnaseq.estimate_expression(samples, run_parallel)

The pipelines often reuse a single set of machines for multiple distributed functions to avoid the overhead of starting up and tearing down machines and clusters.

The run_parallel function returned from the prun.start function enables running on jobs in the parallel on the created machines. The ipython wrapper code contains examples of implementing this. It is a simple function that takes two arguments, the name of the function to run and a set of multiple arguments to pass to that function:

def run(fn_name, items):

The items arguments need to be strings, lists and dictionaries to allow serialization to JSON format. The internals of the run function take care of running all of the code in parallel and returning the results back to the caller function.

In this setup, the main processing code is fully independent from the parallel method used so running on a single multicore machine or in parallel on a cluster return identical results and require no changes to the logical code defining the pipeline.

During re-runs, we avoid the expense of spinning up processing clusters for completed tasks using simple checkpoint files in the checkpoints_parallel directory. The prun.start wrapper writes these on completion of processing for a group of tasks with the same parallel architecture, and on subsequent runs will go through these on the local machine instead of parallelizing. The processing code supports these quick re-runs by checking for and avoiding re-running of tasks when it finds output files.

Plugging new parallelization approaches into this framework involves writing interface code that handles the two steps. First, create a cluster of ready to run machines given the parallel function with expected core and memory utilization:

Second, implement a run_parallel function that handles using these resources to distribute jobs and return results. The multicore wrapper and ipython wrapper are useful starting points for understanding the current implementations.

OVERVIEW

PARALLEL

bcbio calculates callable regions following alignment using goleft depth. These regions determine breakpoints for analysis, allowing parallelization by genomic regions during variant calling. Defining non-callable regions allows bcbio to select breakpoints for parallelization within chromosomes where we won't accidentally impact small variant detection. The callable regions also supplement the variant calls to define positions where not called bases are homozygous reference, as opposed to true no-calls with no evidence. The callable regions plus variant calls is an alternative to gVCF output which explicitly enumerates reference calls in the output variant file.

ABSTRACT

Community Development of Validated Variant Calling Pipelines

Brad Chapman, Rory Kirchner, Oliver Hofmann and Winston Hide Harvard School of Public Health, Bioinformatics Core, Boston, MA, 02115

Translational research relies on accurate identification of genomic variants. However, rapidly changing best practice approaches in alignment and variant calling, coupled with large data sizes, make it a challenge to create reliable and reproducible variant calls. Coordinated community development can help overcome these challenges by sharing testing and updates across multiple groups. We describe bcbio-nextgen, a distributed multi-architecture pipeline that automates variant calling, validation and organization of results for query and visualization. It creates an easily installable, reliable infrastructure from best-practice open source tools with the following goals:

LINKS FROM THE PRESENTATION

This is a teaching orientated example of using bcbio from the Cold Spring Harbor Laboratory's Advanced Sequencing Technology and Applications course. This uses cancer tumor normal data from the ICGC-TCGA DREAM synthetic 3 challenge, subset to exomes on chromosome 6 to reduce runtimes. It demonstrates:

LOADING PRE-RUN ANALYSIS

To save downloading the genome data and running the analysis, we have a pre-prepared AMI with the data and analysis run. Use the AWS Console to launch the pre-built AMI -- search Community AMIs for ami-5e84fe34. Any small instance type is fine for exploring the configuration, run directory and output files. Make sure you associate a public IP and a security group that allows remote ssh.

Once launched, ssh into the remote machine with ssh -i your-keypair ubuntu@public.ip.address to explore the inputs and outputs. The default PATH contains bcbio and third party programs in /usr/local/bin, with the biological data installed in /usr/local/share/bcbio. The run is in a ~/run/cancer-syn3-chr6.

INPUT CONFIGURATION FILE

To run bcbio, you prepare a small configuration file describing your analysis. You can prepare it manually or use an automated configuration method. The example has a pre-written configuration file with tumor/normal data located in the

``

config` directory and this section walks through the settings.

You define the type of analysis (variant calling) along with the input files and genome build:

analysis: variant2
files: [../input/cancer-syn3-chr6-tumor-1.fq.gz, ../input/cancer-syn3-chr6-tumor-2.fq.gz]
genome_build: hg38

Sample description and assignment as a tumor sample, called together with a matched normal:

description: syn3-tumor
metadata:


  batch: syn3


  phenotype: tumor


  sex: female

Next it defines parameters for running the analysis. First we pick our aligner (bwa mem):

algorithm:


  aligner: bwa

Post-alignment, we mark duplicates but do not perform recalibration and realignment:

mark_duplicates: true
recalibrate: false
realign: false

We call variants in exome regions on chromosome 6 using a BED file input, call variants as low as 2% in the tumor sample, and use 3 variant callers along with an ensemble method that combines results for any found in 2 out of 3:

variant_regions: ../input/NGv3-chr6-hg38.bed
min_allele_fraction: 2
variantcaller: [vardict, freebayes, varscan]
ensemble:


  numpass: 2

For structural variant calling, we use two callers and prioritize variants to those found in the CIViC database:

svcaller: [lumpy, manta]
svprioritize: cancer/civic

Call HLA types with OptiType:

hlacaller: optitype

Finally, we validate both the small variants and structural variants. These use pre-installed validation sets that come with bcbio. We limit validation regions to avoid low complexity regions, which cause bias in

``

exclude_regions: [lcr]
validate: dream-syn3-crossmap/truth_small_variants.vcf.gz
validate_regions: dream-syn3-crossmap/truth_regions.bed
svvalidate:


  DEL: dream-syn3-crossmap/truth_DEL.bed


  DUP: dream-syn3-crossmap/truth_DUP.bed


  INV: dream-syn3-crossmap/truth_INV.bed

OUTPUT FILES

Output files are in ~/run/cancer-syn3-chr6/final, extracted from the full work directory in ~/run/cancer-syn3-chr6/work.

The directories with sample information are in syn3-tumor/. Aligned BAMs include a -ready.bam file with all of the original reads (including split and discordants) and separate files with only the split (-sr.bam) and discordant (-disc.bam) reads:

syn3-tumor-ready.bam
syn3-tumor-ready.bam.bai
syn3-tumor-sr.bam
syn3-tumor-sr.bam.bai
syn3-tumor-disc.bam
syn3-tumor-disc.bam.bai

SNP and indel calls for 3 callers, plus combined ensemble calls:

syn3-tumor-ensemble.vcf.gz
syn3-tumor-ensemble.vcf.gz.tbi
syn3-tumor-freebayes.vcf.gz
syn3-tumor-freebayes.vcf.gz.tbi
syn3-tumor-varscan.vcf.gz
syn3-tumor-varscan.vcf.gz.tbi
syn3-tumor-vardict.vcf.gz
syn3-tumor-vardict.vcf.gz.tbi

Structural variant calls for 2 callers, plus a simplified list of structural variants in cancer genes of interest:

syn3-tumor-sv-prioritize.tsv
syn3-tumor-lumpy.vcf.gz
syn3-tumor-lumpy.vcf.gz.tbi
syn3-tumor-manta.vcf.gz
syn3-tumor-manta.vcf.gz.tbi

HLA typing results:

syn3-tumor-hla-optitype.csv

Validation results from comparisons against truth set, including plots:

syn3-tumor-sv-validate.csv
syn3-tumor-sv-validate-DEL.png
syn3-tumor-sv-validate-df.csv
syn3-tumor-sv-validate-DUP.png
syn3-tumor-sv-validate-INV.png
syn3-tumor-validate.png

The top level directory for the project, 2015-11-18_syn3-cshl/ has files relevant to the entire run. There is a consolidated quality control report:

multiqc/multiqc_report.html

Povenance information, with log files of all commands run and program versions used:

bcbio-nextgen.log
bcbio-nextgen-commands.log
programs.txt
data_versions.csv

A top level summary of metrics for alignment, variant calling and coverage that is useful downstream:

project-summary.yaml

PREPARING AND RUNNING

The steps to prepare an AMI from a bare machine and run the analysis. These are pre-done on the teaching AMI to save time:

ssh -i ~/.ec2/your-key.pem ubuntu@public-ip

sudo apt-get update
sudo apt-get install -y build-essential zlib1g-dev wget curl python-setuptools git \


                        openjdk-7-jdk openjdk-7-jre ruby libncurses5-dev libcurl4-openssl-dev \


                        libbz2-dev unzip pigz bsdmainutils
wget https://raw.githubusercontent.com/bcbio/bcbio-nextgen/master/scripts/bcbio_nextgen_install.py
python bcbio_nextgen_install.py /usr/local/share/bcbio --tooldir /usr/local \


       --genomes hg38 --aligners bwa --sudo --isolate -u development

mkdir -p run
cd run
wget https://raw.githubusercontent.com/bcbio/bcbio-nextgen/master/config/teaching/cancer-syn3-chr6-prep.sh
bash cancer-syn3-chr6-prep.sh

cd cancer-syn3-chr6/work
bcbio_nextgen.py ../config/cancer-syn3-chr6.yaml -n 16

https://github.com/bcbio/bcbio-nextgen

SMALL RNA-SEQ

Data was analyzed with bcbio-nextgen (https://github.com/bcbio/bcbio-nextgen) using piDNA to detect the adapter, cutadapt to remove it, STAR/bowtie to align against the genome and seqcluster to detect small RNA transcripts. miRNAs were detected using miraligner tool with miRBase as the reference miRNA database. tRNA profiles were detected using tdrmapper tool. mirdeep2 was used for discovery of novel miRNAs. FastQC was used for QC metrics and multiqc for reporting.

Download BIB format: https://github.com/bcbio/bcbio-nextgen/tree/master/docs/contents/misc/bcbio-smallrna.bib

AUTHOR

Brad Chapman

COPYRIGHT

2019, bcbio-nextgen contributors