Skip to content

Commit 6f88e2d

Browse files
author
Chenghao Zhu
authored
Merge pull request #28 from uclahs-cds/cz-update-output-dir
Enabled input and output directly from and to blcds data storage
2 parents b036799 + e6db3e1 commit 6f88e2d

7 files changed

Lines changed: 638 additions & 114 deletions

File tree

‎.gitattributes‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
*.config linguist-language=groovy

‎README.md‎

Lines changed: 168 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,175 @@
1-
# pipeline-align-DNA
1+
# Call DNA Align Nextflow Pipeline for BWA Alignment of Paired-End Reads
22

3-
For pipeline documentation, please refer [here](https://uclahs.box.com/s/kl4pacq332bprpe9lnfams0l8vglmg30)
3+
1. [Overview](#Overview)
4+
2. [How To Run](#How-To-Run)
5+
3. [Flow Diagram](#Flow-Diagram)
6+
4. [Pipeline Steps](#Pipeline-Steps)
7+
5. [Inputs](#Inputs)
8+
5. [Outputs](#Outputs)
9+
6. [Testing and Validation](#Testing-and-Validation)
10+
7. [References](#References)
411

5-
For docker images, please view our Dockerhub repository [here](https://hub.docker.com/orgs/blcdsdockerregistry/repositories)
12+
## Overview
613

7-
## How to run the pipeline
14+
The align-DNA nextflow pipeline, aligns paired-end data utilizing [BWA-MEM2](https://github.com/bwa-mem2/bwa-mem2), [Picard](https://github.com/broadinstitute/picard) Tools and [Samtools](https://github.com/samtools/samtools). The pipeline has been engineered to run in a 4 layer stack in a cloud-based scalable environment of CycleCloud, Slurm, Nextflow and Docker. Additionally, it has been validated with the SMC-HET dataset and reference GRCh38, where paired-end fastq’s were created with BAM Surgeon.
815

9-
**The pipeline should be run WITH A SINGLE SAMPLE AT TIME IN SINGLE NODE MODE. Otherwise resource allocation and Nextflow errors could cause the pipeline to fail**
16+
The pipeline should be run **WITH A SINGLE SAMPLE AT A TIME**. Otherwise resource allocation and Nextflow errors could cause the pipeline to fail.
1017

11-
1. Update the params section of the .config file
18+
<b><i>Developer's Notes:</i></b>
1219

13-
2. Update the input csv
20+
> For some reads with low mapping qualities, BWA-MEM2 assigns them to different genomic positions when using different CPU-numbers. If you want to 100% reproduce a run, the same CPU-number (`bwa_mem_number_of_cpus`) needs to be set.
1421
15-
3. See the submission script, [here](https://github.com/uclahs-cds/tool-submit-nf), to submit your pipeline
22+
> BWA-MEM2 now only supports five CPU instruction set, AVX, AVX2, AVX512, SSE4.1 and SSE4.2. However we only tested the pipeline on AVX2 and AVX512 CPUs.
23+
24+
> We performed a benchmarking on our SLURM cluster. Using 56 CPUs for alignment (`bwa_mem_number_of_cpus`) gives it the best performance. See [Testing and Validation](#Testing-and-Validation).
25+
26+
---
27+
28+
## How To Run
29+
30+
Pipelines should be run **WITH A SINGLE SAMPLE AT TIME**. Otherwise resource allocation and Nextflow errors could cause the pipeline to fail
31+
32+
1. Make sure the pipeline is already downloaded to yoru machine. You can either download the stable release or the dev version by cloning the repo.
33+
34+
2. Create a config file for input, output, and parameters. An example can be found [here](pipeline/config/align-DNA.config). See [Inputs](#Inputs) for description of each variables in the config file.
35+
36+
3. Update the input csv. See [Inputs](#Inputs) for the columns needed. All columns must exist in order to run the pipeline. An example can be found [here](pipeline/inputs/align-DNA.inputs.csv). The example csv is a single-lane sample, however this pipeline can take multi-lane sample as well, with each record in the csv file representing a lane (a paire of fastq). All records must have the same value in the **sample** column.
37+
38+
4. See the submission script, [here](https://github.com/uclahs-cds/tool-submit-nf), to submit your pipeline
39+
40+
---
41+
42+
## Flow Diagram
43+
44+
A directed acyclic graph of your pipeline.
45+
46+
![align-DNA flow diagram](flowchart-diagram.drawio.svg?raw=true)
47+
48+
---
49+
50+
## Pipeline Steps
51+
52+
### 1A. Alignment
53+
54+
The first step of the pipeline utilizes [BWA-MEM2](https://github.com/bwa-mem2/bwa-mem2) to align paired reads, (see Tools and Infrastructure Section for details). BWA-MEM2 is the successor for the well-known aligner BWA. The bwa-mem2 mem command utilizes the -M option for marking shorter splitsas secondary. This allows for compatibility with Picard Tools in downstream process and in particular prevents the underlying library of Picard Tools from recognizing these splits as duplicate reads (read names). Additionally, the -t option is utilized toincrease the number of threads used for alignment. The number of threads used in this step is by default to allow at least 2.5Gb memory per CPU, because of the large memory usage by BWA-MEM2. This can be overwritten by setting the bwa_mem_number_of_cpus parameter from the config file. For more details regarding the specific command that was run please refer to the Source Code section.
55+
56+
### 1B. Convert Align SAM File to BAM Format
57+
58+
In the sampe step of the pipeline utilizes Samtool’sview command to convert the aligned SAM files into the compressed BAM format, (see Tools and Infrastructure Section for details). The Samtools view command utilizes the -s option for increasing the speed by removing duplicates and outputs the reads as they are ordered in the file. Additionally, the -b option ensures the output is in BAM format and the -@ option is utilized to increase the number of threads. For more details regarding the specific command that was run please refer to the Source Code section.
59+
60+
### 2. Sort BAM Files in Coordinate Order
61+
62+
The next step of the pipeline utilizes Picard Tool’s SortSam command to sort the aligned BAM files in coordinate order that is needed for downstream tools, (see Tools and Infrastructure Section for details). The SortSam command utilizes the VALIDATION_STRINGENCY=LENIENT option to indicate how errors should be handled and keep the process running if possible. Additionally, the SORT_ORDER option is utilized to ensure the file is sorted in coordinate order, opposed to being sorted by read group name. For more details regarding the specific command that was run please refer to the Source Code section.
63+
64+
## 3. Mark Duplicates in BAM Files
65+
66+
The next step of the pipeline utilizes Picard Tool’s MarkDuplicates command to mark duplicates in the BAM files, (see Tools and Infrastructure Section for details). The MarkDuplicates command utilizes the VALIDATION_STRINGENCY=LENIENT option to indicate how errors should be handled and keep the process running if possible. Additionally, the Program_Record_Id is set to “MarkDuplicates”. For more details regarding the specific command that was run please refer to the Source Code section
67+
68+
## 4. Index BAM Files
69+
70+
After marking dup BAM files, the BAM files are then indexed by utilizing Picard Tool’sBuildBamIndex command, (see Tools and Infrastructure Section for details). This utilizes the `VALIDATION_STRINGENCY=LENIENT` option to indicate how errors should be handled and keep the process running if possible. For more details regarding the specific command that was run please refer to the Source Code section.
71+
72+
---
73+
74+
## Inputs
75+
76+
### Input CSV Fields
77+
78+
>The input csv must have all columns below and in the same order. An example of an input csv can be found [here](pipeline/inputs/align-DNA.inputs.csv)
79+
80+
| Field | Type | Description |
81+
|:------|:-----|:------------|
82+
| index | integer | The index of input fastq pairs, starting from 1. |
83+
| read_group_identifier | string | The read group each read blongs to. This is concatenated with the `lane` column (see below) and then passed to the `ID` field of the final BAM. No white space is allowed. For more detail see [here](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups). |
84+
| sequencing_center | string | The sequencing center where the data were produced. This is passed to the `CN` field of the final BAM. No white space is allowed. For more detail see [here](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups) |
85+
| library_identifier | string | The library identifier to be passed to the `LB` field of the final BAM. No white space is allowed. For more detail see [here](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups) |
86+
| platform_technology | string | The platform or technology used to produce the reads. This is passed to the `PL` field of the final BAM. No white space is allowed. For more detail see [here](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups) |
87+
| platform_unit | string | The platform unit to be passed to the `PU` field of the final BAM. No white space is allowed. For more detail see [here](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups) |
88+
| sample | string | The sample name to be passed to the `SM` field of the final BAM. No white space is allowed. For more detail see [here](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups) |
89+
| lane | string | The lane name or index. This is concatenated with the `read_group_identifier` column (see above) and then passed to the `ID` field of the final BAM. No white space is allowed. For more detail see [here](https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups) |
90+
| read1_fastq | path | Absolute path to the R1 fastq file. |
91+
| read2_fastq | path | Absolute path to the R2 fastq file. |
92+
93+
### Config File Parameters
94+
95+
| Input Parameter | Required | Type | Description |
96+
|:----------------|:---------|:-----|:------------|
97+
| `sample_name` | yes | string | The sample name. This is ignored if the output files are directly saved to the Boutros Lab data storage registry, by setting `blcds_registered_dataset_output = true` |
98+
| `input_csv` | yes | path | Absolute path to the input csv. See [here](pipeline/inputs/align-DNA.inputs.csv) for example and above for the detail of required fields. |
99+
| `reference_fasta` | yes | path | Absolute path to the reference genome `fasta` file. The reference genome is used by BWA-MEM2 for alignment. |
100+
| `reference_fasta_index_files` | yes | path | Absolute path to the genome index with a pattern matching. The index must be generated by the `bwa-mem2 index` command using the correct version of BWA-MEM2. |
101+
| `reference_genome_version` | no | string | The genome build version. This is only used when the output files are directly saved to the Boutros Lab data storage registry, by setting `blcds_registered_dataset_output = true`. |
102+
| `output_dir` | yes | path | Absolute path to the directory where the output files to be saved. This is ignored if the output files are directly saved to the Boutros Lab data storage registry, by setting `blcds_registered_dataset_output = true` |
103+
| `temp_dir` | yes | path | Absolute path to the directory where the nextflow's intermediate files are saved. |
104+
| `save_intermediate_files` | yes | boolean | Save intermediate files. If yes, not only the final BAM, but also the unmerged, unsorted, and duplicates unmarked BAM files will also be saved. |
105+
| `cache_intermediate_pipeline_steps` | yes | boolean | Enable cahcing to resume pipeline and the end of the last successful process completion when a pipeline fails (if true the default submission script must be modified). |
106+
| `max_number_of_parallel_jobs` | yes | int | The maximum number of jobs or steps of the pipeline that can be ran in parallel. |
107+
| `bwa_mem_number_of_cpus` | no | int | Number of cores to use for BWA-MEM2. If not set, this will be calculated to ensure at least 2.5Gb memory per core. |
108+
| `blcds_registered_dataset_input` | yes | boolean | Input FASTQs are from the Boutros Lab data registry. |
109+
| `blcds_registered_dataset_output` | yes | boolean | Enable saving final files including BAM and BAM index, and logging files directory to the Boutros Lab Data registry. |
110+
| `blcds_cluster_slurm` | no | boolean | Pipeline is to run on the Slurm cluster. Set to `false` if it is to run on the SGE cluster. This is used only when `blcds_registered_dataset_output = true` and `blcds_registered_dataset_input = false`. It is also ignored if `blcds_mount_dir` is set. |
111+
| `blcds_disease_id` | no | string | The registered disease ID of this dataset from the Boutros Lab data registry. Ignored if `blcds_registered_data_input = true` or `blcds_registered_output = false` |
112+
| `blcds_dataset_id` | no | string | The registered dataset ID of this dataset from the Boutros Lab data registry. Ignored if `blcds_registered_data_input = true` or `blcds_registered_output = false` |
113+
| `blcds_patient_id` | no | string | The registered patient ID of this sample from the Boutros Lab data registry. Ignored if `blcds_registered_data_input = true` or `blcds_registered_output = false` |
114+
| `blcds_sample_id` | no | string | The registered sample ID from the Boutros Lab data registry. Ignored if `blcds_registered_data_input = true` or `blcds_registered_output = false` |
115+
| `blcds_mount_dir` | no | string | The directoyr that the storage is mounted to (e.g., /hot, /data). |
116+
117+
---
118+
119+
## Outputs
120+
121+
| Output | Description |
122+
|:-------|:------------|
123+
| `.bam` | Aligned, sorted, filtered and if needed, merged, BAM file(s) |
124+
| `.bam.bai` | Index file for each BAM file |
125+
| `.bam` files and metrics files | Intermediate outputs for each scientific tool (OPTIONAL) |
126+
| `report.html`, `timeline.html` and `trace.txt` | A Nextflowreport, timeline and trace files |
127+
| `log.command.*` | Process specific logging files created by nextflow. |
128+
129+
---
130+
131+
## Testing and Validation
132+
133+
### Test Data Set
134+
135+
This pipeline was tested using the synthesized SMC-HET dataset as well as a multi-lane real sample CPCG0196-B1, using reference genome version GRCh38. Some benchmarking has been done comparing BWA-MEM2 v2.1, v2.0, and the original BWA. BWA-MEM2 is able to reduce approximately half of the runtime comparing to the original BWA, with the output BAM almost identical. See [here](https://uclahs.app.box.com/file/737244561716) for the benchmarking.
136+
137+
### Validation <version number\>
138+
139+
| metric | Result |
140+
|:-------|:-------|
141+
| raw total sequences | 1.0000000 |
142+
| filtered sequences | NaN |
143+
| sequences | 1.0000000 |
144+
| is sorted | 1.0000000 |
145+
| 1st fragments | 1.0000000 |
146+
| last fragments | 1.0000000 |
147+
| reads mapped | 1.0000000 |
148+
| reads mapped and paired | 1.0000001 |
149+
| reads unmapped | 0.9999950 |
150+
| reads properly paired | 0.9999999 |
151+
| reads paired | 1.0000000 |
152+
| reads duplicated | 0.9999949 |
153+
| reads MQ0 | 1.0000009 |
154+
| reads QC failed | NaN |
155+
| non-primary alignments | 0.9999757 |
156+
| total length | 1.0000000 |
157+
| bases mapped | 1.0000000 |
158+
| bases mapped (cigar) | 1.0000000 |
159+
| bases trimmed | NaN |
160+
| bases duplicated | 0.9999958 |
161+
| mismatches | 0.9999987 |
162+
| error rate | 0.9999987 |
163+
| average length | 1.0000000 |
164+
| maximum length | 1.0000000 |
165+
| average quality | 1.0000000 |
166+
| insert size average | 1.0000000 |
167+
| insert size standard deviation | 1.0000000 |
168+
| inward oriented pairs | 0.9999991 |
169+
| outward oriented pairs | 1.0000477 |
170+
| pairs with other orientation | 0.9999726 |
171+
| pairs on different chromosomes | 1.0000416 |
172+
173+
---
174+
175+
## References

0 commit comments

Comments
 (0)