diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0a750d89..91b10411 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Enhancements & fixes
 
+- [[#327](https://github.com/nf-core/atacseq/issues/327)]Consistently support `.csi` indices as alternative to `.bai` to allow SAMTOOLS_INDEX to be used with the `-c` flag
+
 ## [[2.1.2](https://github.com/nf-core/atacseq/releases/tag/2.1.2)] - 2022-08-07
 
 ### Enhancements & fixes
diff --git a/subworkflows/local/bam_filter_bamtools.nf b/subworkflows/local/bam_filter_bamtools.nf
index 7c0bba66..dbf38045 100644
--- a/subworkflows/local/bam_filter_bamtools.nf
+++ b/subworkflows/local/bam_filter_bamtools.nf
@@ -32,7 +32,7 @@ workflow BAM_FILTER_BAMTOOLS {
     BAMTOOLS_FILTER
         .out
         .bam
-        .branch { 
+        .branch {
             meta, bam ->
                 single_end: meta.single_end
                     return [ meta, bam ]
@@ -49,11 +49,23 @@ workflow BAM_FILTER_BAMTOOLS {
     }
     ch_versions = ch_versions.mix(SAMTOOLS_INDEX.out.versions.first())
 
+    SAMTOOLS_INDEX.out.bai
+        .join(SAMTOOLS_INDEX.out.csi, by: [0], remainder: true)
+        .map {
+            meta, bai, csi ->
+                if (bai) {
+                    [ meta, bai ]
+                } else {
+                    [ meta, csi ]
+                }
+        }
+        .set { ch_index }
+
     //
     // Run samtools stats, flagstat and idxstats on SE BAM
     //
     BAM_STATS_SAMTOOLS (
-        ch_bam.single_end.join(SAMTOOLS_INDEX.out.bai),
+        ch_bam.single_end.join(ch_index),
         ch_fasta
     )
     ch_versions = ch_versions.mix(BAM_STATS_SAMTOOLS.out.versions.first())
@@ -87,6 +99,7 @@ workflow BAM_FILTER_BAMTOOLS {
     name_bam = SAMTOOLS_SORT.out.bam                                                     // channel: [ val(meta), [ bam ] ]
     bam      = BAM_SORT_STATS_SAMTOOLS.out.bam.mix(ch_bam.single_end)                    // channel: [ val(meta), [ bam ] ]
     bai      = BAM_SORT_STATS_SAMTOOLS.out.bai.mix(SAMTOOLS_INDEX.out.bai)               // channel: [ val(meta), [ bai ] ]
+    csi      = BAM_SORT_STATS_SAMTOOLS.out.csi.mix(SAMTOOLS_INDEX.out.csi)               // channel: [ val(meta), [ csi ] ]
     stats    = BAM_SORT_STATS_SAMTOOLS.out.stats.mix(BAM_STATS_SAMTOOLS.out.stats)       // channel: [ val(meta), [ stats ] ]
     flagstat = BAM_SORT_STATS_SAMTOOLS.out.flagstat.mix(BAM_STATS_SAMTOOLS.out.flagstat) // channel: [ val(meta), [ flagstat ] ]
     idxstats = BAM_SORT_STATS_SAMTOOLS.out.idxstats.mix(BAM_STATS_SAMTOOLS.out.idxstats) // channel: [ val(meta), [ idxstats ] ]
diff --git a/workflows/atacseq.nf b/workflows/atacseq.nf
index c76ac139..5daf5839 100644
--- a/workflows/atacseq.nf
+++ b/workflows/atacseq.nf
@@ -161,7 +161,6 @@ workflow ATACSEQ {
     // SUBWORKFLOW: Alignment with BWA & BAM QC
     //
     ch_genome_bam        = Channel.empty()
-    ch_genome_bam_index  = Channel.empty()
     ch_samtools_stats    = Channel.empty()
     ch_samtools_flagstat = Channel.empty()
     ch_samtools_idxstats = Channel.empty()
@@ -176,7 +175,6 @@ workflow ATACSEQ {
                 }
         )
         ch_genome_bam        = FASTQ_ALIGN_BWA.out.bam
-        ch_genome_bam_index  = FASTQ_ALIGN_BWA.out.bai
         ch_samtools_stats    = FASTQ_ALIGN_BWA.out.stats
         ch_samtools_flagstat = FASTQ_ALIGN_BWA.out.flagstat
         ch_samtools_idxstats = FASTQ_ALIGN_BWA.out.idxstats
@@ -198,7 +196,6 @@ workflow ATACSEQ {
                 }
         )
         ch_genome_bam        = FASTQ_ALIGN_BOWTIE2.out.bam
-        ch_genome_bam_index  = FASTQ_ALIGN_BOWTIE2.out.bai
         ch_samtools_stats    = FASTQ_ALIGN_BOWTIE2.out.stats
         ch_samtools_flagstat = FASTQ_ALIGN_BOWTIE2.out.flagstat
         ch_samtools_idxstats = FASTQ_ALIGN_BOWTIE2.out.idxstats
@@ -222,7 +219,6 @@ workflow ATACSEQ {
             []
         )
         ch_genome_bam        = FASTQ_ALIGN_CHROMAP.out.bam
-        ch_genome_bam_index  = FASTQ_ALIGN_CHROMAP.out.bai
         ch_samtools_stats    = FASTQ_ALIGN_CHROMAP.out.stats
         ch_samtools_flagstat = FASTQ_ALIGN_CHROMAP.out.flagstat
         ch_samtools_idxstats = FASTQ_ALIGN_CHROMAP.out.idxstats
@@ -244,7 +240,6 @@ workflow ATACSEQ {
             params.seq_center ?: ''
         )
         ch_genome_bam        = ALIGN_STAR.out.bam
-        ch_genome_bam_index  = ALIGN_STAR.out.bai
         ch_samtools_stats    = ALIGN_STAR.out.stats
         ch_samtools_flagstat = ALIGN_STAR.out.flagstat
         ch_samtools_idxstats = ALIGN_STAR.out.idxstats
@@ -298,7 +293,17 @@ workflow ATACSEQ {
     // SUBWORKFLOW: Filter BAM file
     //
     MERGED_LIBRARY_FILTER_BAM (
-        MERGED_LIBRARY_MARKDUPLICATES_PICARD.out.bam.join(MERGED_LIBRARY_MARKDUPLICATES_PICARD.out.bai, by: [0]),
+        MERGED_LIBRARY_MARKDUPLICATES_PICARD.out.bam
+            .join(MERGED_LIBRARY_MARKDUPLICATES_PICARD.out.bai, by: [0], remainder: true)
+            .join(MERGED_LIBRARY_MARKDUPLICATES_PICARD.out.csi, by: [0], remainder: true)
+            .map {
+                meta, bam, bai, csi ->
+                    if (bai) {
+                        [ meta, bam, bai ]
+                    } else {
+                        [ meta, bam, csi ]
+                    }
+            },
         PREPARE_GENOME.out.filtered_bed.first(),
         PREPARE_GENOME
             .out
@@ -379,7 +384,16 @@ workflow ATACSEQ {
     MERGED_LIBRARY_FILTER_BAM
         .out
         .bam
-        .join(MERGED_LIBRARY_FILTER_BAM.out.bai, by: [0])
+        .join(MERGED_LIBRARY_FILTER_BAM.out.bai, by: [0], remainder: true)
+        .join(MERGED_LIBRARY_FILTER_BAM.out.csi, by: [0], remainder: true)
+        .map {
+            meta, bam, bai, csi ->
+                if (bai) {
+                    [ meta, bam, bai ]
+                } else {
+                    [ meta, bam, csi ]
+                }
+        }
         .set { ch_bam_bai }
 
     if (params.with_control) {
@@ -477,7 +491,16 @@ workflow ATACSEQ {
     MERGED_LIBRARY_MARKDUPLICATES_PICARD
         .out
         .bam
-        .join(MERGED_LIBRARY_MARKDUPLICATES_PICARD.out.bai, by: [0])
+        .join(MERGED_LIBRARY_MARKDUPLICATES_PICARD.out.bai, by: [0], remainder: true)
+        .join(MERGED_LIBRARY_MARKDUPLICATES_PICARD.out.csi, by: [0], remainder: true)
+        .map {
+            meta, bam, bai, csi ->
+                if (bai) {
+                    [ meta, bam, bai ]
+                } else {
+                    [ meta, bam, csi ]
+                }
+        }
         .join(MERGED_LIBRARY_CALL_ANNOTATE_PEAKS.out.peaks, by: [0])
         .set { ch_bam_peaks }