偷偷问一下，关注了吗？
内容获取
1、购买打包合集（《KS科研分享与服务》付费内容打包集合），价格感人，可以加入微信VIP群（答疑交流群，甚至有小伙伴觉得群比代码更好），可以获取建号以来所有内容，群成员专享视频教程，提前更新，其他更多福利！
2、《KS科研分享与服务》公众号有QQ群，进入门槛是20元（完全是为了防止白嫖党，请理解），请考虑清楚。群里有免费推文的注释代码和示例数据（终身拥有），没有付费内容，群成员福利是购买单个付费内容半价！
需要者详情请联系作者（非需要者勿扰，处理太费时间）：
“人有悲欢离合，月游阴晴圆缺，此事古难全”。与scRNA一样，你可能会遇到各种情况的数据load（跟着Cell学单细胞转录组分析(二):单细胞转录组测序文件的读入及Seurat对象构建）。如果自己测序或者自己跑上游的，那么得到的标准的结果，就可以用标准的流程了，可是哪能都那么完美，有时候我们跑不了上游，或者别人提供的数据库没有上游文件或者出错，只有下游文件，但是不标准怎么办呢？我们这里介绍signac分析scATAC下游，signac是seurat团队开发的包，所以数据结构和流程类似scRNA，那么如果你分析过scRNA，就会很好上手。完整内容及视频教程已发布微信VIP群，请自行下载！signac官网有详细的教程，可以阅读！
signac官网教程：https://stuartlab.org/signac/
Signac currently supports the following features:
Calling peaks
Quantifying per-cell counts in different genomic regions
Calculating single-cell QC metrics
Dimensional reduction, visualization, and clustering
Identifying cell-type-specific peaks
Visualizing ‘pseudo-bulk’ coverage tracks
Integration of multiple single-cell datasets
Integration with single-cell RNA-seq datasets
Sequence motif enrichment analysis
Transcription factor footprinting analysis
Linking peaks to correlated genes
Parallelization through the future package
Seamless interface with Seurat, SeuratWrappers, SeuratDisk, and SeuratData functionality
Interoperability with Bioconductor tools
首先安装signac:
setRepositories(ind=1:3) # needed to automatically install Bioconductor dependenciesinstall.packages("Signac")
为什么我们要用其他数据演示流程，官网教程具有迷惑性，会让你失去判断力，可以搜一下signac的教程，90%的博主都是直接运行一遍官网的内容，那有什么意义呢？你都不知道数据怎么来的，自己中间会出现什么问题，以及官网的数据过于完美，导致自己数据判断不足等等。接下来是单个样本标准数据的读取，创建seurat obj，很熟悉吧。
signac输入的数据需要3个：
1、表达矩阵，也有三个标准文件，只不过对于scATAC就是peak矩阵，矩阵也可以读取filtered_peak_bc_matrix.h5.h5文件
2、：fragments.tsv.gz | fragments.tsv.gz.tbi，这是ATAC分析的一个关键文件，是所有cell中片段信息文件，是最大的，一般G起步。
3、singlecell.csv或者per_barcode_metrics.csv（前者可能是大部分遇到的），是细胞metadata，用于质控等。
#加载相关R包library(Signac)library(Seurat)library(ggplot2)library(dplyr)library(harmony)library(hdf5r)library(GenomeInfoDb)
# drwxr-xr-x 5 ks_ts staff 160B Dec 4 19:57 filtered_peak_bc_matrix# -rw-r--r-- 1 ks_ts staff 38M Dec 4 19:58 filtered_peak_bc_matrix.h5# drwxr-xr-x 5 ks_ts staff 160B Dec 4 19:57 filtered_tf_bc_matrix# -rw-r--r-- 1 ks_ts staff 4.5M Dec 4 19:58 filtered_tf_bc_matrix.h5# -rw-r--r-- 1 ks_ts staff 1.5G Dec 4 21:00 fragments.tsv.gz# -rw-r--r-- 1 ks_ts staff 995K Dec 4 19:58 fragments.tsv.gz.tbi# -rw-r--r-- 1 ks_ts staff 9.0M Dec 4 19:57 peak_annotation.tsv# -rw-r--r-- 1 ks_ts staff 33M Dec 4 19:58 peak_motif_mapping.bed# -rw-r--r-- 1 ks_ts staff 3.9M Dec 4 19:57 peaks.bed# -rw-r--r-- 1 ks_ts staff 7.3M Dec 4 19:59 possorted_bam.bam.bai# drwxr-xr-x 5 ks_ts staff 160B Dec 4 19:58 raw_peak_bc_matrix# -rw-r--r-- 1 ks_ts staff 51M Dec 4 19:59 raw_peak_bc_matrix.h5# -rw-r--r-- 1 ks_ts staff 27M Dec 4 19:58 singlecell.csv# -rw-r--r-- 1 ks_ts staff 1.1K Dec 4 19:58 summary.csv# -rw-r--r-- 1 ks_ts staff 29K Dec 4 19:57 summary.json# -rw-r--r-- 1 ks_ts staff 2.7M Dec 4 19:57 web_summary.html#为了方便，这里我们读取矩阵的时候读那个.h5文件HC.counts <- Read10X_h5('./SRR_HC/outs/filtered_peak_bc_matrix.h5')HC.meta <- read.table(file = "./SRR_HC/outs/singlecell.csv",stringsAsFactors = FALSE,sep = ",",header = TRUE,row.names = 1)[-1, ]HC.meta <- HC.meta[HC.meta$passed_filters > 500, ]#对于低计数的cell进行过滤HC.frag <- CreateChromatinAssay(counts = HC.counts,sep = c(":", "-"),fragments = "./SRR_HC/outs/fragments.tsv.gz",min.cells = 10,min.features = 200)#create fragment objectsHC <- CreateSeuratObject(counts = HC.frag,assay = "peaks",meta.data = HC.meta, project = 'HC')#to seurat obj
其余两个样本如法炮制读入：
#读入其余两个样本文件#load AA groupAA.counts <- Read10X_h5('./SRR_AA/outs/filtered_peak_bc_matrix.h5')AA.meta <- read.table(file = "./SRR_AA/outs/singlecell.csv",stringsAsFactors = FALSE,sep = ",",header = TRUE,row.names = 1)[-1, ]AA.meta <- AA.meta[AA.meta$passed_filters > 500, ]#对于低计数的cell进行过滤AA.frag <- CreateChromatinAssay(counts = AA.counts,sep = c(":", "-"),fragments = "./SRR_AA/outs/fragments.tsv.gz",min.cells = 10,min.features = 200)#create fragment objectsAA <- CreateSeuratObject(counts = AA.frag,assay = "peaks",meta.data = AA.meta, project = 'AA')#to seurat obj#load SD groupSD.counts <- Read10X_h5('./SRR_SD/outs/filtered_peak_bc_matrix.h5')SD.meta <- read.table(file = "./SRR_SD/outs/singlecell.csv",stringsAsFactors = FALSE,sep = ",",header = TRUE,row.names = 1)[-1, ]SD.meta <- SD.meta[SD.meta$passed_filters > 500, ]#对于低计数的cell进行过滤SD.frag <- CreateChromatinAssay(counts = SD.counts,sep = c(":", "-"),fragments = "./SRR_SD/outs/fragments.tsv.gz",min.cells = 10,min.features = 200)#create fragment objectsSD <- CreateSeuratObject(counts = SD.frag,assay = "peaks",meta.data = SD.meta, project = 'SD')#to seurat obj
如果没有filtered_peak_bc_matrix.h5，那么矩阵就需要用到filtered_peak_bc_matrix文件夹里面的3个文件了！其他不变！
#来源于signac官网# counts <- Matrix::readMM("filtered_peak_bc_matrix/matrix.mtx")# barcodes <- readLines("filtered_peak_bc_matrix/barcodes.tsv")# peaks <- read.table("filtered_peak_bc_matrix/peaks.bed", sep="\t")# peaknames <- paste(peaks$V1, peaks$V2, peaks$V3, sep="-")## colnames(counts) <- barcodes# rownames(counts) <- peaknames
那么如果提供的数据既没有没有filtered_peak_bc_matrix.h5，也没有三个文件，只有fragments.tsv.gz也可以创建矩阵，用于signac分析！
如果既没有filtered_peak_bc_matrix.h5，也没有三个文件，只有fragments.tsv.gz也可以创建矩阵total_counts <- CountFragments('./SRR_HC/outs/fragments.tsv.gz')barcodes <- total_counts[total_counts$frequency_count > 500, ]$CBfrags <- CreateFragmentObject(path = './SRR_HC/outs/fragments.tsv.gz', cells = barcodes)peaks <- CallPeaks(frags, macs2.path = "/Users/aaa/Library/Python/3.9/bin/macs2")counts <- FeatureMatrix(fragments = frags, features = peaks, cells = barcodes)# Quantify fragments in each peak
关注我们获取精彩内容：
关注不迷路：扫描下面二维码关注公众号！
B站视频号链接：https://space.bilibili.com/471040659?spm_id_from=333.1007.0.0
关注 KS科研分享与服务，
认清正版优质内容和服务！
优质内容持续输出，物超所值！
合作联系：ks_account@163.com
