生物信息學習的正確姿勢

NGS系列文章包括NGS基礎、在線繪圖、轉錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫學設計實驗GEO數據分析 (step-by-step)）、批次效應處理等內容。

ATAC-seq即transposase-accessible chromatin using sequencing，是一種檢測染色體開放區域的技術。該方法由斯坦福大學Greenleaf實驗室在2013年首次發表[1]，兩年后該實驗室又發表基于單細胞的ATAC-seq技術[2]。ATAC-seq相比于之前的FaiRE-seq和DNase-seq，ATAC-seq用Tn5轉座酶對染色體開放區域進行剪切，并加上adaptors序列（圖1的紅藍片段）。最后得到的DNA片段，包括了開放區域的剪切片段，以及橫跨一個或多個核小體的長片段。

圖1. ATAC-seq[1]

根據片段長度可以分為Fragments in nucleosome-free regions(<147 base pairs)和Fragments flanking a single nucleosome (147~294 base pairs), 以及更長的多核片段。片段長度分布如下圖，沒有跨越核小體的小片段最多，其次是單核片段，依次遞減。
? 圖2. DNA片段長度分布

scATAC-seq建庫原理

以10x 建庫方法為例[5]，比較scATAC-seq 和scRNA-seq建庫方法的異同。
二者都用膠珠（GEMs）的方法，不一樣的是ATAC膠珠上的序列中不用UMI，因為基因組只有一對序列，無需像RNA一樣定量。另外序列末端用接頭引物Read 1N代替PolyT。

scRNA-seq通過結合cDNA的PolyA尾進行擴增，而scATAC-seq的DNA片段沒有PolyA尾，取而代之的是Tn5酶轉座剪切時插入的adaptors片段，可以與膠珠上的Read 1N序列互補。

DNA片段接上膠珠后，在另一端加Read2和Sample index序列。在此之前，scRNA-seq需要將cDNA酶切至合適的片段長度，而scATAC-seq的片段不進行打碎，接上Sample index和P7序列后進行擴增。

最后上機測序。scRNAseq如果是3‘單端測序，Read2讀取最近的100bp讀長，而Read1只讀取16bp的細胞barcode序列和10bp的UMI序列，共26bp。scATAC-seq則用雙末端測序，讀長一般不低于45bp[3]。

scATAC-seq最后可以得到4個原始文件：

其中I1/2分別是barcode和sample index，R1/2是目的片段的雙末端。

10x提供cellranger軟件對原始數據進行初步分析，如質控，比對，peak calling等。

$ cellranger-atac count --id=sample345 \ --reference=/opt/refdata-cellranger-atac-GRCh38-1.2.0 \ --fastqs=/home/jdoe/runs/HAWT7ADXX/outs/fastq_path \ --sample=mysample \ --localcores=8 \ --localmem=64

質控方法

重要指標：

Fraction of fragments overlapping any targeted region
覆蓋注釋區域的片段比例。目標區域包括TSS，DNase HS，增強子和啟動子等。一般要求大于55%。

Fraction of transposition events in peaks in cell barcodes
轉座位點位于peaks區域的比例。一般大于25%。如果比例過小，有可能是細胞狀態不好，或者測序深度不夠。

Enrichment score of transcription start sites(TSS)

轉錄起始位點的標準化分數。閾值選擇根據基因組而定。人類一般選擇TSS分數大于5或者6的細胞樣本[3]。如果分數太低說明細胞染色質結構瓦解，或者實驗過程細胞裂解方式不當。

TSS標準化分數的計算方法不同，可能導致閾值偏差。10X的標準化方法是cut sites數除以最小值，而ENCODE的標準是在cut sites數除以兩個末端各100bp的cut site數的平均值。由于一般越靠近中間數值越大，ENCODE的標準化分數整體比10x的分數小一些。

其他指標：

• Fraction of read pairs with a valid barcode > 75%;

• Q30 bases in R1 > 65%;

• Q30 bases in R2 > 65%;

• Q30 bases in barcode > 65%;

• Q30 bases in Sample Index > 90%;

• Estimated Number of Cells 500~10000 +- 20%;

• Median fragments per cell barcode > 500;

• Fragments in nucleosomefree regions >55%;

• Fraction of total read pairs mapped confidently to genome (>30 mapq) >80%;

• Fraction of total read pairs in mitochondria and in cell barcodes < 40%;

下游分析（以Signac為例）

Signac包由Seurat同一團隊開發，獨立于Seurat包，在2020年8月開始發布在GitHub上。目前仍是1.0.0版本[4]。

我們可以用10x官網的PBMC數據做演示。

首先加載所需R包：

library(Signac) library(Seurat) library(GenomeInfoDb) library(EnsDb.Hsapiens.v75) library(ggplot2) library(patchwork) set.seed(1234)

加載peaks, 細胞注釋和片段分布數據，并創建object。這個object和Seurat object類似，只是在assay里多了peaks等信息。

counts<-Read10X_h5(filename= "atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5") metadata <- read.csv(file = "atac_v1_pbmc_10k_singlecell.csv",header = TRUE,row.names = 1 ) chrom_assay <- CreateChromatinAssay(counts = counts,sep = c(":", "-"),genome = 'hg19',fragments = 'atac_v1_pbmc_10k_fragments.tsv.gz',min.cells = 10,min.features = 200 ) pbmc <- CreateSeuratObject(counts = chrom_assay,assay = "peaks",meta.data = metadata )

總共8728個細胞，87405個features。features不是基因，是基因組的注釋區域，如啟動子，增強子等。

pbmc[['peaks']] ## ChromatinAssay data with 87405 features for 8728 cells ## Variable features: 0 ## Genome: hg19 ## Annotation present: FALSE ## Motifs present: FALSE ## Fragment files: 1

加載注釋:

extract gene annotations from EnsDb annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75) # change to UCSC style since the data was mapped to hg19 seqlevelsStyle(annotations) <- 'UCSC' genome(annotations) <- "hg19" # add the gene information to the object Annotation(pbmc) <- annotations

TSS和blacklist比例計算。

# compute nucleosome signal score per cell pbmc <- NucleosomeSignal(object = pbmc) # compute TSS enrichment score per cell pbmc <- TSSEnrichment(object = pbmc, fast = FALSE) # add blacklist ratio and fraction of reads in peaks pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100 pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments pbmc$high.tss <- ifelse(pbmc$TSS.enrichment > 2, 'High', 'Low') TSSPlot(pbmc, group.by = 'high.tss') + NoLegend()

質控。Signac包用5個指標做過濾：

TSS富集分數，大于2；
blacklist比例，小于0.05；
纏繞核小體的片段與非核小體片段（< 147bp）的比例，小于4；
匹配peaks區域的片段比例，大于15%；
匹配peaks區域的片段數，大于3000，小于20000。

pbmc <- subset(x = pbmc,subset = peak_region_fragments > 3000 &peak_region_fragments < 20000 &pct_reads_in_peaks > 15 &blacklist_ratio < 0.05 &nucleosome_signal < 4 &TSS.enrichment > 2 ) pbmc ## An object of class Seurat ## 87405 features across 7060 samples within 1 assay ## Active assay: peaks (87405 features, 0 variable features)

降維聚類。

pbmc <- RunTFIDF(pbmc) pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0') pbmc <- RunSVD(pbmc) pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30) pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30) pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3) DimPlot(object = pbmc, label = TRUE) + NoLegend()

創建基因活性矩陣。之前的聚類區域所用的features是peaks，為了展示不同分群基因活性的差異，首先要創建一個類似RNA表達的矩陣。用基因加上游2000bp區域的比對片段數代表該基因的活性。

gene.activities <- GeneActivity(pbmc) # add the gene activity matrix to the Seurat object as a new assay and normalize it pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities) pbmc <- NormalizeData(object = pbmc,assay = 'RNA',normalization.method = 'LogNormalize',scale.factor = median(pbmc$nCount_RNA) )

展示marker基因活性：

DefaultAssay(pbmc) <- 'RNA'FeaturePlot(object = pbmc,features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'),pt.size = 0.1,max.cutoff = 'q95',ncol = 3 )

與scRNA-seq數據的整合分析。單細胞轉錄組數據地址

# Load the pre-processed scRNA-seq data for PBMCs pbmc_rna <- readRDS("/home/stuartt/github/chrom/vignette_data/pbmc_10k_v3.rds") transfer.anchors <- FindTransferAnchors(reference = pbmc_rna,query = pbmc,reduction = 'cca' ) predicted.labels <- TransferData(anchorset = transfer.anchors,refdata = pbmc_rna$celltype,weight.reduction = pbmc[['lsi']],dims = 2:30 ) pbmc <- AddMetaData(object = pbmc, metadata = predicted.labels)

尋找細胞分群特異的peaks。

# change back to working with peaks instead of gene activities DefaultAssay(pbmc) <- 'peaks' da_peaks <- FindMarkers(object = pbmc,ident.1 = "CD4 Naive",ident.2 = "CD14 Mono",min.pct = 0.2,test.use = 'LR',latent.vars = 'peak_region_fragments' ) plot1 <- VlnPlot(object = pbmc,features = rownames(da_peaks)[1],pt.size = 0.1,idents = c("CD4 Memory","CD14 Mono") ) plot2 <- FeaturePlot(object = pbmc,features = rownames(da_peaks)[1],pt.size = 0.1 ) plot1 | plot2

展示基因在不同細胞類型的開放程度：

# set plotting order levels(pbmc) <- c("CD4 Naive","CD4 Memory","CD8 Naive","CD8 Effector","DN T","NK CD56bright","NK CD56Dim","pre-B",'pro-B',"pDC","DC","CD14 Mono",'CD16 Mono')CoveragePlot(object = pbmc,region = rownames(da_peaks)[1],extend.upstream = 40000,extend.downstream = 20000 )

此外還有其他分析，如TF footprinting等。footprinting顧名思義是指轉錄因子留下的印記，由于Tn5酶不能剪切到TF結合的區域，所以footprinting圖相對與TSS圖，中間有“凹陷”，凹陷的程度根據TF結合的時間確定[6]。

總的來說，Signac包是一個親民實用的工具。雖然有一些不足的地方，如染色體的注釋目前只能選擇UCSC的，不能選Ensemble。

參考文獻：

[1]?Buenrostro, J., Giresi,?P., Zaba, L.?et al.?Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position.?Nat Methods?10,?1213–1218 (2013). https://doi.org/10.1038/nmeth.2688

[2]?Buenrostro, J., Wu, B., Litzenburger, U.?et al.?Single-cell chromatin accessibility reveals principles of regulatory variation.?Nature?523,?486–490 (2015).?https://doi.org/10.1038/nature14590

[3]?https://www.encodeproject.org/atac-seq/

[4]?https://satijalab.org/signac/index.html

[5]https://assets.ctfassets.net/an68im79xiti/Cts31zFxXFXVwJ1lzU3Pc/fe66343ffd3039de73ecee6a1a6f5b7b/CG000202_TechnicalNote_InterpretingCellRangerATACWebSummaryFiles_Rev-A.pdf

[6]?Sung, M., Baek, S. & Hager, G. Genome-wide footprinting: ready for prime time?.?Nat Methods?13,?222–228 (2016). https://doi.org/10.1038/nmeth.3766

往期精品(點擊圖片直達文字對應教程)

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結

以上是生活随笔為你收集整理的scATAC-seq建库原理，质控方法和新R包Signac的使用的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。