SeekOne DD ATAC+GEX Pipeline 用户指南

一、获取数据¶

本次分析示例数据来源 新鲜小鼠脑组织，input细胞量为 8000，ATAC测序数据量 80G，RNA测序数据量140G。本次实验使用的 SeekOne DD ATAC+转录组多组学试剂盒

原始数据下载链接：
https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/rawdata/XYRD-WTSH814-Rawdata.tar
md5sum: c35c347daaf11a03a65bf91e63de01d1
脚本及文件下载链接：
https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/sgrun_demo.sh
https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/code.tar
https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/file.tar
注意：file 目录下的 report_arc.py 和 reportarc 目录需要在同一目录下。

二、参考文件下载¶

Human reference (GRCh38) 下载链接：
https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/reference/refdata-arc-GRCh38.tar.gz
md5sum: 279a3976084d8d32898a0fefafe11d75

Mouse reference (GRCm38) 下载链接：
https://seekgene-public.oss-cn-beijing.aliyuncs.com/sgarc_demo/reference/refdata-arc-GRCm38.tar.gz
md5sum: 3f96e76b93c7dd21c11931b324602249

三、环境准备¶

  1. 请先检查并安装以下软件： seeksoultools 1.2.2
    参考安装网址：http://seeksoul.seekgene.com/zh/v1.2.2/index.html
    fastp 0.23.2
    参考安装网址：https://github.com/OpenGene/fastp
    samtools 1.16.1
    参考安装网址：https://github.com/samtools/samtools/releases
    bwa 0.7.17-r1188
    参考安装网址：https://github.com/lh3/bwa
    bedtools v2.26.0
    参考安装网址：https://bedtools.readthedocs.io/en/latest/content/installation.html
    bgzip 1.16
    参考安装网址：https://github.com/samtools/bcftools/releases
    tabix 1.16
    参考安装网址：https://anaconda.org/bioconda/tabix
    gunzip 1.5
  2. 请先检查并安装以下 python 包：
    click
    matplotlib.pyplot
    scipy.io
    snapatac2
  3. 请先检查并安装以下 R 包：
    BiocParallel
    future
    Signac

四、参数设置¶

In [ ]:
outdir='output dir'
gsample='XYRD-WTSH814-E'
gfq1='XYRD-WTSH814-E_S6_L006_R1_001.fastq.gz'
gfq2='XYRD-WTSH814-E_S6_L006_R2_001.fastq.gz'
asample='XYRD-WTSH814-A'
afq1='XYRD-WTSH814-A_S5_L006_R1_001.fastq.gz'
afq2='XYRD-WTSH814-A_S5_L006_R2_001.fastq.gz'

core='16' # 设置核数
memory='60' # 设置内存使用量
bedtoolspath='bedtools path'

species='human or mouse'
refpath='reference path'
genomefa='reference path/fasta/genome.fa'
genomeDir='reference path/star'
gtf='reference path/genes/genes.gtf'
annords='file/Anno_EnsDb_Mmusculus_v79.rds'

五、数据质控¶

1. 使用 fastp 对 GEX 和 ATAC 文库进行质控，去除低质量的 reads，得到文库质量报告¶

输入：
Rawdata
输出：
fastp 目录

In [ ]:
## 
# 对 gex 文库进行质控
mkdir -p ${outdir}/fastp/${gsample}
fastp \
    -i ${outdir}/Rawdata/${gfq1} \
    -I ${outdir}/Rawdata/${gfq2} \
    -o ${outdir}/fastp/${gsample}/${gfq1} \
    -O ${outdir}/fastp/${gsample}/${gfq2} \
    --reads_to_process 0 --cut_front --cut_front_window_size 4 --cut_front_mean_quality 10 --cut_tail --cut_tail_window_size 1 --cut_tail_mean_quality 3 \
    -j ${outdir}/fastp/${gsample}/${gsample}_fastp.json \
    -h ${outdir}/fastp/${gsample}/${gsample}_fastp.html \
    --thread ${core}

# 对 atac 文库进行质控
mkdir -p ${outdir}/fastp/${asample}
fastp \
    -i ${outdir}/Rawdata/${afq1} \
    -I ${outdir}/Rawdata/${afq2} \
    -o ${outdir}/fastp/${asample}/${afq1} \
    -O ${outdir}/fastp/${asample}/${afq2} \
    --reads_to_process 0 --cut_front --cut_front_window_size 4 --cut_front_mean_quality 10 --cut_tail --cut_tail_window_size 1 --cut_tail_mean_quality 3 \
    -j ${outdir}/fastp/${asample}/${asample}_fastp.json \
    -h ${outdir}/fastp/${asample}/${asample}_fastp.html \
    --thread ${core}

2. 使用 SeekSoulTools 分析 转录组文库数据¶

在这一步，我们使用 SeekSoulTools 分析转录组数据，并得到转录组质控结果。
输入：
fastp/fastq
输出：
outdir

In [ ]:
seeksoultools rna run \
    --fq1 ${outdir}/fastp/${gsample}/${gfq1} \
    --fq2 ${outdir}/fastp/${gsample}/${gfq2} \
    --samplename ${gsample} \
    --outdir ${outdir} \
    --genomeDir ${genomeDir} \
    --gtf ${gtf} \
    --chemistry DD-Q \
    --include-introns \
    --core ${core}

3. 分析 ATAC 文库数据¶

step1：数据预处理¶

① 提取 barcode 和 umi 信息，去除接头序列¶

ATAC 文库构成如下：
wenku.png 在这一步，我们将 barcode 和 umi 信息提取保存到 reads id 中，将barcode 相关统计信息写入 summary.json 文件中，并识别和去除可能测到的接头序列。
输入：
fastp/fastq
输出：
step1/asample_1.fq.gz
step1/asample_2.fq.gz
asample_summary.json

In [ ]:
mkdir -p ${outdir}/${asample}
python /code/barcode.py \
	--fq1 ${outdir}/fastp/${asample}/${afq1} \
	--fq2 ${outdir}/fastp/${asample}/${afq2} \
	--samplename ${asample} \
	--outdir ${outdir}/${asample} \
	--barcode /file/P3CB.barcode.txt.gz \
	--chemistry DD-AG \
	--core ${core}

② 对 step1 处理后的 reads 进行剪切，使剪切后的 R1,R2 仅为插入片段序列¶

R1剪切处理：

  1. 若 R1 序列小于 50 bp 则丢弃；
  2. 匹配 ‘CGTCCGTCGTTGCTCGTAGATGTGTATAAGAGACAG’ 固定序列在 R1 的位置，截取固定序列后的 50bp 序列，若不足 50bp 则全部保留；
  3. 若未匹配到固定序列则截取 43-93bp，若小于 93bp 则 43bp 后所有序列

R2剪切处理：

  1. 若 R2 序列长度 ≥ 50bp，则只保留前 50bp 序列，若不足 50bp，则全部保留

输入：
step1/asample_1.fq.gz
step1/asample_2.fq.gz
asample_summary.json
输出：
cutdata/asample_cutR1.fastq.gz
cutdata/asample_cutR2.fastq.gz
asample_summary.json

In [ ]:
mkdir -p ${outdir}/${asample}/cutdata
python /code/cut_reads.py \
    --afq1 ${outdir}/${asample}/step1/${asample}_1.fq.gz \
    --afq2 ${outdir}/${asample}/step1/${asample}_2.fq.gz \
    --step1json ${outdir}/${asample}/${asample}_summary.json \
    --outdir ${outdir}/${asample}/cutdata \
    --samplename ${asample}

step2：将剪切后的 reads 使用 bwa mem 比对并使用 samtools 排序¶

在这一步，我们将剪切后的 cutdata 目录下的 R1R2 文件，通过 bwa mem 进行双端比对，并使用 samtools 进行排序，输出排序后的 bam 文件
输入：
cutdata/asample_cutR1.fastq.gz
cutdata/asample_cutR2.fastq.gz
输出：
asample_mem_pe_Sort.bam

In [ ]:
mkdir -p ${outdir}/${asample}/step2/bwa_pe
cd ${outdir}/${asample}/step2/bwa_pe
bwa mem \
    -t ${core} -M -R "@RG\tID:$asample\tLB:WGS\tPL:Illumina\tPU:$asample\tSM:$asample" ${genomefa} \
    ${outdir}/${asample}/cutdata/${asample}_cutR1.fastq.gz \
    ${outdir}/${asample}/cutdata/${asample}_cutR2.fastq.gz | samtools sort -@ ${core} -o ${asample}_mem_pe_Sort.bam

step3：使用 snapatac2 分析 bam，输出统计信息和文件。¶

将 bam 文件通过 snapatac2 分析 atac 相关质控指标，使用 gex 分析得到的 cell barcode 与 snapatac2 分析得到的 cell barcode 取交集，联合得到 cell barcode并对 atac adata 和矩阵进行过滤，输出文件保存至目录 asample/step3
qvalue、shift、extsize、min_len 参数设置参考 macs3（https://macs3-project.github.io/MACS/docs/callpeak.html）
输入：
asample_mem_pe_Sort.bam
gsample/step3/counts.xls
gsample/step3/detail.xls
gsample_summary.json
asample_summary.json
输出：
filter_peaks_bc_matrix/ （过滤后的cell-peaks 矩阵）
frag_counts.xls （每个 barcode 中 fragments 及 reads 计数文件）
per_barcode_metrics.csv （每个 barcode 中 各指标计数文件）
raw_peaks_bc_matrix/ （原始的cell-peaks 矩阵）
asample_fragments_size.png （fragments 长度分布图）
asample_fragments.tsv.gz （fragments 文件）
asample_fragments.tsv.gz.tbi （fragments 文件索引）
asample_peaks.bed （peaks 文件）
asample_snapatac2_filter.h5ad （过滤后的 adata 对象）
asample_snapatac2_raw.h5ad （原始 adata 对象）
asample_tss_enrichment.png （TSS 富集分布图）

In [ ]:
mkdir -p ${outdir}/${asample}/step3
python /code/snaprun.py \
    --bam ${outdir}/${asample}/step2/bwa_pe/${asample}_mem_pe_Sort.bam \
    --atacjson ${outdir}/${asample}/${asample}_summary.json \
    --gexjson ${outdir}/${gsample}/${gsample}_summary.json \
    --outdir ${outdir}/${asample}/step3 \
    --samplename ${asample} \
    --countxls ${outdir}/${gsample}/step3/counts.xls \
    --detailxls ${outdir}/${gsample}/step3/detail.xls \
    --species ${species} \
    --refpath ${refpath} \
    --bedtoolspath ${bedtoolspath} \
    --core ${core} \
    --qvalue 0.05 \
    --shift 0 \
    --extsize 400 \
    --min_len 400

对 step3 生成的 asample_fragments.tsv.gz 文件进行排序，并建立索引¶

In [ ]:
cd ${outdir}/${asample}/step3 &&\
gunzip ${asample}_fragments.tsv.gz &&\
bedtools sort -i ${asample}_fragments.tsv > ${asample}_fragments_sort.tsv &&\
bgzip -c ${asample}_fragments_sort.tsv > ${asample}_fragments.tsv.gz &&\
tabix -p bed ${asample}_fragments.tsv.gz &&\
rm ${asample}_fragments.tsv ${asample}_fragments_sort.tsv

step4：使用 Signac 和 Seurat 对 gex 和 atac 的矩阵文件进行下游分析，并计算 links 信息¶

输入：
gsample/step3/raw_feature_bc_matrix
asample/step3/filter_peaks_bc_matrix
asample_fragments.tsv.gz
annords
输出：
filter_peaks_bc_matrix.rds （过滤后的 cell-peaks 矩阵）
linked_feature.xls （计算的 gene 与 peaks 的 links 关系文件）
gex_FindAllMarkers.xls （gex 差异基因文件）
gex_tsne_umi.xls （gex 分群信息）
atac_tsne_umi.xls （atac 分群信息）
joint_peak_link_gene.rds （计算 links 后的 seurat 对象）

In [ ]:
mkdir -p ${outdir}/${asample}/step4
Rscript /code/count_link.R \
    --gex_matrix_dir ${outdir}/${gsample}/step3/raw_feature_bc_matrix \
    --atac_matrix_dir ${outdir}/${asample}/step3/filter_peaks_bc_matrix \
    --fragpath ${outdir}/${asample}/step3/${asample}_fragments.tsv.gz \
    --outdir ${outdir}/${asample}/step4 \
    --species ${species} \
    --anno_rds ${annords} \
    --core ${core} \
    --memory ${memory}

4. 汇总 gex 和 atac 分析结果，输出质控报告¶

输入：
asample_summary.json
gsample_summary.json
输出：
asample_all_summary.csv
asample_report.html

In [ ]:
python /file/report_arc.py \
    --atacjson ${outdir}/${asample}/${asample}_summary.json \
    --gexjson ${outdir}/${gsample}/${gsample}_summary.json \
    --outdir ${outdir} \
    --samplename ${asample}