FAQs
- What is a snap file?
- How to create a snap file from fastq file?
- How to create snap file for 10X fastq file?
- Can I run SnapATAC with CellRanger output?
- How to create snap file from bam or bed file?
- How to analyze multiple samples together?
- How to group reads from any subset of cells?
- How to choose the bin size?
snap (Single-Nucleus Accessibility Profiles) file is a hierarchically structured hdf5 file that is specially designed for storing single nucleus ATAC-seq datasets. A snap file (version 4) contains the following sessions: header (HD), cell-by-bin accessibility matrix (AM), cell-by-peak matrix (PM), cell-by-gene matrix (GM), barcode (BD) and fragment (FM).
- HD session contains snap-file version, created date, alignment and reference genome information.
- BD session contains all unique barcodes and corresponding meta data.
- AM session contains cell-by-bin matrices of different resolutions (or bin sizes).
- PM session contains cell-by-peak count matrix. PM session contains cell-by-gene count matrix.
- FM session contains all usable fragments for each cell. Fragments are indexed for fast search. Detailed information about snap file can be found here.
Step 1. Barcode demultiplexing.
We first de-multicomplex FASTQ file by adding the barcode to the beginning of each read in the following format: "@" + "barcode" + ":" + "read_name". Below is one example of demultiplexed fastq file. Because barcode design can be very different between experiments and platforms, we decide not to include this part in the current analysis pipline. However, this can be easily implemented by awk or python script.
$ wget http://renlab.sdsc.edu/r3fang/share/Fang_2019/MOs_snATAC/fastq/CEMBA180306_2B.demultiplexed.R1.fastq.gz
$ wget http://renlab.sdsc.edu/r3fang/share/Fang_2019/MOs_snATAC/fastq/CEMBA180306_2B.demultiplexed.R2.fastq.gz
$ wget http://renlab.sdsc.edu/r3fang/share/Fang_2019/MOs_snATAC/peaks/all_peak.bed
$ wget http://renlab.sdsc.edu/r3fang/share/Fang_2019/MOs_snATAC/mm10.blacklist.bed.gz
$ wget http://renlab.sdsc.edu/r3fang/share/Fang_2019/MOs_snATAC/genes/gencode.vM16.gene.bed
$ wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/mm10.chrom.sizes
$ wget http://renlab.sdsc.edu/r3fang/share/Fang_2019/MOs_snATAC/genome/mm10.fa
$
$ zcat CEMBA180306_2B.demultiplexed.R1.fastq.gz | head
@AGACGGAGACGAATCTAGGCTGGTTGCCTTAC:7001113:920:HJ55CBCX2:1:1108:1121:1892 1:N:0:0
ATCCTGGCATGAAAGGATTTTTTTTTTAGAAAATGAAATATATTTTAAAG
+
DDDDDIIIIHIIGHHHIIIHIIIIIIHHIIIIIIIIIIIIIIIIIIIIIIStep 2. Index reference genome (snaptools).
Index the reference genome before alignment. Here we show how to index the genome using BWA. User can switch to other aligner by --aligner tag, currently snaptools supports bwa, bowtie2 and minimap2. User also needs to specify the folder that stores the aligner executable binary file. For instance, if bwa is installed under /opt/biotools/bwa/bin/bwa, set --path-to-aligner=/opt/biotools/bwa/bin/.
$ which bwa
/opt/biotools/bwa/bin/bwa
$ snaptools index-genome \
--input-fasta=mm10.fa \
--output-prefix=mm10 \
--aligner=bwa \
--path-to-aligner=/opt/biotools/bwa/bin/ \
--num-threads=5Step 3. Alignment (snaptools).
We next align de-multicomplexed reads to the corresponding reference genome using snaptools with following command. After alignment, reads are automatically sorted by read names (--if-sort). User can set mutiple CPUs to speed up this step by setting (--num-threads).
$ snaptools align-paired-end \
--input-reference=mm10.fa \
--input-fastq1=CEMBA180306_2B.demultiplexed.R1.fastq.gz \
--input-fastq2=CEMBA180306_2B.demultiplexed.R2.fastq.gz \
--output-bam=CEMBA180306_2B.bam \
--aligner=bwa \
--path-to-aligner=/opt/biotools/bwa/bin/ \
--read-fastq-command=zcat \
--min-cov=0 \
--num-threads=5 \
--if-sort=True \
--tmp-folder=./ \
--overwrite=TRUE Step 4. Pre-processing (snaptools).
After alignment, we convert pair-end reads into fragments and for each fragment we check the following attributes: 1) mapping quality score MAPQ; 2) whether two ends are appropriately paired according to the alignment flag; 3) fragment length. We only keep those fragments that are 1) properly paired; 2) whose MAPQ is greater than 30 (--min-mapq); 3) with fragment length less than 1000bp (--max-flen). Because the reads have been sorted based on the barcodes (integrated into read names), fragments belonging to the same barcode are naturally grouped together which allows for removing PCR duplicates separately. After alignment and filtration, we generated a snap-format (Single-Nucleus Accessibility Profiles) file that contains meta data, cell-by-bin count matrices, cell-by-peak count matrix. Detailed information about snap file can be found in here. In the below example, barcodes of fragments less than 100 --min-cov=100 are filtered in order to speed up the process and save memory.
$ fetchChromSizes mm10 > mm10.gs
$ snaptools snap-pre \
--input-file=CEMBA180306_2B.bam \
--output-snap=CEMBA180306_2B.snap \
--genome-name=mm10 \
--genome-size=mm10.chrom.sizes \
--min-mapq=30 \
--min-flen=0 \
--max-flen=1000 \
--keep-chrm=TRUE \
--keep-single=FALSE \
--keep-secondary=FALSE \
--overwrite=True \
--min-cov=100 \
--verbose=TrueStep 5. Cell-by-bin matrix (snaptools).
Using snap file, we next create the cell-by-bin matrix. Snap file allows for storing cell-by-bin matrices of different resolutions. In the below example, three cell-by-bin matrices are created with bin size of 1,000, 5,000 and 10,000. Note that this does not create a new file, all the matrices are stored in CEMBA180306_2B.snap.
$ snaptools snap-add-bmat \
--snap-file=CEMBA180306_2B.snap \
--bin-size-list 1000 5000 10000 \
--verbose=TrueCase 1
First run cellranger-atac mkfastq to generate the following demultiplexed fastq files:
Second, for each library recognize that R1 and R3 are the sequencing reads, and I1 is 16bp i5 (10x Barcode), and R2 is i7 (sample index).
Third, snaptools provide a module dex-fastq to integrate the 10X barcode into the read name.
$ snaptools dex-fastq \
--input-fastq=Library1_1_L001_R1_001.fastq.gz \
--output-fastq=Library1_1_L001_R1_001.dex.fastq.gz \
--index-fastq-list Library1_1_L001_R2_001.fastq.gz
$ snaptools dex-fastq \
--input-fastq=Library1_1_L001_R3_001.fastq.gz \
--output-fastq=Library1_1_L001_R3_001.dex.fastq.gz \
--index-fastq-list Library1_1_L001_R2_001.fastq.gz
$ snaptools dex-fastq \
--input-fastq=Library1_2_L001_R1_001.fastq.gz \
--output-fastq=Library1_2_L001_R1_001.dex.fastq.gz \
--index-fastq-list Library1_2_L001_R2_001.fastq.gz
$ snaptools dex-fastq \
--input-fastq=Library1_2_L001_R3_001.fastq.gz \
--output-fastq=Library1_2_L001_R3_001.dex.fastq.gz \
--index-fastq-list Library1_2_L001_R2_001.fastq.gz
# combine these two library
$ cat Library1_1_L001_R1_001.fastq.gz Library1_2_L001_R1_001.fastq.gz > Library1_L001_R1_001.fastq.gz
$ cat Library1_1_L001_R3_001.fastq.gz Library1_2_L001_R3_001.fastq.gz > Library1_L001_R3_001.fastq.gzFour, run the rest of the pipeline using Library1_L001_R1_001.fastq.dex.gz and Library1_L001_R3_001.fastq.dex.gz.
Case 2
In this example, we have two 10x libraries (each processed through a separate Chromium chip channel) that are multiplexed on a single flow-cell. Note that after running cellranger-atac mkfastq, we run a separate instance of the pipeline on each library as shown in above figure.
Next, for each library integrate the barcode into the read name separately:
$ snaptools dex-fastq \
--input-fastq=Library1_S1_L001_R1_001.fastq.gz \
--output-fastq=Library1_S1_L001_R1_001.dex.fastq.gz \
--index-fastq-list Library1_S1_L001_R2_001.fastq.gz
$ snaptools dex-fastq \
--input-fastq=Library1_S1_L001_R3_001.fastq.gz \
--output-fastq=Library1_S1_L001_R3_001.dex.fastq.gz \
--index-fastq-list Library1_S1_L001_R2_001.fastq.gzFour, run the rest of the pipeline using Library1_S1_L001_R1_001.fastq.dex.gz and Library1_S1_L001_R3_001.fastq.dex.gz.
Yes. There are two entry points
- use the position sorted bam file (recommanded).
$ samtools view atac_v1_adult_brain_fresh_5k_possorted_bam.bam
A00519:218:HJYFLDSXX:2:1216:26458:34976 99 chr1 3000138 60 50M = 3000474 385 TGATGACTGCCTCTATTTCTTTAGGGGAAATGGGACTTTTAGTCCATGAA FFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:50 MC:Z:49M AS:i:50 XS:i:37 CR:Z:GGTTGCGAGCCGCAAA CY:Z:FFFFFFFFFFFFFFFF CB:Z:GGTTGCGAGCCGCAAA-1 BC:Z:AACGGTCA QT:Z:FFFFFFFFGP:i:3000137 MP:i:3000522 MQ:i:40 RG:Z:atac_v1_adult_brain_fresh_5k:MissingLibrary:1:HJYFLDSXX:2The cell barcode is embedded in the tag CB:Z:GGTTGCGAGCCGCAAA-1, you can modify the bam file by add the cell barcode GGTTGCGAGCCGCAAA-1 to the beginning of read
# extract the header file
samtools view atac_v1_adult_brain_fresh_5k_possorted_bam.bam -H > atac_v1_adult_brain_fresh_5k_possorted.header.sam
# create a bam file with the barcode embedded into the read name
cat <( cat atac_v1_adult_brain_fresh_5k_possorted.header.sam ) \
<( samtools view atac_v1_adult_brain_fresh_5k_possorted_bam.bam | awk '{for (i=12; i<=NF; ++i) { if ($i ~ "^CB:Z:"){ td[substr($i,1,2)] = substr($i,6,length($i)-5); } }; printf "%s:%s\n", td["CB"], $0 }' ) \
| samtools view -bS - > atac_v1_adult_brain_fresh_5k_possorted.snap.bam
samtools view atac_v1_adult_brain_fresh_5k_possorted.snap.bam | cut -f 1 | head
GGTTGCGAGCCGCAAA-1:A00519:218:HJYFLDSXX:2:1216:26458:34976
GGTTGCGAGCCGCAAA-1:A00519:218:HJYFLDSXX:2:2256:23194:13823
GGTTGCGAGCCGCAAA-1:A00519:218:HJYFLDSXX:2:2546:5258:31955
CTCAGCTAGTGTCACT-1:A00519:218:HJYFLDSXX:1:2428:8648:18349
CTCAGCTAGTGTCACT-1:A00519:218:HJYFLDSXX:1:2428:8648:18349
GAAGTCTGTAACACTC-1:A00519:218:HJYFLDSXX:3:2546:14968:2331
GAAGTCTGTAACACTC-1:A00519:218:HJYFLDSXX:3:2546:14705:2628
GGTTGCGAGCCGCAAA-1:A00519:218:HJYFLDSXX:2:1216:26458:34976
GGTTGCGAGCCGCAAA-1:A00519:218:HJYFLDSXX:2:2256:23194:13823
GGTTGCGAGCCGCAAA-1:A00519:218:HJYFLDSXX:2:2546:5258:31955
Then sort the bam file by read name:
$ samtools sort -n -@ 10 -m 1G atac_v1_adult_brain_fresh_5k_possorted.snap.bam -o atac_v1_adult_brain_fresh_5k.snap.nsrt.bamThen generate the snap file
$ snaptools snap-pre \
--input-file=atac_v1_adult_brain_fresh_5k.snap.nsrt.bam \
--output-snap=atac_v1_adult_brain_fresh_5k.snap \
--genome-name=mm10 \
--genome-size=mm10.chrom.sizes \
--min-mapq=30 \
--min-flen=50 \
--max-flen=1000 \
--keep-chrm=TRUE \
--keep-single=FALSE \
--keep-secondary=False \
--overwrite=True \
--max-num=20000 \
--min-cov=500 \
--verbose=Trueremove temporary files
rm atac_v1_adult_brain_fresh_5k.snap.snap.bam
rm atac_v1_adult_brain_fresh_5k_possorted.header.sam
- use the fragment tsv file. Fragment file is already filtered, this will disable snaptools to generate quality control metrics.
# decompress the gz file
> gunzip atac_v1_adult_brain_fresh_5k_fragments.tsv.gz
# sort the tsv file using the 4th column (barcode column)
$ sort -k4,4 atac_v1_adult_brain_fresh_5k_fragments.tsv > atac_v1_adult_brain_fresh_5k_fragments.bed
# compress the bed file
$ gzip atac_v1_adult_brain_fresh_5k_fragments.bed
# run snap files using the bed file
$ snaptools snap-pre \
--input-file=atac_v1_adult_brain_fresh_5k_fragments.bed.gz \
--output-snap=atac_v1_adult_brain_fresh_5k.snap \
--genome-name=mm10 \
--genome-size=mm10.chrom.sizes \
--min-mapq=30 \
--min-flen=50 \
--max-flen=1000 \
--keep-chrm=TRUE \
--keep-single=FALSE \
--keep-secondary=False \
--overwrite=True \
--max-num=20000 \
--min-cov=500 \
--verbose=TrueIn many cases, you probably already have indexed the reference genome and finished the alignments. In this case, you can skip the first two steps and directly jump to snaptools pre which takes name-sorted bam or bed file as input. Highly recommend to use unfiltered alignment file as input.
For bam file, the cell barcodes need to be integrated into the read name as below:
$ samtools view demo.bam
AAACTACCAGAAAGACGCAGTT:7001113:968:HMYT2BCX2:1:1215:20520:88475 77 * 0 0 * * 0 0CTATGAGCACCGTCTCCGCCTCAGATGTGTATAAGAGACAGCAGAGTAAC @DDBAI??E?1/<DCGECEHEHHGG1@GEHIIIHGGDGE@HIHEEIIHH1 AS:i:0 XS:i:0
AAACTACCAGAAAGACGCAGTT:7001113:968:HMYT2BCX2:1:1215:20520:88475 141 * 0 0 * * 0 0GGCTTGTACAGAGCAAGTGCTGAAGTCCCTTTCTGATGACGTTCAACAGC 0<000/<<1<D1CC111<<1<1<111<111<<CDCF1<1<DHH<<<<C11 AS:i:0 XS:i:0
AAACTACCAGAAAGACGCAGTT:7001113:968:HMYT2BCX2:1:2201:20009:41468 77 * 0 0 * * 0 0CGGTGCCCCTGTCCTGTTCGTGCCCACCGTCTCCGCCTCAGATGTGTATA DDD@D/D<DHIHEHCCF1<<CCCGH?GHI1C1DHIII0<1D1<111<1<1 AS:i:0 XS:i:0
AAACTACCAGAAAGACGCAGTT:7001113:968:HMYT2BCX2:1:2201:20009:41468 141 * 0 0 * * 0 0GAGCGAGGGCGGCAGAGGCAGGGGGAGGAGACCCGGTGGCCCGGCAGGCT 0<00</<//<//<//111000/<</</<0<1<1<//<<0<DCC/<///<D AS:i:0 XS:i:0
$ snaptools snap-pre \
--input-file=demo.bam \
--output-snap=demo.snap \
--genome-name=mm10 \
--genome-size=mm10.chrom.sizes \
--min-mapq=30 \
--min-flen=0 \
--max-flen=1000 \
--keep-chrm=TRUE \
--keep-single=TRUE \
--keep-secondary=False \
--overwrite=True \
--min-cov=100 \
--verbose=Truefor bed file, the barcode should be the 4th column:
$ zcat demo.bed.gz | head
chr2 74358918 74358981 AACGAGAGCTAAAGACGCAGTT
chr6 134212048 134212100 AACGAGAGCTAAAGACGCAGTT
chr10 93276785 93276892 AACGAGAGCTAAAGACGCAGTT
chr2 128601366 128601634 AACGAGAGCTAAAGCGCATGCT
chr16 62129428 62129661 AACGAGAGCTAACAACCTTCTG
chr8 84946184 84946369 AACGAGAGCTAACAACCTTCTG
$ snaptools snap-pre \
--input-file=demo.bed.gz \
--output-snap=demo.snap \
--genome-name=mm10 \
--genome-size=mm10.chrom.sizes \
--min-mapq=30 \
--min-flen=0 \
--max-flen=1000 \
--keep-chrm=TRUE \
--keep-single=TRUE \
--keep-secondary=False \
--overwrite=True \
--min-cov=100 \
--verbose=TrueGroup reads from one cell ATACAGCCTCGC in snap file sample1.snap.
library(SnapATAC);
snap_files = "sample1.snap";
barcode_sel = "ATACAGCCTCGC";
reads.gr = extractReads(barcode_sel, snap_files);
Group reads from multiple barcodes in one snap file.
library(SnapATAC);
barcode_sel = c("ATACAGCCTCGC", "ATACAGCCTCGG")
snap_files = rep("sample1.snap", 2);
reads.gr = extractReads(barcode_sel, snap_files);
Group reads from multiple barcodes and multiple snap files.
library(SnapATAC);
barcode_sel = rep("ATACAGCCTCGC", 2);
snap_files = c("sample1.snap", "sample2.snap");
reads.gr = extractReads(barcode_sel, snap_files);
Because SnapATAC cluster cells using cell-by-bin matrix which makes it easy to combine multiple samples and perform comparative analysis. But it does require cell-by-bin matrix of the same bin size is created for all samples. Here we are analyzing two samples (PBMC_5K and PBMC_10K) together as an example.
$ R
> library(SnapATAC);
> system("wget http://renlab.sdsc.edu/r3fang/share/Fang_2019/published_scATAC/atac_v1_pbmc_10k_fastqs/atac_v1_pbmc_10k.snap");
> system("wget http://renlab.sdsc.edu/r3fang/share/Fang_2019/published_scATAC/atac_v1_pbmc_5k_fastqs/atac_v1_pbmc_5k.snap");
> file.list = c("atac_v1_pbmc_5k.snap", "atac_v1_pbmc_10k.snap");
> sample.list = c("pbmc.5k", "pbmc.10k");
> x.sp = createSnap(file=file.list, sample=sample.list);createSnap will creates a snap object that contains both snap file name and corresponding barcodes from each file and automatically combines these two samples. Note that there is no need to change the barcode
The rest of the analysis is the same as analyzing one sample. See the full example here
SnapATAC clusters cells based on the cell-by-bin matrix, so how to choose the optimal bin size. Unfortunately, there is no absolute answer to this question.
On one hand, we find that the bin size within the range 5kb-50kb does not significantly change the clustering result (as shown in the below figure) while we did notice that a large bin size usually results in relatively lower number of clusters. The differences in clustering perhaps can be compensated by using Louvain clustering with smaller resolution.
The benefit for using large bin size is save of memory. This is especially useful for large dataset. Here is a subjective suggestion for choosing the bin size.
| # of cells | bin size |
|---|---|
| 0-50k | 5kb |
| 50k-100k | 10kb |
| 100k-1M | 20kb-50kb |