使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第十四章)

本文首發(fā)于我的個(gè)人博客， http://xuzhougeng.top/

往期回顧:

第十四章 ArchR的足跡分析

轉(zhuǎn)錄因子(Transcripts factor, TF)足跡分析使得我們能夠預(yù)測特定位點(diǎn)中TF的精確結(jié)合位置。這是因?yàn)樵撐恢帽籘F結(jié)合避免了轉(zhuǎn)座酶的切割，而TF結(jié)合位點(diǎn)的鄰近位置處于開放狀態(tài)。

footprintingSchematic

理想情況下，TF足跡分析需要在單個(gè)位置上分析從而確定TF的準(zhǔn)確結(jié)合位置。但實(shí)際上，這需要非常高的測序深度，甚至超過混池ATAC-seq或者scATAC-seq的所有數(shù)據(jù)。為了解決這個(gè)問題，我們可以把和待預(yù)測的TF結(jié)合相關(guān)的Tn5插入位置進(jìn)行合并。例如，我們可以提取所有包含CTCF motif的peak，制作一個(gè)全基因組的CTCF的聚合TF足跡。

為了保證足跡的可靠性，我們需要確保能夠可靠的預(yù)測出目標(biāo)TF所對(duì)應(yīng)的結(jié)合位點(diǎn)。ArchR使用自帶的addMotifAnnotations()函數(shù)對(duì)peak區(qū)域進(jìn)行搜索，尋找能夠匹配的DNA序列。考慮到motif的簡并性，無法保證每個(gè)motif都有足夠的peak。添加到ArchRProject的motif注釋以二值矩陣表示(0=無motif, 1=有motif)。一旦你有了這些motif注釋，ArchR使用getFootPrints()函數(shù)分析足跡，它以一個(gè)ArchRProject對(duì)象和一個(gè)GenomicRanges對(duì)象(記錄motif的位置)作為輸入?？梢允褂?code>getPositions()函數(shù)從ArchRProject中提取這些位置。之后足跡可以使用plotFootprints()函數(shù)可視化。

或許更重要的是，ArchR的足跡分析能夠抵消已知的Tn5插入序列偏好性。ArchR使用一個(gè)hexmer位置頻率矩陣和一個(gè)目標(biāo)Tn5插入位置上的k-mer頻率矩陣來實(shí)現(xiàn)該功能。

foot printing Methods

最終，該流程輸出考慮到Tn5插入偏好性的足跡圖。

ArchR支持motif足跡分析和用戶提供特征的足跡分析，在后續(xù)都會(huì)討論。

14.1 motif足跡分析

由于教程用到的數(shù)據(jù)集比較小，因此利用該數(shù)據(jù)集得到足跡并不是那么的清晰。使用更大的數(shù)據(jù)會(huì)得到較小變異的足跡。

在分析足跡時(shí)，我們要先獲取和motif相關(guān)的所有位置。這可以通過getPositions函數(shù)完成。該函數(shù)有一個(gè)可選參數(shù)， name，用于傳入peakAnnotation對(duì)象中我們想要獲取位置的變量名。如果name=NULL, 那么ArchR會(huì)使用peakAnnotation槽(slot)的第一個(gè)條目(entry)。在下面的例子中，我們沒有指定name, ArchR使用的第一個(gè)條目為CIS-BP motifs.

motifPositions <- getPositions(projHeme5)

這會(huì)創(chuàng)建一個(gè)GRangesList對(duì)象，每個(gè)TF motif以不同的GRanges對(duì)象進(jìn)行區(qū)分。

motifPositions
# GRangesList object of length 870:
# $TFAP2B_1
# GRanges object with 16773 ranges and 1 metadata column:
# seqnames ranges strand | score
# |
# [1] chr1 852468-852479 + | 8.17731199746359
# [2] chr1 873916-873927 + | 8.32673820065588
# [3] chr1 873916-873927 - | 8.32673820065588
# [4] chr1 896671-896682 + | 9.96223327271814
# [5] chr1 896671-896682 - | 8.92408377606486
# … … … … . …
# [16769] chrX 153991101-153991112 + | 8.39549159740639
# [16770] chrX 154299568-154299579 + | 8.90119825654299
# [16771] chrX 154664929-154664940 - | 8.16690864294221
# [16772] chrX 154807684-154807695 + | 9.57636587154549
# [16773] chrX 154807684-154807695 - | 10.6117355833828
# ——-
# seqinfo: 23 sequences from an unspecified genome; no seqlengths
#
# …
# <869 more elements>

我們提取部分感興趣的TF motifs用于展示。在我們提取"EBF1"會(huì)附帶"SREBF1", 因此我們需要用%ni%顯式將其過濾。%ni%函數(shù)是R自帶函數(shù)%in%的相反函數(shù)。

motifs <- c("GATA1", "CEBPA", "EBF1", "IRF4", "TBX21", "PAX5")
markerMotifs <- unlist(lapply(motifs, function(x) grep(x, names(motifPositions), value = TRUE)))
markerMotifs <- markerMotifs[markerMotifs %ni% "SREBF1_22"]
markerMotifs
# [1] “GATA1_383” “CEBPA_155” “EBF1_67” “IRF4_632” “TBX21_780” “PAX5_709”

為了準(zhǔn)確找到TF足跡，我們需要大量的reads。因此，細(xì)胞需要進(jìn)行分組生成擬混池ATAC-seq譜才能用于TF足跡分析。這些擬混池譜之前在peak鑒定時(shí)就已經(jīng)保存為分組覆蓋文件。 如果沒有在ArchRProject添加分組覆蓋信息，則運(yùn)行如下命令

projHeme5 <- addGroupCoverages(ArchRProj = projHeme5, groupBy = "Clusters2")

在計(jì)算分組覆蓋度后，我們可以為之前getFootprints()挑選的一組標(biāo)記motif計(jì)算足跡。即便ArchR已經(jīng)優(yōu)化了足跡分析流程，我們也建議先對(duì)一部分motif分析足跡，而不是直接分析所有motif。 我們通過positions參數(shù)來選擇motif。

seFoot <- getFootprints(
  ArchRProj = projHeme5, 
  positions = motifPositions[markerMotifs], 
  groupBy = "Clusters2"
)

當(dāng)我們獲取了這些足跡，我們可以使用plotFootprints()函數(shù)進(jìn)行展示。該函數(shù)能夠同時(shí)以多種方式對(duì)足跡進(jìn)行標(biāo)準(zhǔn)化。下一節(jié)，我們會(huì)討論標(biāo)準(zhǔn)化和實(shí)際的足跡圖。

14.2 Tn5偏好的足跡標(biāo)準(zhǔn)化

使用ATAC-seq數(shù)據(jù)分析TF足跡的一大挑戰(zhàn)就是Tn5轉(zhuǎn)座酶的插入序列偏好性，這會(huì)導(dǎo)致TF足跡的錯(cuò)誤分類。為了降低Tn5插入偏好性的影響，ArchR識(shí)別每個(gè)Tn5插入位置附近的k-mer序列(k由用戶提供，默認(rèn)是6).

對(duì)于該項(xiàng)分析，ArchR為每個(gè)擬混池識(shí)別單堿基分辨率的Tn5插入位點(diǎn)，將這些1-bp位點(diǎn)調(diào)整為k-bp窗口（-k/2和+（k/2-1）bp），然后使用Biostrings包中的oligonucleotidefrequency(w=k, simplify.as="collapse")函數(shù)創(chuàng)建k-mer頻率表。然后，ArchR使用與BSgenome相關(guān)的基因組文件，以相同的函數(shù)計(jì)算出全基因組范圍預(yù)期的k-mers。

為了計(jì)算擬混池足跡的插入偏差，ArchR創(chuàng)建了一個(gè)k-mer頻率矩陣，該矩陣表示為從motif中心到窗口+/-N bp（用戶定義，默認(rèn)為250 bp）的所有可能k-mer。然后，遍歷每個(gè)motif位點(diǎn)，ArchR將定位的k-mer填充到k-mer頻率矩陣中。然后在全基因組范圍內(nèi)計(jì)算每個(gè)motif位置。利用樣本的k-mer頻率表，ArchR可以通過將k-mer位置頻率表乘以觀察/期望 Tn5 k-mer頻率來計(jì)算預(yù)期的Tn5插入。

以上所有這些發(fā)生在plotFootprints()函數(shù)中。

14.2.1 減去Tn5偏好

一個(gè)標(biāo)準(zhǔn)化方式就是從足跡信號(hào)中減去Tn5偏好。該標(biāo)準(zhǔn)化方法通過設(shè)置plotFootprints()的normMethod = "Subtract"實(shí)現(xiàn)

plotFootprints(
  seFoot = seFoot,
  ArchRProj = projHeme5, 
  normMethod = "Subtract",
  plotName = "Footprints-Subtract-Bias",
  addDOC = FALSE,
  smoothWindow = 5
)

默認(rèn)，這些圖保存在ArchRProject的outputDirectory。如果你需要繪制所有motif, 可以將其返回為ggplot2對(duì)象，需要注意這個(gè)ggplot對(duì)象會(huì)非常大。下面是一個(gè)從motif足跡中減去Tn5偏好信號(hào)的結(jié)果

Footprints-Subtract-Bias

14.2.2 除以Tn5偏好

第二種標(biāo)準(zhǔn)化方法就是將足跡除以Tn5偏好信號(hào)。該標(biāo)準(zhǔn)化方法通過設(shè)置plotFootprints()的normMethod = "Divide"實(shí)現(xiàn)

plotFootprints(
  seFoot = seFoot,
  ArchRProj = projHeme5, 
  normMethod = "Divide",
  plotName = "Footprints-Divide-Bias",
  addDOC = FALSE,
  smoothWindow = 5
)

下面是一個(gè)從motif足跡中除以Tn5偏好信號(hào)的結(jié)果

Footprints-Divide-Bias

14.2.3 無Tn5偏好標(biāo)準(zhǔn)化的足跡

盡管我們高度推薦將足跡根據(jù)Tn5序列插入偏好性進(jìn)行標(biāo)準(zhǔn)化，當(dāng)然你可以通過設(shè)置plotFootprints()的normMethod = "None"來省去標(biāo)準(zhǔn)化。

plotFootprints(
  seFoot = seFoot,
  ArchRProj = projHeme5, 
  normMethod = "None",
  plotName = "Footprints-No-Normalization",
  addDOC = FALSE,
  smoothWindow = 5
)

下面是一個(gè)未經(jīng)標(biāo)準(zhǔn)化的motif足跡結(jié)果

Footprints-No-Normalization

14.3 特征足跡

除了motif足跡分析，ArchR還允許用戶分析任意定義的特征數(shù)據(jù)集。為了對(duì)功能進(jìn)行說明，我們將會(huì)使用plotFootprints()函數(shù)創(chuàng)建TSS插入譜(在之前數(shù)據(jù)質(zhì)量控制一節(jié)中引入)。一個(gè)TSS插入譜本質(zhì)上就是特殊的足跡。

我們?cè)谥靶」?jié)討論過，足跡會(huì)用到來源于擬混池重復(fù)的分組覆蓋文件。我們?cè)谥皃eak鑒定時(shí)創(chuàng)建過這些文件。如果你沒有在ArchRProject加入分組覆蓋信息, 那么需要運(yùn)行如下代碼

projHeme5 <- addGroupCoverages(ArchRProj = projHeme5, groupBy = "Clusters2")

我們接著創(chuàng)建一個(gè)沒有經(jīng)過Tn5偏好性校正的TSS插入譜。和之前分析一個(gè)主要不同是，我們?cè)O(shè)置了flank=2000, 將TSS向前向后分別延伸2000 bp.

seTSS <- getFootprints(
  ArchRProj = projHeme5, 
  positions = GRangesList(TSS = getTSS(projHeme5)), 
  groupBy = "Clusters2",
  flank = 2000
)

我們接著用plotFootprints()對(duì)每個(gè)細(xì)胞分組繪制TSS插入譜。

plotFootprints(
  seFoot = seTSS,
  ArchRProj = projHeme5, 
  normMethod = "None",
  plotName = "TSS-No-Normalization",
  addDOC = FALSE,
  flank = 2000,
  flankNorm = 100
)
TSS-No-Normalization
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

友情鏈接更多精彩內(nèi)容