diff --git a/README.md b/README.md
index b9837d5..f75fa28 100644
--- a/README.md
+++ b/README.md
@@ -222,7 +222,7 @@ Please follow the instructions at the [Singularity website](https://docs.sylabs.
You can first check the available options and parameters by running:
```bash
-nextflow run Juke34/AliNe -r v1.0.1 --help
+nextflow run Juke34/AliNe -r v1.1.0 --help
```
### Profile
@@ -234,7 +234,7 @@ To run the workflow you must select a profile according to the container platfor
The command will look like that:
```bash
-nextflow run Juke34/AliNe -r v1.0.1 -profile docker <rest of paramaters>
+nextflow run Juke34/AliNe -r v1.1.0 -profile docker <rest of paramaters>
```
Another profile is available (/!\\ actually not yet implemented):
@@ -244,7 +244,7 @@ Another profile is available (/!\\ actually not yet implemented):
The use of the `slurm` profile will give a command like this one:
```bash
-nextflow run Juke34/AliNe -r v1.0.1 -profile singularity,slurm <rest of paramaters>
+nextflow run Juke34/AliNe -r v1.1.0 -profile singularity,slurm <rest of paramaters>
```
### Example
@@ -254,13 +254,15 @@ Here, we use the docker container platform, remote read and genome files, specif
```bash
nextflow run Juke34/AliNe \
- -r v1.0.1 \
+ -r v1.1.0 \
-profile docker \
--reads https://github.com/Juke34/AliNe/raw/refs/heads/main/test/illumina/yeast_R1.fastq.gz \
--genome https://raw.githubusercontent.com/Juke34/AliNe/refs/heads/main/test/yeast.fa \
--read_type short_single \
--aligner bbmap,bowtie2,bwaaln,bwamem,bwasw,graphmap2,hisat2,minimap2,ngmlr,nucmer,star,subread,sublong \
--trimming_fastp \
+ --fastqc \
+ --samtools_statss \
--star_options "--genomeSAindexNbases 9"
```
@@ -271,25 +273,25 @@ Test data are included in the AliNe repository in the `test` folder.
Test with short single reads:
```bash
-nextflow run -profile docker,test_illumina_single Juke34/AliNe -r v1.0.1
+nextflow run -profile docker,test_illumina_single Juke34/AliNe -r v1.1.0
```
Test with short paired reads:
```bash
-nextflow run -profile docker,test_illumina_paired Juke34/AliNe -r v1.0.1
+nextflow run -profile docker,test_illumina_paired Juke34/AliNe -r v1.1.0
```
Test with ont reads:
```bash
-nextflow run -profile docker,test_ont Juke34/AliNe -r v1.0.1
+nextflow run -profile docker,test_ont Juke34/AliNe -r v1.1.0
```
Test with pacbio reads:
```bash
-nextflow run -profile docker,test_pacbio Juke34/AliNe -r v1.0.1
+nextflow run -profile docker,test_pacbio Juke34/AliNe -r v1.1.0
```
On success you should get a message looking like this:
@@ -329,6 +331,7 @@ On success you should get a message looking like this:
Extra steps
--trimming_fastp run fastp for trimming (default: false)
--fastqc run fastqc on raw and aligned reads (default: false)
+ --samtools_stats run samtools stats on aligned reads (default: false)
--multiqc_config path to the multiqc config file (default: config/multiqc_conf.yml)
Aligner specific options
@@ -389,6 +392,12 @@ Here the description of typical ouput you will get from AliNe:
│ ├── sample1_seqkit_trim_aligner2_sorted.log # Contains the log of the aligner
│ └── sample1_seqkit_trim_aligner2_sorted.bam # Sorted bam output
│
+ ├─── samtools_stats # Samtools stats folder
+ │ ├── aligner1 # Folder with Samtools stats result for aligner1
+ │ │ └── sample1_seqkit_trim_aligner1_sorted.txt # Samtools stats file for sample1
+ │ └── aligner2 # Folder with Samtools stats result for aligner2
+ │ └── sample1_seqkit_trim_aligner2_sorted.txt # Samtools stats file for sample1
+ |
├── fastqc # FastQC statistics folder
│ ├── raw # Folder with FastQC result for raw data
│ │ └── fastqc_sample1_raw_logs # Folder with FastQC result for raw sample1 data
@@ -414,9 +423,9 @@ Here the description of typical ouput you will get from AliNe:
### Statistics
-In order to compare several aligner output you should activate the `--fastqc` parameter. AliNe will run the FastQC program for each output file, thereafter all FastQC file will be gathered in an html file using MultiQC. The resulting file called `multiqc_report.html` can be found in <outdir>/MultiQC ( <outdir> by default is called `alignment_results` and can be setup using `--outdir` AliNe parameter).
+#### FastQC
-To compare the output of multiple aligners, you should enable the `--fastqc` parameter. AliNe will execute the FastQC program for each output file. Subsequently, all FastQC reports will be gathered into an HTML file using MultiQC. The resulting file, named `multiqc_report.html`, can be found in the `<output_directory>/MultiQC directory`. By default, the output directory is named `alignment_results`, but this can be customized using the `--outdir` parameter in AliNe.
+To compare the output of multiple aligners, you can enable the `--fastqc` parameter. AliNe will execute the FastQC program for each output file. Subsequently, all FastQC reports will be gathered into an HTML file using MultiQC. The resulting file, named `multiqc_report.html`, can be found in the `<output_directory>/MultiQC directory`. By default, the output directory is named `alignment_results`, but this can be customized using the `--outdir` parameter in AliNe.
FastQC collects the following information:
* Sequence Counts
@@ -430,9 +439,27 @@ FastQC collects the following information:
* Overrepresented sequences
* Adapter Content
+Among these metrics, "Sequence Duplication Levels", "Per Sequence GC Content" and "Sequence Count" are reported at the top of the `multiqc_report.html` file in a table called `General Statistics` as "% Dups", "%GC" and "M Seqs" accordingly ([see below](#general-statistics)).
+
+
+#### Samtools stats
+
+To compare the output of multiple aligners, you can enable the `--samtools_stats` parameter. AliNe will execute the Samtools stats program for each alignment file. Subsequently, all Samtools stats output will be gathered into an HTML file using MultiQC. The resulting file, named `multiqc_report.html`, can be found in the `<output_directory>/MultiQC directory`. By default, the output directory is named `alignment_results`, but this can be customized using the `--outdir` parameter in AliNe.
+
+Samtools stats produces comprehensive statistics, see [here](http://www.htslib.org/doc/samtools-stats.html) for details.
+
+Among all the produceds metrics, "Error rate", "Non-primary", "Reads mapped", "% Mapped", "Total seqs"" are reported at the top of the `multiqc_report.html` file in a table called `General Statistics` ([see below](#general-statistics)).
+These fields correspond to the following information :
+
+ * Error rate - The percentage of mismatches between the aligned reads and the reference genome (mismatches (NM) / bases mapped (CIGAR))
+ * Non-primary - Non-primary alignments is the number of reads that are aligned to multiple locations in the reference genome.
+ * Reads mapped - The number of reads that are successfully aligned to the reference genome.
+ * % Mapped - The percentage of total reads that are mapped. Calculated from reads mapped / total sequences
+ * Total seqs - The total number of reads in the BAM
+
#### General Statistics
-"Sequence Duplication Levels", "Per Sequence GC Content" and "Sequence Count" are reported at the top of the `multiqc_report.html` file in a table called `General Statistics` as "% Dups", "%GC" and "M Seqs" accordingly (see below).
+Some information produced via FastQC or Samtools stats are reported at the top of the `multiqc_report.html` file in a table called `General Statistics` (see below):
<img src="img/multiqc.png" />
@@ -443,61 +470,66 @@ In order to facilitate the reading of this `General Statistics` you can export t
install.packages("dplyr")
install.packages("stringr")
install.packages("tidyr")
+install.packages("knitr")
# Load necessary libraries
library(dplyr)
library(stringr)
library(tidyr)
+library(knitr)
# Read the TSV file
file_path <- "general_stats_table.tsv"
df <- read.delim(file_path, check.names = FALSE)
+# clean sample name to remove suffix _*_samtoolsstats
+df$Sample <- df$Sample |> stringr::str_remove_all("_\\d+_samtoolsstats")
+
# sample name as row name
rownames(df) <- df$Sample
# remove Sample column and clean up the column names
-tmp <- cbind(ID = rownames(df), stack(df[-1])) |>
+tableout <- cbind(ID = rownames(df), stack(df[-1])) |>
transform(ind = as.character(ind) |> stringr::str_remove_all("\\.\\d+"))
# remove na values
-tmp <- tmp[!is.na(tmp$values),]
+tableout <- tableout[!is.na(tableout$values),]
# remove . values
-tmp$values <- tmp$values |> stringr::str_remove_all("^\\.$")
+tableout$values <- tableout$values |> stringr::str_remove_all("^\\.$")
# pivot data
-tmp <- tmp |> pivot_wider(id_cols = ID , names_from = ind, values_from = values,
+tableout <- tableout |> pivot_wider(id_cols = ID , names_from = ind, values_from = values,
values_fn = \(x) paste(unique(x), collapse = ""))
-# print only 4 first columns
-tmp[1:4]
+# print with nice output
+knitr::kable(tmp)
```
You will get a table similar to this one:
```
- ID Dups GC Seqs
- <chr> <chr> <chr> <chr>
- 1 yeast_R1 73.01 43 0.01
- 2 yeast_R1_seqkit_trim 69.21661409043114 44 0.007607999999999999
- 3 yeast_R1_seqkit_trim_STAR_sorted 69.54090258811289 44 0.007689
- 4 yeast_R1_seqkit_trim_bbmap_sorted 69.21661409043114 44 0.007607999999999999
- 5 yeast_R1_seqkit_trim_bowtie2_sorted 69.21661409043114 44 0.007607999999999999
- 6 yeast_R1_seqkit_trim_bwaaln_sorted 69.21661409043114 44 0.007607999999999999
- 7 yeast_R1_seqkit_trim_bwamem_sorted 69.18933123111286 44 0.007611
- 8 yeast_R1_seqkit_trim_bwasw_sorted 69.23279033105622 44 0.007612
- 9 yeast_R1_seqkit_trim_graphmap2_sorted 69.21661409043114 44 0.007607999999999999
-10 yeast_R1_seqkit_trim_hisat2_sorted 69.21661409043114 44 0.007607999999999999
-11 yeast_R1_seqkit_trim_kallisto_sorted 69.21661409043114 44 0.007607999999999999
-12 yeast_R1_seqkit_trim_minimap2_sorted 69.63058976020739 44 0.007715
-13 yeast_R1_seqkit_trim_nucmer.fixed_sorted 64.70588235294117 42 0.00016999999999999999
-14 yeast_R1_seqkit_trim_sublong 53.23343848580442 44 0.007607999999999999
-15 yeast_R1_seqkit_trim_subread 69.21661409043114 44 0.007607999999999999
+|ID |Dups |GC |Seqs |Error rate |Non-primary |Reads mapped |% Mapped |Total seqs |
+|:-----------------------------------|:-----------------|:--|:--------------------|:-------------------|:----------------------|:-----------------------|:------------------|:----------|
+|yeast_R1 |73.01 |43 |0.01 | | | | | |
+|yeast_R1_seqkit_STAR_sorted |73.22686241444302 |43 |0.010081 |0 |0.00008099999999999999 |0.000096 |0.96 |0.01 |
+|yeast_R1_seqkit_bbmap_sorted |73.01 |43 |0.01 |2020.4080000000001 |0 |0.000099 |0.9900000000000001 |0.01 |
+|yeast_R1_seqkit_bowtie2_sorted |73.01 |43 |0.01 |11.46977 |0 |0.00008599999999999999 |0.86 |0.01 |
+|yeast_R1_seqkit_bowtie_sorted |73.01 |43 |0.01 |49.57576 |0 |0.000032999999999999996 |0.33 |0.01 |
+|yeast_R1_seqkit_bwaaln_sorted |73.01 |43 |0.01 |0 |0 |0 |0 |0.01 |
+|yeast_R1_seqkit_bwamem_sorted |72.98080767692923 |43 |0.010003999999999999 |0.37602579999999997 |0 |0.001914 |19.139999999999997 |0.01 |
+|yeast_R1_seqkit_bwasw_sorted |73.0288797841511 |43 |0.010007 |0.2639272 |0 |0.0018369999999999999 |18.357149995003496 |0.010007 |
+|yeast_R1_seqkit_graphmap2_sorted |73.01 |43 |0.01 |0 |0 |0 |0 |0.01 |
+|yeast_R1_seqkit_hisat2_sorted |73.01 |43 |0.01 |0 |0 |0.000004 |0.04 |0.01 |
+|yeast_R1_seqkit_kallisto_sorted |73.01 |43 |0.01 |0 |0 |0.00181 |18.099999999999998 |0.01 |
+|yeast_R1_seqkit_minimap2_sorted |73.28848436881678 |43 |0.010107999999999999 |0.3211603 |0.000108 |0.000252 |2.52 |0.01 |
+|yeast_R1_seqkit_nucmer.fixed_sorted |64.73988439306359 |42 |0.000173 |0.2809365 |0 |0.000169 |100 |0.000169 |
+|yeast_R1_seqkit_sublong_sorted |45.26 |43 |0.01 |0 |0 |0.0012339999999999999 |12.34 |0.01 |
+|yeast_R1_seqkit_subread_sorted |73.01 |43 |0.01 |0.1395645 |0 |0.00156 |15.6 |0.01 |
```
## Integrating AliNe in another nf pipeline
-In nextflow it is possible to call external workflow like AliNe from another workflow.
+In Nextflow it is possible to call external workflow like AliNe from another workflow.
This require to write a dedicated process that will call AliNe and get back the results.
A complete guide how to do so is axailable [here](https://github.com/mahesh-panchal/nf-cascade?tab=readme-ov-file)
diff --git a/aline.nf b/aline.nf
index a544330..d6b81bc 100644
--- a/aline.nf
+++ b/aline.nf
@@ -64,6 +64,7 @@ params.sublong_options = '-X'// -X turn on the RNA-seq mode.
// Report params
params.fastqc = false
+params.samtools_stats = false
params.multiqc_config = "$baseDir/config/multiqc_conf.yml"
// other
@@ -337,9 +338,10 @@ General Parameters
outdir : ${params.outdir}
Report Parameters
- MultiQC parameters
- fastqc : ${params.fastqc}
- multiqc_config : ${params.multiqc_config}
+ fastqc : ${params.fastqc}
+ samtools_stats : ${params.samtools_stats}
+ multiqc_config : ${params.multiqc_config}
+
"""
log.info printAlignerOptions(aligner_list, annotation_file, params.star_index_options)
@@ -376,6 +378,11 @@ include {samtools_sort as samtools_sort_bbmap; samtools_sort as samtools_sort_bo
samtools_sort as samtools_sort_hisat2; samtools_sort as samtools_sort_minimap2; samtools_sort as samtools_sort_ngmlr;
samtools_sort as samtools_sort_novoalign; samtools_sort as samtools_sort_nucmer;
samtools_sort as samtools_sort_sublong } from "$baseDir/modules/samtools.nf"
+include {samtools_stats as samtools_stats_ali_bbmap; samtools_stats as samtools_stats_ali_bowtie ; samtools_stats as samtools_stats_ali_bowtie2 ;
+ samtools_stats as samtools_stats_ali_bwaaln; samtools_stats as samtools_stats_ali_bwamem; samtools_stats as samtools_stats_ali_bwasw; samtools_stats as samtools_stats_ali_graphmap2 ;
+ samtools_stats as samtools_stats_ali_hisat2; samtools_stats as samtools_stats_ali_kallisto; samtools_stats as samtools_stats_ali_minimap2; samtools_stats as samtools_stats_ali_ngmlr;
+ samtools_stats as samtools_stats_ali_novoalign ; samtools_stats as samtools_stats_ali_nucmer; samtools_stats as samtools_stats_ali_star; samtools_stats as samtools_stats_ali_subread ;
+ samtools_stats as samtools_stats_ali_sublong } from "$baseDir/modules/samtools.nf"
include {seqtk_sample} from "$baseDir/modules/seqtk.nf"
include {subread_index; subread; sublong_index; sublong} from "$baseDir/modules/subread.nf"
include {prepare_star_index_options; star_index; star; star2pass} from "$baseDir/modules/star.nf"
@@ -609,6 +616,10 @@ workflow align {
fastqc_ali_bbmap(bbmap_ali, "fastqc/bbmap", "bbmap")
logs.concat(fastqc_ali_bbmap.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_bbmap(bbmap_ali, genome.collect(), "samtools_stats/bbmap", "bbmap")
+ logs.concat(samtools_stats_ali_bbmap.out).set{logs} // save log
+ }
}
// ------------------- BOWTIE -----------------
@@ -628,6 +639,10 @@ workflow align {
fastqc_ali_bowtie(bowtie_ali, "fastqc/bowtie", "bowtie")
logs.concat(fastqc_ali_bowtie.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_bowtie(bowtie_ali, genome.collect(), "samtools_stats/bowtie", "bowtie")
+ logs.concat(samtools_stats_ali_bowtie.out).set{logs} // save log
+ }
}
// ------------------- BOWTIE2 -----------------
@@ -649,6 +664,10 @@ workflow align {
fastqc_ali_bowtie2(bowtie2_ali, "fastqc/bowtie2", "bowtie2")
logs.concat(fastqc_ali_bowtie2.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_bowtie2(bowtie2_ali, genome.collect(), "samtools_stats/bowtie2", "bowtie2")
+ logs.concat(samtools_stats_ali_bowtie2.out).set{logs} // save log
+ }
}
// ------------------- BWA -----------------
@@ -671,6 +690,10 @@ workflow align {
fastqc_ali_bwaaln(bwaaln_ali, "fastqc/bwaaln", "bwaaln")
logs.concat(fastqc_ali_bwaaln.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_bwaaln(bwaaln_ali, genome.collect(), "samtools_stats/bwaaln", "bwaaln")
+ logs.concat(samtools_stats_ali_bwaaln.out).set{logs} // save log
+ }
}
if ("bwamem" in aligner_list){
// align
@@ -688,6 +711,10 @@ workflow align {
fastqc_ali_bwamem(bwamem_ali, "fastqc/bwamem", "bwamem")
logs.concat(fastqc_ali_bwamem.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_bwamem(bwamem_ali, genome.collect(), "samtools_stats/bwamem", "bwamem")
+ logs.concat(samtools_stats_ali_bwamem.out).set{logs} // save log
+ }
}
if ("bwasw" in aligner_list){
// align
@@ -705,6 +732,10 @@ workflow align {
fastqc_ali_bwasw(bwasw_ali, "fastqc/bwasw", "bwasw")
logs.concat(fastqc_ali_bwasw.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_bwasw(bwasw_ali, genome.collect(), "samtools_stats/bwasw", "bwasw")
+ logs.concat(samtools_stats_ali_bwasw.out).set{logs} // save log
+ }
}
}
@@ -727,6 +758,10 @@ workflow align {
fastqc_ali_graphmap2(graphmap2_ali, "fastqc/graphmap2", "graphmap2")
logs.concat(fastqc_ali_graphmap2.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_graphmap2(graphmap2_ali, genome.collect(), "samtools_stats/graphmap2", "graphmap2")
+ logs.concat(samtools_stats_ali_graphmap2.out).set{logs} // save log
+ }
}
// ------------------- HISAT2 -----------------
@@ -748,6 +783,10 @@ workflow align {
fastqc_ali_hisat2(hisat2_ali, "fastqc/hisat2", "hisat2")
logs.concat(fastqc_ali_hisat2.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_hisat2(hisat2_ali, genome.collect(), "samtools_stats/hisat2", "hisat2")
+ logs.concat(samtools_stats_ali_hisat2.out).set{logs} // save log
+ }
}
// ------------------- KALLISTO -----------------
@@ -765,6 +804,10 @@ workflow align {
fastqc_ali_kallisto(kallisto_ali, "fastqc/kallisto", "kallisto")
logs.concat(fastqc_ali_kallisto.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_kallisto(kallisto_ali, genome.collect(), "samtools_stats/kallisto", "kallisto")
+ logs.concat(samtools_stats_ali_kallisto.out).set{logs} // save log
+ }
}
// ------------------- minimap2 -----------------
if ("minimap2" in aligner_list ){
@@ -785,6 +828,10 @@ workflow align {
fastqc_ali_minimap2(minimap2_ali, "fastqc/minimap2", "minimap2")
logs.concat(fastqc_ali_minimap2.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_minimap2(minimap2_ali, genome.collect(), "samtools_stats/minimap2", "minimap2")
+ logs.concat(samtools_stats_ali_minimap2.out).set{logs} // save log
+ }
}
// --------------------- NGMLR --------------------
if ("ngmlr" in aligner_list ){
@@ -799,10 +846,14 @@ workflow align {
// save aligned reads
sorted_bam.concat(ngmlr_ali).set{sorted_bam}
// stat on aligned reads
- if(params.fastqc){
+ if (params.fastqc){
fastqc_ali_ngmlr(ngmlr_ali, "fastqc/ngmlr", "ngmlr")
logs.concat(fastqc_ali_ngmlr.out).set{logs} // save log
}
+ if (params.samtools_stats){
+ samtools_stats_ali_ngmlr(ngmlr_ali, genome.collect(), "samtools_stats/ngmlr", "ngmlr")
+ logs.concat(samtools_stats_ali_ngmlr.out).set{logs} // save log
+ }
}
// ------------------- novoalign -----------------
@@ -819,10 +870,14 @@ workflow align {
// save aligned reads
sorted_bam.concat(novoalign_ali).set{sorted_bam}
// stat on aligned reads
- if(params.fastqc){
+ if (params.fastqc){
fastqc_ali_novoalign(novoalign_ali, "fastqc/novoalign", "novoalign")
logs.concat(fastqc_ali_novoalign.out).set{logs} // save log
}
+ if (params.samtools_stats){
+ samtools_stats_ali_novoalign(novoalign_ali, genome.collect(), "samtools_stats/novoalign", "novoalign")
+ logs.concat(samtools_stats_ali_novoalign.out).set{logs} // save log
+ }
}
// ------------------- nucmer (mummer4) -----------------
@@ -838,10 +893,14 @@ workflow align {
// save aligned reads
sorted_bam.concat(nucmer_ali).set{sorted_bam}
// stat on aligned reads
- if(params.fastqc){
+ if (params.fastqc){
fastqc_ali_nucmer(nucmer_ali, "fastqc/nucmer", "nucmer")
logs.concat(fastqc_ali_nucmer.out).set{logs} // save log
}
+ if (params.samtools_stats){
+ samtools_stats_ali_nucmer(nucmer_ali, genome.collect(), "samtools_stats/nucmer", "nucmer")
+ logs.concat(samtools_stats_ali_nucmer.out).set{logs} // save log
+ }
}
// ------------------- STAR -----------------
@@ -869,6 +928,10 @@ workflow align {
fastqc_ali_star(star_ali, "fastqc/star", "star")
logs.concat(fastqc_ali_star.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_star(star_ali, genome.collect(), "samtools_stats/star", "star")
+ logs.concat(samtools_stats_ali_star.out).set{logs} // save log
+ }
}
// ---------------- subread -----------------
@@ -885,6 +948,10 @@ workflow align {
fastqc_ali_subread(subread_ali, "fastqc/subread", "subread")
logs.concat(fastqc_ali_subread.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_subread(subread_ali, genome.collect(), "samtools_stats/subread", "subread")
+ logs.concat(samtools_stats_ali_subread.out).set{logs} // save log
+ }
}
// ---------------- sublong -----------------
@@ -903,6 +970,10 @@ workflow align {
fastqc_ali_sublong(sublong_ali, "fastqc/sublong", "sublong")
logs.concat(fastqc_ali_sublong.out).set{logs} // save log
}
+ if(params.samtools_stats){
+ samtools_stats_ali_sublong(sublong_ali, genome.collect(), "samtools_stats/sublong", "sublong")
+ logs.concat(samtools_stats_ali_sublong.out).set{logs} // save log
+ }
}
// ------------------- MULTIQC -----------------
@@ -978,6 +1049,7 @@ def helpMSG() {
Extra steps
--trimming_fastp run fastp for trimming (default: false)
--fastqc run fastqc on raw and aligned reads (default: false)
+ --samtools_stats run samtools stats on aligned reads (default: false)
--multiqc_config path to the multiqc config file (default: config/multiqc_conf.yml)
Aligner specific options
diff --git a/config/multiqc_conf.yml b/config/multiqc_conf.yml
index ca837fe..b5ff552 100644
--- a/config/multiqc_conf.yml
+++ b/config/multiqc_conf.yml
@@ -3,86 +3,140 @@ title: "AliNe - Nextflow Alignment pipeline (https://github.com/Juke34/AliNe)"
run_modules:
- fastqc
- fastp
- - bowtie2
- - hisat2
- - star
- - kallisto
+ - samtools
module_order:
- fastqc:
name: "FastQC (raw)"
path_filters:
- - "*raw_logs*"
+ - "*raw_logs/*"
- fastp
- fastqc:
name: "FastQC (trimmed)"
path_filters:
- - "*trimmed_logs*"
+ - "*trimmed_logs/*"
# issue #2866 - bbmap
- fastqc:
name: "FastQC (bbmap)"
path_filters:
- - "*bbmap_logs*"
- - bowtie2:
- name: "bowtie2"
+ - "*bbmap_logs/*"
+ - samtools:
+ name: "Samtools stats (bbmap)"
+ path_filters:
+ - "*bbmap_*.txt"
+ - fastqc:
+ name: "FastQC (bowtie)"
+ path_filters:
+ - "*bowtie_logs/*"
+ - samtools:
+ name: "Samtools stats (bowtie)"
+ path_filters:
+ - "*bowtie_*.txt"
- fastqc:
name: "FastQC (bowtie2)"
path_filters:
- - "*bowtie2_logs*"
+ - "*bowtie2_logs/*"
+ - samtools:
+ name: "Samtools stats (bowtie2)"
+ path_filters:
+ - "*bowtie2_*.txt"
- fastqc:
name: "FastQC (bwaaln)"
path_filters:
- - "*bwaaln_logs*"
+ - "*bwaaln_logs/*"
+ - samtools:
+ name: "Samtools stats (bwaaln)"
+ path_filters:
+ - "*bwaaln_*.txt"
- fastqc:
name: "FastQC (bwamem)"
path_filters:
- - "*bwamem_logs*"
+ - "*bwamem_logs/*"
+ - samtools:
+ name: "Samtools stats (bwamem)"
+ path_filters:
+ - "*bwamem_*.txt"
- fastqc:
name: "FastQC (bwasw)"
path_filters:
- - "*bwasw_logs*"
+ - "*bwasw_logs/*"
+ - samtools:
+ name: "Samtools stats (bwasw)"
+ path_filters:
+ - "*bwasw_*.txt"
- fastqc:
name: "FastQC (graphmap2)"
path_filters:
- - "*graphmap2_logs*"
- - fastqc:
- name: "FastQC (hisat)"
+ - "*graphmap2_logs/*"
+ - samtools:
+ name: "Samtools stats (graphmap2)"
path_filters:
- - "*hisat_logs*"
- - hisat2
+ - "*graphmap2_*.txt"
- fastqc:
name: "FastQC (hisat2)"
path_filters:
- - "*hisat2_logs*"
- - kallisto
+ - "*hisat2_logs/*"
+ - samtools:
+ name: "Samtools stats (hisat2)"
+ path_filters:
+ - "*hisat2_*.txt"
- fastqc:
name: "FastQC (kallisto)"
path_filters:
- - "*kallisto_logs*"
+ - "*kallisto_logs/*"
+ - samtools:
+ name: "Samtools stats (kallisto)"
+ path_filters:
+ - "*kallisto_*.txt"
- fastqc:
name: "FastQC (minimap2)"
path_filters:
- - "*minimap2_logs*"
+ - "*minimap2_logs/*"
+ - samtools:
+ name: "Samtools stats (minimap2)"
+ path_filters:
+ - "*minimap2_*.txt"
- fastqc:
name: "FastQC (nucmer)"
path_filters:
- - "*nucmer_logs*"
+ - "*nucmer_logs/*"
+ - samtools:
+ name: "Samtools stats (nucmer)"
+ path_filters:
+ - "*nucmer*.txt"
- fastqc:
- name: "FastQC (subread)"
+ name: "FastQC (STAR)"
+ path_filters:
+ - "*star_logs/*"
+ - "*STAR_logs/*"
+ - samtools:
+ name: "Samtools stats (STAR)"
path_filters:
- - "*subread_logs*"
+ - "*star_*.txt"
+ - "*STAR_*.txt"
- fastqc:
- name: "FastQC (sublong)"
+ name: "FastQC (STARlong)"
+ path_filters:
+ - "*starlong_logs/*"
+ - "*STARlong_logs/*"
+ - samtools:
+ name: "Samtools stats (STARlong)"
path_filters:
- - "*sublong_logs*"
- - star
+ - "*starlong_*.txt"
+ - "*STARlong_*.txt"
- fastqc:
- name: "FastQC (STAR)"
+ name: "FastQC (subread)"
+ path_filters:
+ - "*subread_logs/*"
+ - samtools:
+ name: "Samtools stats (subread)"
path_filters:
- - "*star_logs*"
- - "*STAR_logs*"
+ - "*subread_*.txt"
- fastqc:
- name: "FastQC (STARlong)"
+ name: "FastQC (sublong)"
+ path_filters:
+ - "*sublong_logs/*"
+ - samtools:
+ name: "Samtools stats (sublong)"
path_filters:
- - "*starlong_logs*"
- - "*STARlong_logs*"
\ No newline at end of file
+ - "*sublong_*.txt"
\ No newline at end of file
diff --git a/img/multiqc.png b/img/multiqc.png
index b42cc36..b97f10d 100644
Binary files a/img/multiqc.png and b/img/multiqc.png differ
diff --git a/modules/samtools.nf b/modules/samtools.nf
index a3f49d7..a5f5cd2 100644
--- a/modules/samtools.nf
+++ b/modules/samtools.nf
@@ -63,4 +63,30 @@ process samtools_sort {
"""
samtools sort -@ ${task.cpus} ${bam} -o ${bam.baseName}_sorted.bam
"""
-}
\ No newline at end of file
+}
+
+
+/*
+http://www.htslib.org/doc/samtools-stats.html
+Produces comprehensive statistics from alignment file
+*/
+process samtools_stats {
+ label 'samtools'
+ tag "$sample"
+ publishDir "${params.outdir}/${outpath}", mode: 'copy'
+
+ input:
+ tuple val(sample), path(bam)
+ path(genome_fasta)
+ val outpath
+ val suffix
+
+ output:
+ path ("*.txt"), emit: bam_samtools_stats
+
+ script:
+
+ """
+ samtools stats --reference ${genome_fasta} ${bam} > ${bam.baseName}.txt
+ """
+}
diff --git a/modules/subread.nf b/modules/subread.nf
index 8412b0e..8d21d2e 100644
--- a/modules/subread.nf
+++ b/modules/subread.nf
@@ -53,7 +53,7 @@ process subread {
}
// remove fastq.gz
- def fileName = fastq[0].baseName.replace('.fastq','') + "_subread"
+ def fileName = fastq[0].baseName.replace('.fastq','') + "_subread_sorted"
// prepare index name
def index_prefix = genome.baseName + "_index"
@@ -72,7 +72,7 @@ process subread {
}
}
"""
- subread-align -T ${task.cpus} ${read_orientation} -i ${index_prefix} ${input} -o ${fileName}.bam --sortReadsByCoordinates ${params.subread_options} > ${fileName}_subread.log
+ subread-align -T ${task.cpus} ${read_orientation} -i ${index_prefix} ${input} -o ${fileName}.bam --sortReadsByCoordinates ${params.subread_options} > ${fileName}_subread_sorted.log
"""
}
diff --git a/paper/paper.md b/paper/paper.md
index 8dc924d..0109110 100644
--- a/paper/paper.md
+++ b/paper/paper.md
@@ -20,14 +20,14 @@ bibliography: paper.bib
Alignment of sequencing reads against a reference genome is a fundamental step in many bioinformatics workflows. Aligner performance varies by speed, memory efficiency, and accuracy, with some tailored to specific sequencing technologies and others more versatile, which makes the selection of an appropriate tool context-dependent. To streamline this process, we present AliNe (Alignment in Nextflow), a flexible and efficient read alignment pipeline built on the Nextflow framework [@nextflow]. AliNe contains a broad range of commonly used aligners, and is designed to accommodate any high-throughput sequencing projects.
-AliNe supports short reads (both paired-end and single-end) as well as long reads generated by PacBio and Oxford Nanopore Technologies (ONT). It currently supports 16 widely used alignment tools, including BBMap [@bbmap], Bowtie [@bowtie], Bowtie2 [@bowtie2], BWA [@bwaaln], BWA-MEM [@bwamem], BWA-SW [@bwasw], GraphMap2 [@graphmap2], HISAT2 [@hisat2], Kallisto [@kallisto], Minimap2 [@minimap2], ngmlr [@ngmlr], novoAlign [@novoalign], nucmer [@nucmer], STAR (single or two-pass mode) [@star], subread [@subread] and sublong [@subread]. These aligners are integrated into a single, easy-to-use workflow, providing a unified entry point for any project requiring alignment and giving users the flexibility to choose the best tool for their specific data and objectives. AliNe is designed to minimize user inputs and avoid common parameter mistakes ( e.g. scoring system, strandedness, orientation).
+AliNe supports short reads (both paired-end and single-end) as well as long reads generated by PacBio and Oxford Nanopore Technologies (ONT). It is designed to minimize user inputs and avoid common parameter mistakes ( e.g. scoring system, strandedness, orientation). It currently supports 16 widely used alignment tools, including BBMap [@bbmap], Bowtie [@bowtie], Bowtie2 [@bowtie2], BWA [@bwaaln], BWA-MEM [@bwamem], BWA-SW [@bwasw], GraphMap2 [@graphmap2], HISAT2 [@hisat2], Kallisto [@kallisto], Minimap2 [@minimap2], ngmlr [@ngmlr], novoAlign [@novoalign], nucmer [@nucmer], STAR (single or two-pass mode) [@star], subread [@subread] and sublong [@subread]. These aligners are integrated into a single, easy-to-use workflow, providing a unified entry point for any project requiring alignment and giving users the flexibility to choose the tool that best suits their specific data and objectives. To facilitate aligner comparison, AliNe provides a comprehensive set of metrics, including computational resource usage (CPU, memory, runtime, and disk I/O) reported by Nextflow, as well as alignment statistics generated by FastQC [@fastqc] and Samtools stats [@samtools].
# Statement of Need
The complexity of sequencing data analysis, especially with the advent of next-generation sequencing (NGS), has created a demand for robust, scalable, and reproducible workflows that can handle large datasets while allowing users to flexibly choose their preferred analysis tools. Numerous bioinformatics analyses rely on read alignment steps that can play a major role in the quality of the downstream results. Each alignment tool may behave differently and have its own strengths. Users may need to try out several alignment tools in order to achieve optimal results specific to a project's requirements. AliNe empowers researchers with a single-entry point for alignment, eliminating the need to juggle between disparate tools. By integrating multiple aligners into a unified pipeline, AliNe simplifies the process of alignment algorithm selection and establishes a standardized workflow that ensures both reproducibility and transparency across analyses
The majority of pipelines integrating read alignment steps are specialized for particular sequencing technologies and types of analyses [@tronflow] [@nfcorernaseq] [@sarek] [@readmapping]. There remains a lack of unified, versatile workflows capable of integrating aligners into a single, user-friendly platform.
-AliNe addresses this challenge by providing a comprehensive, streamlined solution for alignment tasks across a wide range of sequencing technologies, whatever the field of application. Moreover, AliNe can be integrated as a sub-workflow into other Nextflow pipeline [@nfcascade] and provide a standardize step that offer users a wider choice of alignment tools.
+AliNe addresses this challenge by providing a comprehensive, streamlined solution for alignment tasks across a wide range of sequencing technologies, whatever the field of application. Moreover, AliNe can be integrated as a sub-workflow into other Nextflow pipeline and provide a standardize step that offer users a wider choice of alignment tools.
# Workflow Overview
diff --git a/profiles/test_illumina_single.config b/profiles/test_illumina_single.config
index 9cb4ca1..029b046 100644
--- a/profiles/test_illumina_single.config
+++ b/profiles/test_illumina_single.config
@@ -10,6 +10,8 @@ params {
params.read_type = "short_single"
aligner = 'bbmap,bowtie,bowtie2,bwaaln,bwamem,bwasw,graphmap2,hisat2,kallisto,minimap2,ngmlr,nucmer,star,subread,sublong'
trimming_fastp = true
+ fastqc = true
+ samtools_stats = true
star_options = "--genomeSAindexNbases 9" // the default 14 is too large for the genome size=1351857
multiqc_config = "$baseDir/config/multiqc_conf.yml"
}
\ No newline at end of file