超能餅乾|SnapATAC分析單細胞ATAC-seq數據(一)

SnapATAC簡介

SnapATAC (Single Nucleus Analysis Pipeline for ATAC-seq) 是一個能夠快速、準確和全面分析單細胞ATAC-seq數據的R包，它可以對單細胞ATAC-seq數據進行常規的數據降維、聚類和批次校正分析，鑑定遠端調控元件並預測其調控的靶基因，調用chromVAR軟體進行motif分析，同時還可以將scRNA-seq和scATAC-seq數據進行整合分析等。

SnapATAC軟體安裝RequirementsPre-print

Rongxin Fang, Sebastian Preissl, Xiaomeng Hou, Jacinta Lucero, Xinxin Wang, Amir Motamedi, Andrew K. Shiau, Eran A. Mukamel, Yanxiao Zhang, M. Margarita Behrens, Joseph Ecker, Bing Ren. Fast and Accurate Clustering of Single Cell Epigenomes Reveals Cis-Regulatory Elements in Rare Cell Types. bioRxiv 615179; doi: https://doi.org/10.1101/615179

Installation

SnapATAC軟體主要由以下兩個組件組成：Snaptools和SnapATAC.

# Install snaptools from PyPI. 
# NOTE: Please use python 2.7 if possible.
# 使用pip安裝snaptools模塊
pip install snaptools

# Install SnapATAC R pakcage (development version).
# 安裝一些依賴R包
install.packages(c('raster','bigmemory','doSNOW','plot3D'))

# 使用devtools安裝SnapATAC包
library(devtools)
install_github("r3fang/SnapATAC")

SnapATAC常見問題匯總1）What is a snap file?

snap (Single-Nucleus Accessibility Profiles)格式文件是一個層級結構的hdf5文件，它可以用來存儲single nucleus ATAC-seq數據集。該文件（version 4）主要由以下幾個部分組成：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: 包含snap文件的版本、創建日期、比對信息和參考基因組信息。

BD session: 包含所有unique細胞的barcodes和相應的meta data信息。

AM session: 包含不同解析度（bin size）下的cell-by-bin數據矩陣。

PM session: 包含cell-by-peak的計數矩陣和cell-by-gene的計數矩陣。

FM session: 包含每個細胞中可用的所有frangments片段信息。

2）How to create a snap file from fastq file?Step 1. Barcode demultiplexing.

首先，我們將barcode信息以"@" + "barcode" + ":" + "read_name"的格式添加到每條read的開頭，用以拆分FASTQ文件。
下面是一個用於拆分fastq文件的示例，我們可以通過awk或python腳本輕鬆實現。

# 下載示例數據
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

# 查看fastq文件
zcat CEMBA180306_2B.demultiplexed.R1.fastq.gz | head 
@AGACGGAGACGAATCTAGGCTGGTTGCCTTAC:7001113:920:HJ55CBCX2:1:1108:1121:1892 1:N:0:0
ATCCTGGCATGAAAGGATTTTTTTTTTAGAAAATGAAATATATTTTAAAG
+
DDDDDIIIIHIIGHHHIIIHIIIIIIHHIIIIIIIIIIIIIIIIIIIIII

Step 2. Index reference genome (snaptools).

接下來，我們將對參考基因組構建索引用於後續的比對分析，這裡我們展示了如何使用BWA來構建參考基因組的索引。用戶可以通過—-aligner參數指定要使用的比對軟體，目前snaptools可以支持bwa, bowtie2和minimap2比對軟體。同時，我們還需要指定比對軟體所在文件夾的路徑，例如，如果bwa安裝在/opt/biotools/bwa/bin/bwa下，我們需要指定--path-to-aligner=/opt/biotools/bwa/bin/

# 查看bwa軟體所在的路徑
which bwa
/opt/biotools/bwa/bin/bwa 

# 使用snaptools構建參考基因組索引
snaptools index-genome  \
    --input-fasta=mm10.fa  \
    --output-prefix=mm10  \
    --aligner=bwa  \
    --path-to-aligner=/opt/biotools/bwa/bin/  \
    --num-threads=5

Step 3. Alignment (snaptools).

構建好參考基因組索引後，我們使用snaptools將拆分後的FASTQ reads序列比對到相應的參考基因組上。比對後，將比對好的bam文件按reads名稱進行排序。我們還可以通過設置--num-threads參數指定多個CPU加快比對的速度。

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

比對完之後，我們將比對好的雙端reads轉換為fragments片段，並查看每個 fragment片段的以下屬性:
1)比對質量評分MAPQ;
2)比對上的reads兩端是否按比對的flag值正確配對;
3)fragments片段的長度：
我們根據以下條件進行過濾篩選，只保留滿足條件的fragments片段;
1）兩端正確配對的fragments片段；
2）MAPQ值大於30的fragments片段(-min-mapq);
3）長度小於1000bp的fragments片段 (-max-flen)。

# 提取參考基因染色體長度信息
fetchChromSizes mm10.fa > mm10.chrom.sizes

# 使用snaptools進行數據預處理，生成snap格式文件
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=True

Step 5. Cell-by-bin matrix (snaptools).

最後，我們使用snap文件創建不同解析度下的cell-by-bin矩陣文件。在下面的示例中，我們分別創建了1,000、5,000和10,000 bin size下的三個cell-by-bin矩陣，並將所有的矩陣都存儲在cemba180306_2B.snap文件中。

# 使用snaptools創建cell-by-bin矩陣
snaptools snap-add-bmat    \
    --snap-file=CEMBA180306_2B.snap \
    --bin-size-list 1000 5000 10000 \
    --verbose=True

3）How to create snap file from 10X dataset?

Case 1

1）首先，運行cellranger-atac mkfastq生成拆庫後的fastq文件；
2）接下來，對於每個測序文庫，識別相應的測序文件，其中R1和R3是測序的reads，I1是16bp i5 (10x Barcode), R2是i7 (sample index)。
3）最後，使用snaptools提供的dex-fastq子程序，將10X barcode信息添加到reads的名稱中。

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.gz

Case 2 

在本示例中，我們有兩個10x的測序文庫(每個文庫都通過單獨的Chromium chip channel進行處理)。請注意，在運行完cellranger-atac mkfastq拆分數據之後，我們對每個文庫進行單獨的cellranger-atac count處理。

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.gz

4）Can I run SnapATAC with CellRanger outcome?

Yes. There are two entry points

（1）use the position sorted bam file (recommanded).

# 查看比對的bam文件信息
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:2

在比對的read名稱前添加barcode信息
The 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

# 查看添加好barcode信息的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

對bam文件按read的名稱進行排序
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.bam

使用snaptools進行數據預處理生成snap文件
Then 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=True

刪除中間文件
remove temporary files

rm atac_v1_adult_brain_fresh_5k.snap.snap.bam 
rm atac_v1_adult_brain_fresh_5k_possorted.header.sam

（2）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=True

5）How to create a snap file from bam or bed file?

在很多情況下，我們可以直接使用snaptools pre子程序將比對好的、按read名稱進行排序的bam或bed文件作為輸入，生成snap格式文件。強烈建議使用未經過濾的比對文件作為輸入。
（1）對於bam文件，我們需要在read的名稱前添加細胞的barcode信息，如下所示：

samtools view demo.bam|head

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將bam文件轉換為snap格式文件
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=True

（2）對於bed格式的文件，應將barcode信息添加到第四列中，如下所示：

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將bed文件轉換為snap格式文件
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=True

6）How to group reads?

（1）Group 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);

（2）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);

（3）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);

7）How to analyze multiple samples together?

因為SnapATAC軟體使用cell-by-bin矩陣對細胞進行聚類分群，這使他很容易將多個樣本進行結合併執行比較分析。它需要將所有的樣本創建相同bin size大小的cell-by-bin矩陣。在這裡，我們以PBMC_5K和PBMC_10K數據為例進行分析。

$ R
# 加載SnapATAC包
> 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");

# 使用createSnap函數將兩個數據結合在一起
> x.sp = createSnap(file=file.list, sample=sample.list);

createSnap函數將創建一個snap對象，該對象中包含了每個snap文件的名稱和相應的barcodes信息。

8）How to choose bin size?

SnapATAC軟體是基於cell-by-bin矩陣對細胞進行聚類分群的，因此選擇不同的bin size可能對細胞的聚類分群產生較大的影響。如何選擇最優的bin size，這個問題沒有絕對的答案。

一方面，我們發現在5kb-50kb範圍內的bin size大小的改變，並沒有明顯地改變細胞聚類分群的結果(如下圖所示)。而另一方面，我們確實注意到了一個大的bin size通常會生成相對較少的cluster。這種聚類的差異，可以使用解析度較小的Louvain聚類算法進行彌補。

使用較大bin size的好處是可以節省一些內存，這對於一些大型數據集尤為有用。這裡提供了一個bin size大小選擇的主觀建議。

參考來源：https://github.com/r3fang/SnapATAC/wiki/FAQs#what-is-a-snap-file

為了使更多的科研工作者快速做出高端大氣的科研圖表，並能夠加入到hiplot開發小組和成長計劃中。經團隊商量決定：從今日起，hiplot網站開放免費註冊，不再需要邀請碼，小夥伴們還在等什麼，快點註冊吧！

註冊網址：

https://hiplot.com.cn/signup

掃描下方二維碼快速註冊

科研繪圖神器—hiplot，是2020年7月19日openbiox聯合科研貓鄭重推出的全網首個開源繪圖平臺，目前提供基於R語言的60餘種基礎可視化和50餘種進階繪圖的功能，同時還部署了多個 openbiox社區項目（如bget下載文獻附錄、UCSCXenaShiny 等）。

截止目前，網站的總訪問量大約13萬餘次，在免註冊使用大部分可視化插件，以及僅開放教育用戶和邀請註冊的情況下，已有正式註冊用戶千餘人，日均訪問量千餘次。

https://hiplot.com.cn

點擊圖片進入Hiplot平臺介紹