寻因单细胞转录组+甲基化双组学分析教程
1. 背景介绍
单细胞转录组+甲基化双组学技术是当前单细胞多组学分析的重要发展方向，能够同时获取单个细胞的基因表达信息和DNA甲基化信息，为理解细胞异质性、发育轨迹和表观遗传调控机制提供重要工具。
1.1 技术原理与优势
单细胞转录组+甲基化双组学技术通过特殊的文库构建方法，在同一个细胞中同时捕获RNA和DNA，实现：
- 转录组分析：检测基因表达水平，反映细胞的功能状态
- 甲基化分析：检测DNA甲基化水平，反映表观遗传调控状态
- 整合分析：将两种数据类型结合，揭示基因表达与表观遗传调控的关系
IMPORTANT
单细胞转录组+甲基化双组学分析能够揭示传统单组学分析无法发现的生物学机制，如基因表达与DNA甲基化的动态关系、细胞分化过程中的表观遗传重编程等。
1.2 应用价值
该技术在多个研究领域具有重要应用价值：
- 发育生物学：追踪胚胎发育过程中的基因表达和表观遗传变化
- 干细胞研究：研究干细胞分化过程中的表观遗传重编程
- 癌症研究：发现肿瘤异质性的表观遗传基础
- 神经科学：理解神经元分化和功能特化的表观遗传机制
- 免疫学：研究免疫细胞分化和激活的表观遗传调控
2. 流程介绍
寻因单细胞转录组+甲基化双组学流程适用于分析北京寻因生物有限公司（以下简称寻因）单细胞转录组+甲基化双组学试剂盒（以下简称试剂盒）产品的分析流程。其中，单细胞转录组分析使用寻因开发的SeekSoul Tools v1.2.2软件。甲基化分析则基于bismark和ALLCools等软件，能够做到检测单个细胞的全基因的甲基化水平，构建特定大小bins内的平均甲基化水平矩阵，并根据矩阵进行分群聚类等下游分析。
NOTE
本教程基于寻因生物的单细胞转录组+甲基化双组学试剂盒，采用特定的文库结构和分析流程。使用其他平台的数据可能需要调整相应的参数和流程。
3. 软件下载及安装
3.1 软件下载及解压
TIP
建议在下载前检查网络连接和磁盘空间，确保有足够的存储空间用于软件和数据文件。
# 复制粘贴下面命令到shell并运行。将对软件进行下载和解压。解压后可以直接运行，无需安装。
mkdir methy_pipe
cd methy_pipe
wget -dc -O database.tar.gz "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Fdatabase.tar.gz?Expires=2110214572&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=LQGg30uNtrobYWL9%2BE%2FDGfXT4Ig%3D"
wget -dc -O database.tar.gz.md5 "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Fdatabase.tar.gz.md5?Expires=2110214656&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=EG9ZBDhfEZzp3IeIxJ1w1Y2tqak%3D"
wget -dc -O script.tar.gz "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Fscript.tar.gz?Expires=2110214677&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=SNy%2FVOJVdC%2FLTq1PiVVV93BAUYs%3D"
wget -dc -O script.tar.gz.md5 "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Fscript.tar.gz.md5?Expires=2110214688&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=mScGtJ6%2BHLnAbKz5M90fBm9bvds%3D"
wget -dc -O seeksoultools.1.2.2.tar.gz "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Fseeksoultools.1.2.2.tar.gz?Expires=2110214707&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=AF5BhIkxXcz0E29eRX8QEXglgkY%3D"
wget -dc -O seeksoultools.1.2.2.tar.gz.md5 "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Fseeksoultools.1.2.2.tar.gz.md5?Expires=2110214718&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=Qzd8c1zY25C41PIo8AKMZ71sUW0%3D"
# 注意：解压后，methy_pipe目录下应该包含database，seeksoultools.1.2.2和script等3个子目录
tar -xzf database.tar.gz
mkdir -p seeksoultools.1.2.2
tar -xzf seeksoultools.1.2.2.tar.gz -C seeksoultools.1.2.2
tar -xzf script.tar.gzIMPORTANT
解压完成后，请确认methy_pipe目录下包含database、seeksoultools.1.2.2和script三个子目录，这是后续分析的必要组件。
4. Demo数据下载
NOTE
本教程使用样本XYRD-WTJW880作为演示数据，该数据来源于人的PBMC，单细胞转录组的测序数据量为3G，单细胞甲基化的数据量为500G。
4.1 单细胞转录组数据下载
# 单细胞转录组数据
wget -dc -O XYRD-WTJW880-E_S1_L005_R1_001.fastq.gz "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Ffastq%2FXYRD-WTJW880-E_S1_L005_R1_001.fastq.gz?Expires=2110217167&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=%2BH6wD2ok1tZ%2FcSVdaV5xKilSC0Q%3D"
wget -dc -O XYRD-WTJW880-E_S1_L005_R2_001.fastq.gz "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Ffastq%2FXYRD-WTJW880-E_S1_L005_R2_001.fastq.gz?Expires=2110217177&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=%2BoCn07K6WKht9yArdEDYuB7EG14%3D"
wget -dc -O XYRD-WTJW880-E_S1_L005_R1_001.fastq.gz.md5 "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Ffastq%2FXYRD-WTJW880-E_S1_L005_R1_001.fastq.gz.md5?Expires=2110217198&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=0DY4vPJPqDweSEQjPTpp5neNwWs%3D"
wget -dc -O XYRD-WTJW880-E_S1_L005_R2_001.fastq.gz.md5 "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Ffastq%2FXYRD-WTJW880-E_S1_L005_R2_001.fastq.gz.md5?Expires=2110217188&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=sS8XYhQNjvyB%2Bx%2Bps7sulBs1rpM%3D"4.2 单细胞甲基化数据下载
# 单细胞甲基化数据
wget -dc -O XYRD-WTJW880-MET_S04_L004_R1_001.fastq.gz "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Ffastq%2FXYRD-WTJW880-MET_S04_L004_R1_001.fastq.gz?Expires=2110214755&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=ooflvzELvsOvz5vrJhlilPUoxPk%3D"
wget -dc -O XYRD-WTJW880-MET_S04_L004_R2_001.fastq.gz "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Ffastq%2FXYRD-WTJW880-MET_S04_L004_R2_001.fastq.gz?Expires=2110214765&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=QrvJ1llxxr8J9djIwvxyapac9Fw%3D"
wget -dc -O XYRD-WTJW880-MET_S04_L004_R1_001.fastq.gz.md5 "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Ffastq%2FXYRD-WTJW880-MET_S04_L004_R1_001.fastq.gz.md5?Expires=2110215374&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=Yprqhrj8VzOTHyhUVXMETsrO%2Frc%3D"
wget -dc -O XYRD-WTJW880-MET_S04_L004_R2_001.fastq.gz.md5 "http://seekgene-public.oss-cn-beijing.aliyuncs.com/methy_demo%2Fmethy_exp%2Ffastq%2FXYRD-WTJW880-MET_S04_L004_R2_001.fastq.gz.md5?Expires=2110215405&OSSAccessKeyId=LTAIxk81dBP5kdTu&Signature=bbTOujzXu9zaoBPEdzmzvmTzJeo%3D"5. 软件使用
5.1 使用方法
使用sc_methy_workflow.sh脚本完成整个分析，其输入包含7个参数，依次为：样本名称，单细胞转录组Read1的fastq文件，单细胞转录组Read2的fastq文件，单细胞甲基化Read1的fastq文件，单细胞甲基化Read2的fastq文件，结果输出路径，cpu数量。
# methy_pipe修改为你软件解压后的实际路径
sh ${methy_pipe}/script/sc_methy_workflow.sh \
XYRD-WTJW880 \
XYRD-WTJW880-E_S1_L005_R1_001.fastq.gz \
XYRD-WTJW880-E_S1_L005_R2_001.fastq.gz \
XYRD-WTJW880-MET_S04_L004_R1_001.fastq.gz \
XYRD-WTJW880-MET_S04_L004_R2_001.fastq.gz \
./ \
645.2 系统资源要求
IMPORTANT
单细胞转录组+甲基化双组学分析对计算资源要求较高，请确保满足以下要求：
- CPU要求：建议使用64核CPU，最低要求32核
- 内存要求：建议128GB内存
- 存储空间：根据数据量大小，建议预留足够的磁盘空间
5.3 运行时间预估
在CPU=64，内存=128G的机器上，XYRD-WTJW880样本的运行时间约为4天。
WARNING
实际运行时间可能因数据质量、机器性能、网络状况等因素而有所不同。建议在开始分析前预留充足的时间。
6. 软件代码说明文档
完整流程代码可查看${methy_pipe}/script/sc_methy_workflow.sh，下面对部分代码进行解释说明。
6.1 数据质控
使用fastp软件对raw data进行质控，过滤掉质量低的reads，得到clean data。
NOTE
fastp具体参数解释，请使用--help或详见fastp官网。
# fastp and rmdup
mkdir -p ${outdir}/fastp/${sample}
exp_fq1_clean=${outdir}/fastp/$(basename ${exp_fq1})
exp_fq2_clean=${outdir}/fastp/$(basename ${exp_fq2})
fastp \
-i ${exp_fq1_local} \
-I ${exp_fq2_local} \
-o ${exp_fq1_clean} \
-O ${exp_fq2_clean} \
--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/$(basename ${exp_fq1})_fastp.json \
-h ${outdir}/fastp/$(basename ${exp_fq1})_fastp.html \
--length_required 40 --thread ${core} -D
methy_fq1_clean=${outdir}/fastp/$(basename ${methy_fq1})
methy_fq2_clean=${outdir}/fastp/$(basename ${methy_fq2})
fastp \
-i ${methy_fq1_local} \
-I ${methy_fq2_local} \
-o ${methy_fq1_clean} \
-O ${methy_fq2_clean} \
--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/$(basename ${methy_fq1})_fastp.json \
-h ${outdir}/fastp/$(basename ${methy_fq1})_fastp.html \
--length_required 40 --thread ${core} -D6.2 转录组文库分析
使用SeekSoul Tools v1.2.2对转录组文库进行分析，后续甲基化文库的细胞基于转录组文库判定的细胞barcode。
exp_outdir=${outdir}/${sample}_exp
mkdir -p ${exp_outdir}
seeksoultools rna run \
--fq1 ${exp_fq1_clean} \
--fq2 ${exp_fq2_clean} \
--samplename ${sample} \
--genomeDir ${genomeDir} \
--gtf $gtf \
--chemistry DDV2 \
--core ${core} \
--include-introns \
--outdir ${exp_outdir}
gexcb=$exp_outdir/${sample}/step3/filtered_feature_bc_matrix/barcodes.tsv.gz6.3 甲基化文库分析
6.3.1 Step1：Barcode提取和CT转化率统计
根据Read1的结构设计和参数对barcode/UMI进行提取和处理，统计C-T转化率，对Read1和Read2进行过滤，输出新的fastq文件。
输入：去除了低质量reads的fq1，fq2
输出：step1目录
文库结构
结构说明：
- SP1/SP2：接头序列
- 7F：连接序列
- 17L和ME：固定序列
- 9bp：为插入片段的延伸序列
- B：barcode碱基
- U：UMI碱基
barcode和UMI描述以字母和数字描述Read1的基本结构，字母描述碱基含义，数字描述碱基长度。
Barcode识别和矫正
根据结构设计，确定barcode所在序列位置，取出相应序列。当取出的barcode序列在白名单中时，我们认为它是有效barcode，计入有效barcode的reads数量；当barcode不在白名单中时，我们认为它是无效barcode。
测序过程中，有一定几率发生测序错误。在提供有白名单的情况下，SeekSoul Tools可以尝试barcode矫正。在启用矫正时，当无效barcode一个碱基错配（一个hamming distance）的序列存在于白名单中：
- 只有唯一一个序列存在于白名单中：我们将这个无效barcode矫正为白名单中barcode
- 有多个序列存在于白名单中：我们将这个无效barcode改为read支持数量最多的序列
统计CT转化率
每条R1上理论上都包含下面的A或者B序列之一，A和B的结构从前到后依次是7L序列（7bp），17L序列（17bp）和ME序列（19bp）。其中，A序列代表非甲基化C碱基转化为T碱基的情况，B序列代表甲基化C的情况。
A：TTGTTGTTGTTTGTTGTTGTTTGTAGATGTGTATAAGAGATAG
B：TTGCTGTCGTCCGTCGTTGCTCGTAGATGTGTATAAGAGACAG
为了计算CT转化率，我们首先需要筛选出A序列，条件包括：包含正确的7F序列，第22和41位碱基是T（区分A和B）。然后统计A的第8,11,12,15和20位，共5个位置上的"T"碱基的数量。
CT转化率计算公式：
# 提取具有有效barcode的reads，剪切接头。并将reads和统计信息写入outdir/sample/sample_summary.json中。
methy_dir=${outdir}/${sample}_methy
########################## step1 ##########################
mkdir -p ${methy_dir}
mkdir -p ${methy_dir}/step1
python ${script_path}/barcode_cs_multi.py \
--fq1 ${methy_fq1_clean} \
--fq2 ${methy_fq2_clean} \
--samplename ${sample} \
--outdir ${methy_dir} \
--barcode ${U3CB_methylation} \
--chemistry DD-M \
--core ${core}6.3.2 Step2：比对和去重
将step1得到的clean reads进行以下分析：
- 在read名称末尾添加UMI，方便后续使用deduplicate_bismark结合UMI信息进行去重
- 使用bismark比对，得到bam文件
- 得到的bam文件，根据read name进行排序，然后使用deduplicate_bismark进行去重
输入：${sample}_1.fq.gz、${sample}_2.fq.gz
输出：
${sample}_1_rename_bismark_bt2_pe.bam${sample}_1_rename_bismark_bt2_PE_report.txt${sample}_1_rename_bismark_bt2_pe_sortbyname.bam${sample}_1_rename_bismark_bt2_pe_sortbyname.deduplicated.bam${sample}_1_rename_bismark_bt2_pe_sortbyname.deduplication_report.txt
WARNING
重要注意事项：
${bismark_core}=${core}/8，如果bismark--parallel参数设置为8，确保机器的CPU数量至少为64- bismark比对时间较长，如果数据量大，可能会运行时间较长，建议使用配置较高的机器
- 例如，500G数据量，
--parallel设置为8，64个CPU机器大概运行30个小时
bismark具体参数解释，请使用--help参数或详见bismark官网。
mkdir -p ${methy_dir}/step2/bismark
# 向read name中添加UMI信息
python ${script_path}/add_umi_to_fastq_end_of_read_name.py \
--fq1 ${methy_dir}/step1/${sample}_1.fq.gz \
--fq2 ${methy_dir}/step1/${sample}_2.fq.gz \
--samplename ${sample} \
--outdir ${methy_dir}/step2/bismark
# bismark比对，实际使用的CPU应该是设置的core的4-5倍
bismark \
--parallel ${bismark_core} \
--genome ${bismark_genome} \
-1 ${methy_dir}/step2/bismark/${sample}_1_rename.fq.gz \
-2 ${methy_dir}/step2/bismark/${sample}_2_rename.fq.gz \
-o ${methy_dir}/step2/bismark/ \
--temp_dir ${methy_dir}/step2/bismark/ \
--non_directional
# deduplicate_bismark要求bam文件按照read名称排序，使用samtools或者sambamba进行排序
samtools sort \
-@ ${core} \
-n \
-o ${methy_dir}/step2/bismark/${sample}_1_rename_bismark_bt2_pe_sortbyname.bam \
${methy_dir}/step2/bismark/${sample}_1_rename_bismark_bt2_pe.bam;
# deduplicate_bismark 去重
deduplicate_bismark \
-p \
--output_dir ${methy_dir}/step2/bismark/ \
--barcode \
--bam ${methy_dir}/step2/bismark/${sample}_1_rename_bismark_bt2_pe_sortbyname.bam# 查看bismark比对的输出结果文件
cat ${methy_dir}/${sample}_1_rename_bismark_bt2_PE_report.txt计算比对率为：
唯一比对率为：
6.3.3 ALLCools统计样本中所有的甲基化位点和甲基化水平
ALLCools具体参数解释，请使用--help参数或详见ALLCools官网。
输入：${sample}_1_rename_bismark_bt2_pe_sortbyname.deduplicated.bam
输出：${sample}_allc.gz
# allcools要求bam文件按照染色体位置排序，需要使用samtools或者sambamba进行排序
samtools sort \
-@ ${core} \
-o ${methy_dir}/step2/bismark/${sample}_1_deduplicated_sort.bam \
${methy_dir}/step2/bismark/${sample}_1_rename_bismark_bt2_pe_sortbyname.deduplicated.bam
# allcools统计CNN位点的methylation reads和coverage reads
allcools bam-to-allc \
--bam_path ${methy_dir}/step2/bismark/${sample}_1_deduplicated_sort.bam \
--convert_bam_strandness \
--reference_fasta ${genomefa} \
--output_path ${methy_dir}/step2/bismark/${sample}_allc# 查看输出文件
# 第1列：chromosome
# 第2列：position(1-based)
# 第3列：strand
# 第4列：sequence context
# 第5列：count of reads supporting methylation
# 第6列：read coverage
# 第7列：indicator of significant methylation (1 if no test is performed)
less ${methy_dir}/step2/bismark/${sample}_allc.gzALLC文件格式示例：
6.3.4 使用featureCounts对bam进行注释
输入：${sample}_1_rename_bismark_bt2_pe_sortbyname.bam
输出：${sample}_1_rename_bismark_bt2_pe_sortbyname.bam.featureCounts.bam
CAUTION
注意：-t参数应当根据所选参考基因组而决定
mkdir -p ${methy_dir}/step2/wgs
cd ${methy_dir}/step2/wgs
featureCounts \
-T ${core} -t gene -s 0 -M -O -g gene_id --fracOverlap 0.5 \
-a ${gtf} \
-p \
--countReadPairs \
-o ${methy_dir}/step2/wgs/wgs_gene_counts.txt \
-R BAM \
${methy_dir}/step2/bismark/${sample}_1_rename_bismark_bt2_pe_sortbyname.bam
samtools sort \
-@ ${core} \
-o ${methy_dir}/step2/wgs/${sample}_sort.bam \
${methy_dir}/step2/bismark/${sample}_1_rename_bismark_bt2_pe_sortbyname.bam6.3.5 统计覆盖度
统计样本在基因组的覆盖度，并将比对信息和覆盖度信息写入json文件中。
输入：${sample}_1_rename_bismark_bt2_pe_sortbyname.bam
输出：bedtools_genomecoverage.txt
samtools sort \
-@ ${core} \
-o ${methy_dir}/step2/wgs/${sample}_sort.bam \
${methy_dir}/step2/bismark/${sample}_1_rename_bismark_bt2_pe_sortbyname.bam
bedtools genomecov \
-ibam ${methy_dir}/step2/wgs/${sample}_sort.bam > ${methy_dir}/step2/wgs/bedtools_genomecoverage.txt
python ${script_path}/wgs_coverage_depth.py \
--outdir ${methy_dir}/step2/wgs/ \
--samplename ${sample} \
--align_summary ${methy_dir}/step2/bismark/${sample}_1_rename_bismark_bt2_PE_report.txt \
--summary_json ${methy_dir}/${sample}_summary.json \
--coveragefile ${methy_dir}/step2/wgs/bedtools_genomecoverage.txt6.4 Step3：细胞判定和单细胞分析
6.4.1 细胞判定
统计每个barcode含有的基因数（若无featureCounts注释的基因，则以染色体+起始位置作为feature），umi数，reads数。
输入：${sample}_1_rename_bismark_bt2_pe_sortbyname.bam.featureCounts.bam
输出：counts.xls、detail.xls、umi.xls
samtools sort \
-n -O BAM -@ ${core} \
-o ${methy_dir}/step2/wgs/${sample}_featureCounts_SortByName.bam \
${methy_dir}/step2/wgs/${sample}_1_rename_bismark_bt2_pe_sortbyname.bam.featureCounts.bam
python ${script_path}/step3dnam3.py \
--bam ${methy_dir}/step2/wgs/${sample}_featureCounts_SortByName.bam \
--outdir ${methy_dir}/step3# 查看输出文件
# 第1列是barcode
# 第2列是feature
# 第3列是umi数量
# 第4列是reads数量
less ${methy_dir}/step3/counts.xlsCounts文件格式示例：
甲基化数据根据转录组判定的细胞barcode ${gexcb}来过滤得到细胞的umicounts
输入：counts.xls
输出：raw_umicounts.xls、filter_umicounts.xls
# --expected 为期望捕获细胞数
# raw_umicounts.xls 为所有barcode按umi数排序的结果
# filter_umicounts.xls 为判定为细胞的barcode按umi数排序的结果
python ${script_path}/wgs_umi_count_cs.py \
--infile ${methy_dir}/step3/counts.xls \
--rawcsv ${methy_dir}/step3/raw_umicounts.xls \
--filtercsv ${methy_dir}/step3/filter_umicounts.xls \
--gexcb ${gexcb} \
--cbcsv ${cbcsv} \
--outdir ${methy_dir} \
--samplename ${sample}6.4.2 拆分bam
将判断为细胞的barcode输出到filtered_barcode。将去重后的bam按照filtered_barcode进行拆分为单个细胞的bam文件。
输入：filtered_barcode、${sample}_1_rename_bismark_bt2_pe_sortbyname.deduplicated.bam
输出：split_bams目录
# 将判断为细胞的 barcode 输出到 filtered_barcode
cat ${methy_dir}/step3/filter_umicounts.xls |cut -f 1 |grep -v 'Barcode' > ${methy_dir}/step3/filtered_barcode
# 拆分
python ${script_path}/split_bams.py \
--bam ${methy_dir}/step2/bismark/${sample}_1_rename_bismark_bt2_pe_sortbyname.deduplicated.bam \
--outdir ${methy_dir}/step3 \
--samplename ${sample} \
--filtered_barcode ${methy_dir}/step3/filtered_barcode6.4.3 单细胞甲基化水平统计
对每个细胞的bam文件，使用samtools进行排序，使用ALLCools统计每个细胞每个CNN位点的甲基化水平。
输入：split_bams目录，filtered_barcode
输出：allcools目录，allcools_generate_datasets目录
python ${script_path}/step3_run_allcools_and_generate_datasets.py \
--indir ${methy_dir}/step3/split_bams \
--samplename ${sample} \
--outdir ${methy_dir}/step3/ \
--genomefa ${genomefa} \
--chrom_size_path ${chrom_size_path} \
--filtered_barcode ${methy_dir}/step3/filtered_barcode统计每个细胞的覆盖度
输入：allcools目录
输出：bed目录、cb_cov目录、${sample}_cb_genomecov.xls
python ${script_path}/wgs_sort_ml_anno_dev.py \
--filtercb ${methy_dir}/step3/filtered_barcode \
--core ${core} \
--outdir ${methy_dir}/step3/ \
--samplename ${sample} \
--genomefa ${genomefa} \
--allcpath ${methy_dir}/step3/allcools# 查看sample_cb_genomecov.xls文件
# 第1列为barcode
# 第2列为1X深度下的覆盖度
less ${methy_dir}/step3/${sample}_cb_genomecov.xls基因组覆盖度示例：
统计每个细胞的CpG数量
输入：allcools目录，bed目录，filtered_barcode文件，${sample}_allc.gz
输出：filter_gcov_umi_sort.xls
python ${script_path}/wgs_mk_martix_bins.py \
--cbfile ${methy_dir}/step3/filtered_barcode \
--outdir ${methy_dir}/step3 \
--samplename ${sample} \
--outcsv ${methy_dir}/step3/${sample}_CPGnum_cb.xls \
--allcpgfile ${methy_dir}/step2/bismark/${sample}_allc.gz \
--indir ${methy_dir}/step3/bed \
--summary_json ${methy_dir}/${sample}_summary.json生成最终的质控表
输入：summary.json
输出：${sample}_wgs_summary.csv
python ${script_path}/wgs_summary_yf.py \
--outdir ${methy_dir} \
--samplename ${sample} \
--summary_json ${methy_dir}/${sample}_summary.json6.5 Step4：单细胞甲基化数据降维聚类
利用ALLCools软件进行单细胞甲基化数据降维聚类
输入：allcools_generate_datasets目录
输出：step4目录，包含tsne、umap降维聚类图，h5ad文件。
mkdir -p ${methy_dir}/step4
python ${script_path}/step4_allcools_PCA_cluster.py \
--mcds_path ${methy_dir}/step3/allcools_generate_datasets/${sample}.mcds \
--samplename ${sample} \
--var_dim chrom1M \
--filtered_barcode_file ${methy_dir}/step3/filtered_barcode \
--outdir ${methy_dir}/step4单细胞甲基化数据聚类结果示例：
7. 注意事项与最佳实践
7.1 数据质量要求
IMPORTANT
为确保分析结果的准确性，请确保：
- 测序数据质量良好，Q30比例>80%
- 文库构建质量符合预期，CT转化率在预期范围内
- 细胞数量和质量符合实验设计
7.2 计算资源优化
TIP
为提高分析效率，建议：
- 使用SSD存储加速I/O操作
- 合理设置并行参数，避免过度占用系统资源
- 定期监控系统资源使用情况
7.3 结果解读
NOTE
分析完成后，请重点关注：
- 细胞捕获效率和质控指标
- 甲基化数据的覆盖度和深度
- 转录组和甲基化数据的一致性
- 聚类结果的生物学意义
8. 常见问题与解决方案
8.1 数据下载问题
如果数据下载失败，请检查：
- 网络连接是否正常
- 磁盘空间是否充足
- 下载链接是否过期
8.2 分析流程问题
如果分析过程中出现错误，请检查：
- 输入文件路径是否正确
- 系统资源是否充足
- 软件版本是否兼容
8.3 结果质量问题
如果结果质量不理想，请检查：
- 原始数据质量
- 参数设置是否合理
- 质控步骤是否完整
9. 总结
本教程详细介绍了寻因单细胞转录组+甲基化双组学分析的完整流程，包括软件安装、数据准备、分析执行和结果解读等各个环节。通过本教程的学习，用户应该能够独立完成单细胞转录组+甲基化双组学数据的分析工作。
TIP
建议在实际分析前，先使用本教程提供的demo数据进行练习，熟悉整个分析流程，然后再进行实际数据的分析。
随着单细胞多组学技术的不断发展，我们期待更多的分析方法和工具能够为科学研究提供更强大的支持。