SeekArc Tools v1.0.0
SeekArc Tools v1.0.0 是寻因生物开发的一套处理单细胞ATAC多组学数据的软件，该软件分别对转录组文库和ATAC文库进行处理，如 Barcode 提取与矫正，比对。对每个 Barcode 进行定量，进行联合判定细胞，得到下游分析的细胞表达矩阵，可进行下游聚类和差异分析，计算基因和 Peak 的连接关系。
软件下载
SeekArc Tools v1.0.0
Download-SeekArc Tools - md5: af6e054480ca082a554b499fe79c3a48
wget下载方式
mkdir seekarctools_v1.0.0
cd seekarctools_v1.0.0
wget -c -O seekarctools_v1.0.0.tar.gz "https://seekgene-public.oss-cn-beijing.aliyuncs.com/software/seekarctools/seekarctools_v1.0.0.tar.gz"curl下载方式
mkdir seekarctools_v1.0.0
cd seekarctools_v1.0.0
curl -C - -o seekarctools_v1.0.0.tar.gz "https://seekgene-public.oss-cn-beijing.aliyuncs.com/software/seekarctools/seekarctools_v1.0.0.tar.gz"软件安装
安装：
# decompress
tar zxf seekarctools_v1.0.0.tar.gz
# install
source ./bin/activate
./bin/conda-unpack
# export path in bashrc
export PATH=`pwd`/bin:$PATH
echo "export PATH=$(pwd)/bin:\$PATH" >> ~/.bashrc确认安装：
seekarctools_py --version使用教程
数据准备
测试数据下载
测试数据 - md5: 45a48572f5ffd00dd6e6d79e3ad495d4（物种：小鼠）
wget下载方式
wget -c -O demo_data.tar.gz "https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/rawdata/demo_data.tar.gz"
# decompress
tar -zxf demo_data.tar.gzcurl下载方式
curl -C - -o demo_data.tar.gz "https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/rawdata/demo_data.tar.gz"
# decompress
tar -zxf demo_data.tar.gz参考基因组下载
Download-human-reference-GRCh38 - md5: a9ce266ead3d89bc0aa77cfa646ca309
Download-mouse-reference-GRCm38 - md5: c4f2244f1d0560d5d3b4c71da06fd08f
wget下载方式
# human
wget -c -O refdata-arc-GRCh38.tar.gz "https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/reference/refdata-arc-GRCh38.tar.gz"
# mouse
wget -c -O refdata-arc-GRCm38.tar.gz "https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/reference/refdata-arc-GRCm38.tar.gz"
# decompress
tar -zxvf refdata-arc-GRCh38.tar.gzcurl下载方式
# human
curl -C - -o refdata-arc-GRCh38.tar.gz "https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/reference/refdata-arc-GRCh38.tar.gz"
# mouse
curl -C - -o refdata-arc-GRCm38.tar.gz "https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/reference/refdata-arc-GRCm38.tar.gz"
# decompress
tar -zxvf refdata-arc-GRCh38.tar.gz参考基因组的构建可以参考如何构建参考基因组？
运行 SeekArc Tools v1.0.0
运行示例
示例1：基本用法
为分析设置必要的配置文件，包括样本数据的路径、参考基因组路径等。使用以下命令运行 SeekArc Tools v1.0.0：
seekarctools_py arc run \
--rnafq1 /path/to/demo/demo_GE_S1_L001_R1_001.fastq.gz \
--rnafq2 /path/to/demo/demo_GE_S1_L001_R2_001.fastq.gz \
--atacfq1 /path/to/demo/demo_arc_S2_L002_R1_001.fastq.gz \
--atacfq2 /path/to/demo/demo_arc_S2_L002_R2_001.fastq.gz \
--samplename demo \
--outdir /path/to/outdir \
--refpath /path/to/reference/GRCh38 \
--include-introns \
--core 16示例2：设置阈值重新进行 call peaks 或 call cell
若觉得默认参数下判定出的 peaks 或 cells 不符合需求，可调整参数跳过比对等之前的步骤，重新运行以节省时间
seekarctools_py arc retry \
--samplename demo \
--outdir /path/to/outdir \
--refpath /path/to/reference/GRCh38 \
--core 16 \
--qvalue 0.01 \
--snapshift 0 \
--extsize 200 \
--min_len 200 \
--min_atac_count 1000 \
--min_gex_count 500NOTE
需注意保持之前运行的文件没有删除或移动。--outdir为SeekArc Tools v1.0.0第一次分析后的目录路径。--min_atac_count需和--min_gex_count一起使用，否则不生效。
参数说明
| 参数 | 参数说明 |
|---|---|
| --rnafq1 | 表达文库 R1 fastq数据路径。 |
| --rnafq2 | 表达文库 R2 fastq数据路径。 |
| --atacfq1 | ATAC文库 R1 fastq数据路径。 |
| --atacfq2 | ATAC文库 R2 fastq数据路径。 |
| --samplename | 样本名称。 |
| --outdir | 结果输出目录。默认值：./。 |
| --skip_misB | barcode不允许碱基错配，默认允许一个碱基错配。 |
| --skip_misL | linker不允许碱基错配，默认允许一个碱基错配。 |
| --skip_multi | 舍弃能纠错为多个白名单barocde的reads，默认纠错为比例最高的barcode。 |
| --skip_len | 在适配器筛选后跳过筛选短读取，将使用短读取。 |
| --core | 分析使用的线程数。 |
| --include-introns | 不启用时，只会选择exon reads⽤于定量；启用时，intron reads也会⽤于定量。 |
| --refpath | 参考基因组路径。 |
| --star_path | 外部 STAR 软件路径。如果参考基因组中的索引由其他STAR构建，请指定该路径。 |
| --qvalue | 峰值检测的最小FDR（q值）截止值。默认值：0.05。 |
| --nolambda | 若启用则 MACS3 将使用背景lambda作为本地lambda。这意味着macs3不会考虑峰值候选区域的局部偏差。 |
| --snapshift | MACS3峰值检测的移动幅度。默认值：0。 |
| --extsize | MACS3峰值检测的拓展幅度。默认值：400。 |
| --min_len | 峰的最小长度。如果无，则设置为"extsize"。默认值：400。 |
| --broad | 若启用则进行广泛的峰值调用，以UCSC gappedPeak格式生成结果，该格式封装了峰值的嵌套结构。 |
| --broad_cutoff | 广泛的峰值调用的阈值。默认值：0.1。 |
| --min_atac_count | 定义细胞条形码峰中ATAC转置事件的最小数量（ATAC计数）。 |
| --min_gex_count | 定义细胞条形码中UMI的最小数量（GEX计数）。 |
| -h，--help | 显示参数说明 |
结果文件说明
以下是输出目录结构：每一行代表一个文件或文件夹，用 "├──" 表示，数字表示重要的输出文件。
./
├── outs
│ ├── filtered_feature_bc_matrix 1
│ │ ├── barcodes.tsv.gz
│ │ ├── features.tsv.gz
│ │ └── matrix.mtx.gz
│ ├── filtered_peaks_bc_matrix 2
│ │ ├── barcodes.tsv.gz
│ │ ├── features.tsv.gz
│ │ └── matrix.mtx.gz
│ ├── raw_feature_bc_matrix 3
│ │ ├── barcodes.tsv.gz
│ │ ├── features.tsv.gz
│ │ └── matrix.mtx.gz
│ ├── raw_peaks_bc_matrix 4
│ │ ├── barcodes.tsv.gz
│ │ ├── features.tsv.gz
│ │ └── matrix.mtx.gz
│ ├── demo_A_fragments.tsv.gz 5
│ ├── demo_A_fragments.tsv.gz.tbi 6
│ ├── demo_A_mem_pe_Sort.bam 7
│ ├── demo_A_peaks.bed 8
│ ├── demo_E_SortedByName.bam 9
│ ├── demo.rds 10
│ ├── demo_report.html 11
│ └── demo_summary.csv 12
└── analysis- 转录组细胞稀疏矩阵
- ATAC 细胞稀疏矩阵
- 转录组原始稀疏矩阵
- ATAC 原始稀疏矩阵
- fragments 文件
- fragments 文件索引
- ATAC 的 bam 文件
- peaks 文件
- 转录组的 bam 文件
- rds格式的样本结果
- 样本的html报告
- 样本的csv格式质控信息
处理过程
表达文库分析
estep1: barcode/UMI提取
转录组文库结构： 
barcode 处理：
根据结构设计，取出相应序列。当取出的barcode序列在白名单中时，我们认为它是有效barcode，计入有效barcode的reads数量；当barcode不在白名单中时，我们认为它是无效barcode。
TIP
测序过程中，有一定几率发生测序错误。在提供有白名单的情况下，SeekArc Tools v1.0.0可以尝试barcode矫正。在启用矫正时，当无效barcode一个碱基错配（一个hamming distance）的序列存在于白名单中：
- 只有唯一一个序列存在于白名单中：我们将这个无效barcode矫正为白名单中barcode
- 有多个序列存在于白名单中：我们将这个无效barcode改为read支持数量最多的序列
接头处理：
在转录组中，Read2的末端有可能会出现 7F 的反向互补序列，建库时可能引入的接头序列。我们对这些污染序列进行切除，剪切完的read2长度需要大于设定的最小长度来保证有足够的信息，准确比对到基因组的位置。如果剪切完成后的read2长度小于最小长度，我们认为这条read为无效read。
estep2: 进行比对并找到比对基因
序列比对
- SeekArc Tools v1.0.0使用STAR比对软件将处理后的R2比对到参考基因组上。
- SeekArc Tools v1.0.0使用QualiMap软件和转录本注释文件GTF，统计reads比对外显子、内含子和基因间区等的比例。
- 使用featureCounts将比对上的read注释到基因上，可以选择不同的注释规则，如链方向性和定量的feature。当使用外显子定量时，当read超过50%碱基比对到外显子区域时，认为该read来源于此外显子以及外显子对应的基因；当使用转录本定量时，当read超过50%碱基比对到转录本区域时，认为该read来源于此转录本以及转录本对应的基因。
经过上述处理后，有如下数据指标：
- Reads Mapped to Genome: 比对到参考基因组上的reads占所有reads的比例（包括只有一个比对位置和多个比对位置的reads）
- Reads Mapped Confidently to Genome: 在参考基因组上只有一个比对位置的reads占所有reads的比例
- Reads Mapped to Intergenic Regions：比对到基因间隔区reads占所有reads的比例
- Reads Mapped to Intronic Regions：比对到内含子reads占所有reads的比例
- Reads Mapped to Exonic Regions：比对到外显子reads占所有reads的比例
estep3: 定量
UMI定量：
SeekArc Tools v1.0.0以barcode为单位提取featureCounts输出的bam数据，统计注释到基因的UMI和UMI对应的reads数：
- 过滤掉对应UMI为单个重复碱基的reads， 例如UMI为TTTTTTTT
- 当read注释到多个基因，在有唯一外显子的注释时，为有效read，其他情形下都过滤掉
UMI矫正：
NOTE
测序过程中，UMI也有一定概几率出现测序错误。SeekArc Tools v1.0.0默认使用UMI-tools的adjacency方法对UMI进行矫正。
Image source: https://umi-tools.readthedocs.io/en/latest/the_methods.html
ATAC 文库分析
astep1: barcode 提取
ATAC 文库结构： 
barcode 处理：
同estep1的处理方式一样。
接头和 reads 处理：
在ATAC 文库中，Read1 带有固定序列，末端有可能会出现接头序列；Read2的末端有可能会出现 固定序列的反向互补序列和接头序列。我们对这些污染序列进行切除，剪切完的read1或read2长度需要大于设定的最小长度来保证有足够的信息，准确比对到基因组的位置。如果剪切完成后的read1或read2长度小于最小长度，我们认为这条read为无效read；如果剪切完成后的read1或read2长度大于50，则截取read1或read2序列的前50bp；如果剪切完成后的read1或read2长度小于50且大于最小长度，则保留read1或read2的全部序列。
astep2: 进行比对
序列比对：
使用 bwa mem 将处理后的read1和read2比对到参考基因组上。
astep3: 定量
使用 SnapATAC2 v2.8.0 处理 astep2 输出的 bam 文件。
- 读取 bam ，获取ATAC 的 Sequencing 和 Mapping 相关指标。
- 由 bam 生成原始 fragments 文件。过滤掉线粒体相关的 fragments 后得到过滤后的 fragments 文件。
- 由过滤后的 fragments 文件生成 Anndata 对象。
- 绘制 TSS 得分分布图。
- 使用 MACS3 进行 call peak，并过滤重复的 peaks，输出到 peaks.bed。
- 根据 peaks 生成 raw_peaks_bc_matrix。
- 统计每个 barcode 的相关指标，输出到 per_barcode_metrics.csv。
- 绘制插入片段分布图。
- 细胞判定
- 生成 filtered_feature_bc_matrix 和 filtered_peaks_bc_matrix。
- 统计ATAC 的 Cells 和 Targeting 相关指标。
细胞判定
IMPORTANT
若指定 --min_atac_count 和 --min_gex_count：
- 则判定 barcode 中 UMI 数量大于 min_gex_count 值且 events 数量大于 min_atac_count 值的barcode为细胞。
若不指定：
- 首先，保留 barcode 中 UMI 数量大于1 或 events in peaks 数量大于 1的 barcode。
- 然后，过滤 "the fraction of fragments overlapping called peaks" 小于 "the fraction of genome in peaks" 的 barcode。
- 去除重复值。将具有相同 UMI 值和 events in peaks 值的 barcode 过滤。
- 将过滤后的 barcode 的UMI 值和 events in peaks 值进行 log10 转换。
- 根据转换后的UMI 值和 events in peaks 值，使用 KMeans 算法将 barcode 分为两群其中均值高的群为细胞，另一群则不是细胞。
- 将去重过滤掉的 barcode 通过判定的 barcode 映射为细胞或非细胞。
astep4: 后续分析
经过上述步骤，得到细胞的表达矩阵和ATAC 矩阵后，我们可以进行下一步的分析。
使用 Signac v1.14.0 将细胞的表达矩阵和ATAC 矩阵合并为一个 seurat object，过滤不在数据库中的 feature，分别对表达矩阵和ATAC 矩阵进行归一化、寻找高变基因、降维聚类之后寻找差异基因，进行 WNN 分群，最后计算 feature link。
发布说明
v1.0.0
首次发布
- 发布稳定版本，支持ATAC多组学数据标准分析流程