SCOPE是一个统计框架，用于从全基因组单细胞DNA测序读取深度调用拷贝数变异(CNVs)。SCOPE的显著特征包括:
利用细胞特异性基尼系数进行质量控制和鉴定正常/二倍体细胞。在大多数单细胞癌症基因组学研究中，二倍体细胞不可避免地从邻近的正常组织中提取用于测序，因此可以作为正常对照进行读深归一化。然而，并非所有的平台/实验都在scDNA-seq之前允许或采用基于流排序的技术，因此细胞倍性和病例对照标记并不总是容易获得。基尼系数能够从整个细胞群中索引二倍体细胞，并作为识别细胞异常值的良好代理。
采用EM算法来模拟GC含量偏差，这说明沿基因组不同的拷贝数状态。SCOPE基于泊松潜在因子模型进行跨样本归一化，借用跨区域和跨样本的信息来估计偏差项。
结合多样本分割程序，以识别来自相同遗传背景的细胞间共享的断点
我们在这里演示预处理生物信息管道是如何工作的。的皮卡德工具包是开源和免费的所有用途。的split_script.pyPython脚本被预先存储在包中，用于解复用。
有两种类型的scDNA-seq数据源:来自NCBI Sequence Read Archive的公共数据和来自10X Genomics的数据。对于NCBI SRA数据，请从SRA文件开始。FASTQ -dump获取FASTQ文件。将FASTQ序列与NCBI hg19参考基因组对齐并转换为BAM文件。对于10X基因组数据集，从原始的集成BAM文件进行处理。每次读取的错误校正铬细胞条形码信息存储为CB标签字段。只有包含CB标签并且在感兴趣的条形码列表中的读取才会通过Python脚本进行多路解复用。在对齐/解复用的bam上排序、添加读组和重复数据删除。使用重删的BAM文件作为输入。
#公共数据从NCBI序列读取存档SRR =SRRXXXXXXX金=/松/ scr / r / u /如金/ Kim_Navin_et_al_Cell_2018金fastq_dir = $/ fastq金align_dir = $/对齐将FASTQ序列与NCBI hg19参考基因组比对单端测序单元格只有1个FASTQ文件;#配对端测序将生成两个FASTQ文件，#带后缀“_1”和“_2”)cdfastq_dir美元bwamem -M -t 16 \｀ls|grep＂SRR美元＂|tr' \ n '＇ ＇｀>align_dir美元/＂SRR美元＂.sam#将。sam转换为。bamcdalign_dir美元samtools视图废话＂SRR美元＂.sam>＂SRR美元＂本演讲#排序java-Xmx30G -jar /proj/yuchaojlab/bin/picard.jar SortSam \＂SRR美元＂本输出=＂SRR美元＂.sorted。bam \#添加读组java-Xmx40G -jar /proj/yuchaojlab/bin/picard.jar AddOrReplaceReadGroups \＂SRR美元＂.sorted。bam O =＂SRR美元＂.sorted.rg.bamRGID=＂SRR美元＂＼＂SRR美元＂samtools指数＂SRR美元＂.sorted.rg.bam# Dedupjava-Xmx40G -jar /proj/yuchaojlab/bin/picard.jar markduplicate \＂SRR美元＂.sorted.rg.bamO=＂SRR美元＂.sorted.rg.dedup.bam \＂SRR美元＂.sorted.rg.dedup.metrics.txt \java-jar /proj/yuchaojlab/bin/picard.jar BuildBamIndex \＂SRR美元＂.sorted.rg.dedup.bam# 10X基因组学XGenomics =/松/ scr / r / u /如金/ 10 xgenomics数据集=breast_tissue_A_2koutput_dir = $ XGenomics/元数据集/输出align_dir = $ XGenomics/元数据集/对齐#分工cdoutput_dir美元samtools视图${数据集}_possorted_bam.bam|pythonXGenomics美元/ split_script.py#为解复用的bam文件添加头以便进一步处理cdXGenomics美元samtools视图- h元数据集/输出/${数据集}_possorted_bam.bam>＼元数据集/ header.txt条形码=AAAGATGGTGTAAAGT猫header.txtalign_dir美元/美元的条形码/美元的条形码-1.山姆>＼align_dir美元/美元的条形码/美元的条形码-1. header.sam#将。sam转换为。bamcdalign_dir美元/美元的条形码samtools视图废话＂美元的条形码＂-1. header.sam>＂美元的条形码＂本演讲SCOPE允许在下游分析之前通过指定参数重建用户定义的具有固定间隔长度的全基因组连续bin基因组而且决议在get_bam_bed ()函数。确保所有染色体的命名一致并与本演讲文件。SCOPE处理整个基因组。使用功能get_bam_bed ()完成预准备步骤。缺省情况下，SCOPE为hg19设计，bin-length = 500kb固定。
计算每个bin的GC内容和可映射性。默认情况下，SCOPE用于hg19参考基因组。为了计算hg19的可映射性，我们使用了ENCODE项目的100-mers可映射性跟踪(wgEncodeCrgMapabilityAlign100mer.bigwig从链接)，如果多个ENCODE区域与同一bin重叠，则计算可映射性得分的加权平均。对于SCOPE，人类hg19装配的全基因组可映射性跟踪作为包的一部分存储。
人类hg38装配的全基因组可映射性跟踪也存储在SCOPE包中。有关可映射性计算的详细信息，请参见CODEX2 for hg38．加载hg38参考包，指定参数Hgref = "hg38"在get_mapp ()而且get_gc ()．默认情况下,hg19使用。
## [1] 0.1672858 0.5042187 0.9529839 0.9065575 0.9826570 0.9100177get_gc(ref_rawhgref =“hg38”）值(ref_raw) < -cbind（值(ref_raw),DataFrame攻读硕士学位(gc))小鼠基因组的CODEX2)．指定参数Hgref = "mm10"在get_bam_bed ()，get_mapp ()，get_gc ()而且get_coverage_scDNA ()．GC内容和可映射性的计算需要从默认值(hg19)修改。对于mm10，有两种解决方案:1)将所有的可映射性设置为1，以避免大量的计算;2)采用基于标注结果的QC流程，如过滤掉“黑名单”区域内的箱子，这些箱子的可映射性一般较低。具体来说，我们获得并预存储了小鼠基因组“黑名单”区域Amemiya等人，科学报告，2019年．
对于未预先计算可映射性轨迹的未知参考程序集，请参阅CODEX2:可映射性预计算．
获得单端或配对端测序读深度矩阵。在默认情况下，SCOPE采用固定的binning方法来计算覆盖深度，同时删除映射到多个基因组位置和“黑名单”区域的读取。接下来是一个附加的质量控制步骤，以删除具有极端可映射性的箱子，以避免错误检测。具体来说，“黑名单”垃圾箱，包括节段复制区而且参考组件中的间隙来自端粒、着丝粒和/或异染色质区域。
示例1的覆盖率:AAAGCAATCTGACGCG…示例2:CTCGTCACAGGTTAAA…示例3:GCAGTTACACTGTATG…get_samp_QC ()用于对单个cell执行QC步骤，返回总读数/占比、总读数/映射读数占比、总读数/映射非重复读数占比、映射质量大于20的读数/占比。使用perform_qc ()进一步去除低比例映射读取的样本/细胞，具有极端GC含量(小于20%和大于80%)和低可映射性(小于0.9)的容器，以减少伪影。
##获取样品1的样品质量控制度量##获取样品2的样品质量控制度量##获取样品3的样品质量控制度量由于库准备失败，删除了0个样本。##由于未能满足最小覆盖要求，删除了0个样本。由于映射读取的比例低，删除了0个样本。由于垃圾收集内容过多，排除了216个垃圾箱。由于可映射性低，排除了355个箱子。删除了0个样本，因为在##库大小计算中有过多的零读计数。QC步骤后有3个样品和5029个箱子。SCOPE的一个特点是使用基尼指数来识别正常/二倍体细胞。每个单元格的基尼系数是洛伦兹曲线与对角线之间面积的2倍。基尼指数的值在0到1之间变化，其中0是最均匀的，1是最极端的。建议排除基尼系数极高(大于0.5)的细胞。建立一个识别二倍体/正常细胞的Gini阈值(例如Gini小于0.12)。我们演示预先存储的玩具数据集如下。
正常细胞指数由基尼系数或先验知识决定。SCOPE利用来自第一次通过归一化运行的倍性估计来确保快速收敛并避免局部最优。指定SoS。情节=真实在initialize_ploidy ()使倍性预估可视化。归一化过程包括泊松广义线性模型中的期望最大化算法。注意，通过设置ploidyInt在归一化函数中，SCOPE允许用户利用先验知识倍体作为输入，并手动调优一些拟合不佳的单元格。
#第一次通过CODEX2运行时没有潜在因素normalize_codex2_ns_noK（Y_qc =Y_sim,gc_qc =ref_sim＄gc,norm_index =哪一个(基尼< =0.12）)没有潜在因素的计算归一化##迭代1 beta diff =0.00882 f(GC) diff =4.28e-07##迭代2 beta diff =2.21e-06 f(GC) diff =9.35e-12在迭代2中停止。初始化单元格1的倍性#如果使用高性能集群，并行计算是#简单，提高计算效率。简单地使用# normalize_scope_foreach()代替normalize_scope()。#所有参数相同。normalize_scope_foreach（Y_qc =Y_sim,gc_qc =ref_sim＄gc,K =1，ploidyInt =ploidy.sim,norm_index =哪一个(基尼< =0.12)，T =1：5，beta0 =normObj.sim＄beta.hat,nCores =2）##初始化…## 1 2 3 4 5用k = 1个潜在因子计算归一化…## k = 1##迭代1 beta diff =3.43e-05 f(GC) diff =4.05e-07## hhat diff =0.37## hhat diff =9.83e-06在迭代1中停止。## aic1 = 50521100.706## bic1 = 50510332.725## rss1 = 2441.606# normalize_scope(Y_qc = Y_sim, gc_qc = ref_sim$gc，# K = 1, ploidyInt = ploidy.sim，# norm_index = which(基尼系数<=0.12)，T = 1:5，# beta0 = beta.hat.noK.sim)normObj.scope.sim＄Yhat [[which.max(normObj.scope.sim＄BIC)]]normObj.scope.sim＄fGC.hat [[which.max(normObj.scope.sim＄BIC)]]可视化每个单元格的选择结果。缺省情况下，使用BIC选择最优的CNV组。
plot_EM_fit（Y_qc =Y_sim,gc_qc =ref_sim＄gc,norm_index =哪一个(基尼< =0.12)，T =1：5，ploidyInt =ploidy.sim,beta0 =normObj.sim＄beta.hat,文件名=“plot_EM_fit_demo.pdf”）在完成了SCOPE的默认归一化和分割之后，SCOPE包含了基于归一化z分数矩阵、估计拷贝数或估计变更点的聚类单元格的选项。对于推断出的子克隆，SCOPE可以选择执行第二轮组级倍性初始化和规范化
#组型倍性初始化c（“正常”，“tumor1”，“正常”，“tumor1”，“tumor1”）initialize_ploidy_group（Y =Y_sim,Yhat =Yhat.noK.sim,ref =ref_sim,组=克隆)#分组规范化normalize_scope_group（Y_qc =Y_sim,gc_qc =ref_sim＄gc,K =1，ploidyInt =ploidy.sim.group,norm_index =哪一个(克隆= =“正常”)，组=克隆,T =1：5，beta0 =beta.hat.noK.sim)normObj.scope.sim.group＄Yhat [[which.max（＄BIC)]]normObj.scope.sim.group＄fGC.hat [[which.max（＄BIC)]]SCOPE提供了跨样本分割，它在来自同一克隆的单元格之间输出共享断点。这一步对整个基因组染色体逐一进行处理。将返回共享断点和整数副本编号配置文件。
SCOPE提供了推断的整数拷贝数配置文件的热图，其中细胞通过分层聚类进行聚类。
## R版本4.2.1(2022-06-23)##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 20.04.5 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack。所以## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_GB LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4 stats graphics grDevices utils datasets methods ##[8]基础## ##其他附加包:[1] BSgenome.Hsapiens.UCSC.hg38_1.4.4 WGSmapp_1.9.0 ## [5] BSgenome_1.66.0 rtracklayer_1.58.0 ## [7] Rsamtools_2.14.0 Biostrings_2.66.0 ## [9] XVector_0.38.0 GenomicRanges_1.50.0 ## [13] S4Vectors_0.36.0 BiocGenerics_0.44.0 ## ##通过命名空间加载(并且没有附加):[1] MatrixGenerics_1.10.0 mvtnorm_1.1-3 ## [25] proxy_0.4-27 cachem_1.0.6 ## [9] expm_0.999-6 gld_2.6.6 ## [11] lmom_2.9 GenomeInfoDbData_1.2.9 ## [13] cellranger_1.1.0 yaml_2.3.6 ## [15] lattice_0.20-45 digest_0.6.30 ## [17] RColorBrewer_1.1-3 htmltools_0.5.3 ## [19] Matrix_1.5-1 XML_3.99-0.12 ## [23] rootSolve_1.8.2.3 [25] proxy_0.4-27 cachem_1. 32.0 ## # [27][39] BiocIO_1.8.0 matrixStats_0.62.0 ## [41] tools_4.2.1 data.table_1.14.4 ## [43] DelayedArray_0.24.0 compiler_4.2.1 ## [45] jquerylib_0.1.4 e1071_1.7-12 ## [47] caTools_1.18.2 rlang_1.0.6 ## [49] grid_4.2.1 RCurl_1.98-1.9 ## b[51] iterators_1.0.14 rstudioapi_0.14 ## [39]rjson_0.2.21 bitops_1.0-7 ## [55] rmarkdown_2.17 boot_1.3-28 ## [57] DNAcopy_1.72.0 DescTools_0.99.47 ## [59] restfulr_0.0.15 codetools_0.2-18 ## [61] R6_2.5.1 GenomicAlignments_1.34.0 ## [63] knitr_1.40 fastmap_1.1.0 ## [65] KernSmooth_2.23-20 stringi_1.7.8 ## [67] parallel_4.2.1 Rcpp_1.0.9 ## [69] xfun_0.34