NAME¶

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¶

AUTOMATED¶

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¶

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

Tune the upgrade with these options:

CUSTOMIZING DATA INSTALLATION¶

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¶

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¶

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¶

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 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¶

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-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¶

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¶

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¶

SINGLE-CELL RNA-SEQ¶

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¶

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¶

STANDARD¶

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¶

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¶

EXAMPLE PIPELINES¶

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¶

$ 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.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¶

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¶

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¶

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¶

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¶

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¶

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¶

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¶

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¶

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 -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¶

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¶

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

IPYTHON PARALLEL¶

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¶

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¶

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¶

execd_params                 S_DESCRIPTORS=20000

CURRENT STATUS¶

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_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¶

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¶

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)¶

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)¶

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)¶

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)¶

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¶

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)¶

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)¶

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¶

VARIANT CALLING¶

RNA-SEQ¶

SMALL RNA-SEQ¶

DOWNSTREAM ANALYSIS¶

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¶

OVERVIEW¶

DEVELOPMENT INFRASTRUCTURE¶

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¶

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¶

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¶

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¶

ABSTRACT¶

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¶

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¶

``

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¶

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¶

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¶

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

AUTHOR¶

COPYRIGHT¶