科研-生信学习

ArchR官网教程学习笔记3:创建ArchRProject

2020-11-20  本文已影响0人  生信start_site

系列回顾：
ArchR官网教程学习笔记1：Getting Started with ArchR
ArchR官网教程学习笔记2:基于ArchR推测Doublet

ArchRProject允许我们将多个 Arrow files组合到一个项目中。这个ArchRProject很小，并且存储在内存中。通过与ArchRProject交互，我们可以快速地从 Arrow files中抓取数据。因此，它构成了几乎所有ArchR函数和分析工作流程的基础。此外，ArchRProject对象可以在保存后重新加载，提供了分析的连续性，并方便合作者共享分析项目。本章描述了如何创建ArchRProject对象并与之交互。

在开始前，要先注意，不是所有的官网教程的代码在你电脑上都没问题，比如这一章的代码，由于使用的系统不一样，代码的效果也会很不一样，甚至有些会报错！这就是为什么要自己跟着官网走一遍，这样才知道哪一步是你要注意的地方！

（一）创建ArchRProject

首先，我们必须通过提供一个Arrow files列表和一些其他参数来创建ArchRProject。这里的outputDirectory表示所有下游分析和绘图的保存文件夹路径。ArchR会自动将之前提供的geneAnnotation和genomeAnnotation与新的ArchRProject关联起来。

> projHeme1 <- ArchRProject(
  ArrowFiles = ArrowFiles, 
  outputDirectory = "HemeTutorial",
  copyArrows = TRUE #This is recommened so that if you modify the Arrow files you have an original copy for later usage.
)
#运行时的输出：
Using GeneAnnotation set by addArchRGenome(Hg19)!
Using GeneAnnotation set by addArchRGenome(Hg19)!
Validating Arrows...
Getting SampleNames...
1 2 3 
Copying ArrowFiles to Ouptut Directory! If you want to save disk space set copyArrows = FALSE
1 2 3 
Getting Cell Metadata...
1 2 3 
Merging Cell Metadata...
Initializing ArchRProject...

我们把这个ArchRProject为“projHeme1”。在本演练中，我们将修改和更新这个ArchRProject，并通过迭代项目号(即“projHeme2”)来跟踪我们使用的版本。

我们可以检查我们的ArchRProject的内容:

> projHeme1
    
class: ArchRProject 
outputDirectory: D:\yanfang\HemeTutorial 
samples(3): scATAC_BMMC_R1 scATAC_CD34_BMMC_R1 scATAC_PBMC_R1
sampleColData names(1): ArrowFiles
cellColData names(15): Sample TSSEnrichment ... DoubletEnrichment BlacklistRatio
numberOfCells(1): 10661
medianTSS(1): 16.832
medianFrags(1): 3050

从上面我们可以看到，ArchRProject已经初始化了几个重要的属性:

1.指定的输出文件夹。
2.从 Arrow files中获得的每个样品的样品名。
3.一个名为sampleColData的矩阵，它包含与每个样本相关的数据。
4.一个名为cellColData的矩阵，它包含与每个细胞关联的数据。因为我们已经使用addDoubletScores()计算了doublet富集分数，它将这些值添加到Arrow files中的每个细胞中，所以我们可以看到cellColData矩阵中与“DoubletEnrichment”和“DoubletScore”对应的列。
5.我们project中的细胞总数，代表了doublets鉴定和去除后的所有样本。
6.TSS富集分数的中位数和所有细胞和所有样本的fragments数量的中位数。

我们还可以检查储存ArchRProject用了多少内存：

> paste0("Memory Size = ", round(object.size(projHeme1) / 10^6, 3), " MB")
[1] "Memory Size = 37.477 MB"

我们也可以询问ArchRProject哪一个矩阵是可以用来进行下游分析的：

> getAvailableMatrices(projHeme1)
[1] "GeneScoreMatrix" "TileMatrix"

（二）对ArchRProject进行操作

现在我们已经创建了一个ArchRProject对象，现在我们可以做很多事情来访问数据，并对相关数据进行操作。

（1）例一：$符号允许直接访问cellColData

我们可以访问细胞名称：

> head(projHeme1$cellNames)
[1] "scATAC_BMMC_R1#TTATGTCAGTGATTAG-1" "scATAC_BMMC_R1#AAGATAGTCACCGCGA-1"
[3] "scATAC_BMMC_R1#GCATTGAAGATTCCGT-1" "scATAC_BMMC_R1#TATGTTCAGGGTTCCC-1"
[5] "scATAC_BMMC_R1#TCCATCGGTCCCGTGA-1" "scATAC_BMMC_R1#AGTTACGAGAACGTCG-1"

我们可以访问每个细胞相对的样品名：

> head(projHeme1$Sample)
[1] "scATAC_BMMC_R1" "scATAC_BMMC_R1" "scATAC_BMMC_R1" "scATAC_BMMC_R1" "scATAC_BMMC_R1"
[6] "scATAC_BMMC_R1"

我们还可以访问每个细胞的TSS富集分数：

> quantile(projHeme1$TSSEnrichment)
    0%    25%    50%    75%   100% 
 4.027 13.922 16.832 19.937 41.782
（2）例二：按照细胞来取ArchRProject子集

有多种方法：

我们可以根据行的数量来取子集：

> projHeme1[1:100, ]
    
class: ArchRProject 
outputDirectory: D:\yanfang\HemeTutorial 
samples(3): scATAC_BMMC_R1 scATAC_CD34_BMMC_R1 scATAC_PBMC_R1
sampleColData names(1): ArrowFiles
cellColData names(15): Sample TSSEnrichment ... DoubletEnrichment BlacklistRatio
numberOfCells(1): 100
medianTSS(1): 10.7725
medianFrags(1): 10200.5

还可以根据细胞名来取子集：

> projHeme1[projHeme1$cellNames[1:100], ]
    
class: ArchRProject 
outputDirectory: D:\yanfang\HemeTutorial 
samples(3): scATAC_BMMC_R1 scATAC_CD34_BMMC_R1 scATAC_PBMC_R1
sampleColData names(1): ArrowFiles
cellColData names(15): Sample TSSEnrichment ... DoubletEnrichment BlacklistRatio
numberOfCells(1): 100
medianTSS(1): 10.7725
medianFrags(1): 10200.5

我们可以所有的细胞都对应于一个特定的样本:

> idxSample <- BiocGenerics::which(projHeme1$Sample %in% "scATAC_BMMC_R1")
> cellsSample <- projHeme1$cellNames[idxSample]
> projHeme1[cellsSample, ]

class: ArchRProject 
outputDirectory: D:\yanfang\HemeTutorial 
samples(3): scATAC_BMMC_R1 scATAC_CD34_BMMC_R1 scATAC_PBMC_R1
sampleColData names(1): ArrowFiles
cellColData names(15): Sample TSSEnrichment ... DoubletEnrichment BlacklistRatio
numberOfCells(1): 4932
medianTSS(1): 15.254
medianFrags(1): 2771

还可以根据TSS富集分数来取子集：

> idxPass <- which(projHeme1$TSSEnrichment >= 8)
> cellsPass <- projHeme1$cellNames[idxPass]
> projHeme1[cellsPass, ]
    
class: ArchRProject 
outputDirectory: D:\yanfang\HemeTutorial 
samples(3): scATAC_BMMC_R1 scATAC_CD34_BMMC_R1 scATAC_PBMC_R1
sampleColData names(1): ArrowFiles
cellColData names(15): Sample TSSEnrichment ... DoubletEnrichment BlacklistRatio
numberOfCells(1): 10500
medianTSS(1): 16.9275
medianFrags(1): 3042
（3）例三：向ArchRProject里添加数据

我们可以在cellColData里添加列，来储存任何与我们project相关的细胞特异的metadata。比如，在cellColData添加列，通过删除原始样本名称中多余的信息使样本名称更易读:

> bioNames <- gsub("_R2|_R1|scATAC_","",projHeme1$Sample)
> head(bioNames)
[1] "BMMC" "BMMC" "BMMC" "BMMC" "BMMC" "BMMC"
> projHeme1$bioNames <- bioNames

或者，我们可以使用addCellColData()往cellColData添加列，ArchR允许只针对某些细胞添加多于的列。

> bioNames <- bioNames[1:10]
> cellNames <- projHeme1$cellNames[1:10]
> projHeme1 <- addCellColData(ArchRProj = projHeme1, data = paste0(bioNames),
    cells = cellNames, name = "bioNames2")

对于剩下的没有填充的数值，Arch会用NA来填充：

> getCellColData(projHeme1, select = c("bioNames", "bioNames2"))
DataFrame with 10661 rows and 2 columns
                                     bioNames   bioNames2
                                  <character> <character>
scATAC_BMMC_R1#TTATGTCAGTGATTAG-1        BMMC        BMMC
scATAC_BMMC_R1#AAGATAGTCACCGCGA-1        BMMC        BMMC
scATAC_BMMC_R1#GCATTGAAGATTCCGT-1        BMMC        BMMC
scATAC_BMMC_R1#TATGTTCAGGGTTCCC-1        BMMC        BMMC
scATAC_BMMC_R1#TCCATCGGTCCCGTGA-1        BMMC        BMMC
...                                       ...         ...
scATAC_PBMC_R1#GCTGCGAAGATCCGAG-1        PBMC          NA
scATAC_PBMC_R1#GCAGCTGGTGGCCTTG-1        PBMC          NA
scATAC_PBMC_R1#GCAGATTGTACGCAAG-1        PBMC          NA
scATAC_PBMC_R1#TTCGTTACATTGAACC-1        PBMC          NA
scATAC_PBMC_R1#CGCTATCGTGAGGTCA-1        PBMC          NA
（4）例四：获取cellColData里的列

ArchR提供了getCellColData()来提取metadata列。
比如：我们要根据列的名字来提取列，比如每个细胞里唯一nuclear的fragments数量：

> df <- getCellColData(projHeme1, select = "nFrags")
> df
DataFrame with 10661 rows and 1 column
                                     nFrags
                                  <numeric>
scATAC_BMMC_R1#TTATGTCAGTGATTAG-1     26189
scATAC_BMMC_R1#AAGATAGTCACCGCGA-1     20648
scATAC_BMMC_R1#GCATTGAAGATTCCGT-1     18991
scATAC_BMMC_R1#TATGTTCAGGGTTCCC-1     18296
scATAC_BMMC_R1#TCCATCGGTCCCGTGA-1     17458
...                                     ...
scATAC_PBMC_R1#GCTGCGAAGATCCGAG-1      1038
scATAC_PBMC_R1#GCAGCTGGTGGCCTTG-1      1037
scATAC_PBMC_R1#GCAGATTGTACGCAAG-1      1033
scATAC_PBMC_R1#TTCGTTACATTGAACC-1      1033
scATAC_PBMC_R1#CGCTATCGTGAGGTCA-1      1002

你还可以对给定的列名进行操作，输出结果：

> df <- getCellColData(projHeme1, select = c("log10(nFrags)", "nFrags - 1"))
> df
DataFrame with 10661 rows and 2 columns
                                  log10(nFrags) nFrags - 1
                                      <numeric>  <numeric>
scATAC_BMMC_R1#TTATGTCAGTGATTAG-1       4.41812      26188
scATAC_BMMC_R1#AAGATAGTCACCGCGA-1       4.31488      20647
scATAC_BMMC_R1#GCATTGAAGATTCCGT-1       4.27855      18990
scATAC_BMMC_R1#TATGTTCAGGGTTCCC-1       4.26236      18295
scATAC_BMMC_R1#TCCATCGGTCCCGTGA-1       4.24199      17457
...                                         ...        ...
scATAC_PBMC_R1#GCTGCGAAGATCCGAG-1       3.01620       1037
scATAC_PBMC_R1#GCAGCTGGTGGCCTTG-1       3.01578       1036
scATAC_PBMC_R1#GCAGATTGTACGCAAG-1       3.01410       1032
scATAC_PBMC_R1#TTCGTTACATTGAACC-1       3.01410       1032
scATAC_PBMC_R1#CGCTATCGTGAGGTCA-1       3.00087       1001
（5）绘制QC矩阵-- log10(Unique Fragments) vs TSS富集分数

重复上面的例子，我们可以很容易地获得用于单个细胞质量控制的标准scATAC-seq矩阵。我们发现，质量控制最可靠的度量标准是TSS富集分数(ATAC-seq数据中signal-to-background的度量)和唯一核fragments的数量。

> df <- getCellColData(projHeme1, select = c("log10(nFrags)", "TSSEnrichment"))
> df
DataFrame with 10661 rows and 2 columns
                                  log10(nFrags) TSSEnrichment
                                      <numeric>     <numeric>
scATAC_BMMC_R1#TTATGTCAGTGATTAG-1       4.41812         7.149
scATAC_BMMC_R1#AAGATAGTCACCGCGA-1       4.31488         7.911
scATAC_BMMC_R1#GCATTGAAGATTCCGT-1       4.27855         4.505
scATAC_BMMC_R1#TATGTTCAGGGTTCCC-1       4.26236         6.946
scATAC_BMMC_R1#TCCATCGGTCCCGTGA-1       4.24199         4.799
...                                         ...           ...
scATAC_PBMC_R1#GCTGCGAAGATCCGAG-1       3.01620        24.356
scATAC_PBMC_R1#GCAGCTGGTGGCCTTG-1       3.01578        22.537
scATAC_PBMC_R1#GCAGATTGTACGCAAG-1       3.01410        20.146
scATAC_PBMC_R1#TTCGTTACATTGAACC-1       3.01410        30.198
scATAC_PBMC_R1#CGCTATCGTGAGGTCA-1       3.00087        21.485

现在让我们通过TSS富集分数来绘制唯一核fragments的数量(log10)。这种类型的图是鉴别高质量细胞的关键。你要注意到，我们之前(通过filterTSS和filterFrags)创建Arrow files时指定的截断部分已经删除了低质量的细胞。但是，如果我们注意到之前应用的QC过滤不适用于这个样本，我们可以根据这个图进一步调整我们的cutoff，或者根据需要重新生成Arrow files。

> p <- ggPoint(
  x = df[,1], 
  y = df[,2], 
  colorDensity = TRUE,
  continuousSet = "sambaNight",
  xlabel = "Log10 Unique Fragments",
  ylabel = "TSS Enrichment",
  xlim = c(log10(500), quantile(df[,1], probs = 0.99)),
  ylim = c(0, quantile(df[,2], probs = 0.99))
) + geom_hline(yintercept = 4, lty = "dashed") + geom_vline(xintercept = 3, lty = "dashed")

> p

保存图片为pdf：

> plotPDF(p, name = "TSS-vs-Frags.pdf", ArchRProj = projHeme1, addDOC = FALSE)
Plotting Ggplot!

（三）绘制来自ArchRProject的统计信息

在处理单个集成数据集中的多个不同样本时，比较所有样本中的各种指标很重要。ArchR为分组数据提供了两种主要的绘图机制:山脊图和小提琴图。它们都可以通过plotGroups()函数访问。当然，这种绘图类型并不局限于样本数据，还可以用于绘制group等下游信息，比如clusters。

（1）为每个样本绘制TSS富集分数的山脊图
> p1 <- plotGroups(
  ArchRProj = projHeme1, 
  groupBy = "Sample", 
  colorBy = "cellColData", 
  name = "TSSEnrichment",
  plotAs = "ridges"
)
> p1
（2）为每个样品绘制TSS富集分数小提琴图
> p2 <- plotGroups(
  ArchRProj = projHeme1, 
  groupBy = "Sample", 
  colorBy = "cellColData", 
  name = "TSSEnrichment",
  plotAs = "violin",
  alpha = 0.4,
  addBoxPlot = TRUE
)
> p2
（3）为每个样品绘制log10（唯一核fragment）山峦图
> p3 <- plotGroups(
  ArchRProj = projHeme1, 
  groupBy = "Sample", 
  colorBy = "cellColData", 
  name = "log10(nFrags)",
  plotAs = "ridges"
)
> p3
（4）为每个样品绘制log10（唯一核fragment）小提琴图
> p4 <- plotGroups(
    ArchRProj = projHeme1, 
    groupBy = "Sample", 
    colorBy = "cellColData", 
    name = "log10(nFrags)",
    plotAs = "violin",
    alpha = 0.4,
    addBoxPlot = TRUE
   )
> p4

保存图片：

> plotPDF(p1,p2,p3,p4, name = "QC-Sample-Statistics.pdf", ArchRProj = projHeme1, addDOC = FALSE, width = 4, height = 4)

（四）绘制样品的Fragments大小分布图以及TSS富集情况

由于数据的存储和访问方式，ArchR可以非常快速地从Arrow files中计算fragment 大小分布和TSS富集。

绘制所有样本的fragments大小分布，我们使用plotfragmentsize()函数。在ATAC-seq中，fragments大小分布在不同的样本、细胞类型和批次之间变化很大。如下所示的细微差异是常见的，并不一定与数据质量的差异相关:

> p1 <- plotFragmentSizes(ArchRProj = projHeme1)
> p1

绘制TSS富集文件，我们使用plotTSSEnrichment()函数。TSS富集应在中心显示清晰的峰，中心右侧较小的峰，这是由于+1核小体的定位造成的。

> p2 <- plotTSSEnrichment(ArchRProj = projHeme1)
> p2

注意！注意！！在这一步骤中，我遇到了一个很奇怪的问题：我的电脑生成的图只画出了一个峰，另外两个样品没有峰。我google了一下，发现有这个问题的人还不在少数，大部分人在这一步有问题都是用的windows系统，参考：https://github.com/GreenleafLab/ArchR/issues/252这里的解决方法，其中一个开发者建议安装最新版本的ArchR，是一个月前左右刚release的（然而我试过以后并不work）。但是另一个开发者说“Windows is not supported and its not possible for us to troubleshoot this. ”（这话说的，windows难道就不配用这个R包了吗！）。

之后就不得已的又重启电脑切换成了我的Linux系统，在Rstudio里把前面的过程重新过了一遍，然后这一步果然没有问题！（看来真的是系统的问题。。。）见下图：

保存图片：

> plotPDF(p1,p2, name = "QC-Sample-FragSizes-TSSProfile.pdf", ArchRProj = projHeme1, addDOC = FALSE, width = 5, height = 5)

（五）保存并加载ArchRProject对象

ArchR提供了一种简便的方法来保存ArchRProject对象，以便稍后重新加载或与其他用户共享。基本上，ArchRProject指向一组Arrow files文件。因此，使用saveArchRProject()函数保存ArchRProject的过程将:

1.将当前Arrow files复制到指定的输出文件夹，以便它们与新的ArchRProject对象关联。
2.在outputDirectory中保存指定的ArchRProject的副本。

例如，我们可以使用saveArchRProject()来保存projHeme1，我们可以在以后的章节中使用这个对象。

> saveArchRProject(ArchRProj = projHeme1, outputDirectory = "Save-ProjHeme1", load = FALSE)
Copying ArchRProject to new outputDirectory : D:\yanfang\Save-ProjHeme1
Copying Arrow Files...
Copying Arrow Files (1 of 3)
Copying Arrow Files (2 of 3)
Copying Arrow Files (3 of 3)
Getting ImputeWeights
No imputeWeights found, returning NULL
Copying Other Files...
Copying Other Files (1 of 1): Plots
Saving ArchRProject...

这将复制Arrow files文件，并将projHeme1 ArchRProject对象的. rds文件保存在指定的outputDirectory中。非常重要! 此进程不会在当前R会话中自动更新正在使用的ArchRProject对象。具体来说，当前R会话中名为projHeme1的对象仍然指向Arrow files的原始位置，而不是位于指定outputDirectory中的复制的箭头文件。如果我们想这样做，我们要指定load = TRUE，这将导致saveArchRProject()函数返回保存的ArchRProject对象，你可以分配该对象来覆盖原始的ArchRProject对象。这样就可以有效地保存并从它的新位置加载ArchRProject。

（六）从ArchRProject里过滤doublets

在使用addDoubletScores()添加了关于预测的doublets的信息之后，可以使用filterDoublets()删除这些doublets。这一步最关键的要素就是filterRatio，filterRatio允许你在多个不同的样本上应用一致的filter，这些样本可能有不同的doublets百分比。filterRatio越高，被过滤掉的可能的doublets数量就越多。

首先，我们过滤doublets。我们将它保存为一个新的ArchRProject，以供本步教程使用，但你也可以覆盖你原来的ArchRProject对象：

> projHeme2 <- filterDoublets(projHeme1) #this step has problem for windows,using the codes below:
Filtering 410 cells from ArchRProject!
    scATAC_BMMC_R1 : 243 of 4932 (4.9%)
    scATAC_CD34_BMMC_R1 : 107 of 3275 (3.3%)
    scATAC_PBMC_R1 : 60 of 2454 (2.4%)

注意！！！这一步如果你按照官网教程运行，仍然可能出错！！！报错信息:Error in as.POSIXct.numeric(e) : 'origin' must be supplied。我google了报错原因：https://github.com/GreenleafLab/ArchR/issues/217，根据这网页里的讨论内容，这一步在windows里运行基本上都会出错。有人提出了直接修改filterDoublets()函数的代码，使用新function代码，再运行上面的步骤，就不会报错了！具体代码我不贴上来了，很长，大家可以在网站里自行copy运行。（然而我在Linux系统里运行，使用原代码就没有任何报错。。。）

> projHeme2

class: ArchRProject 
outputDirectory: D:\yanfang\HemeTutorial 
samples(3): scATAC_BMMC_R1 scATAC_CD34_BMMC_R1 scATAC_PBMC_R1
sampleColData names(1): ArrowFiles
cellColData names(17): Sample TSSEnrichment ... bioNames bioNames2
numberOfCells(1): 10251
medianTSS(1): 16.845
medianFrags(1): 2998

如果你想过滤更多的细胞，要使用更多的filterRatio值：

> projHemeTmp <- filterDoublets(projHeme1, filterRatio = 1.5)
Filtering 614 cells from ArchRProject!
    scATAC_BMMC_R1 : 364 of 4932 (7.4%)
    scATAC_CD34_BMMC_R1 : 160 of 3275 (4.9%)
    scATAC_PBMC_R1 : 90 of 2454 (3.7%)

因为这个Tmp只是举个例子，所以我们要把它删掉：

> rm(projHemeTmp)
上一篇下一篇

猜你喜欢

热点阅读