Socrates是一款用于分析snATAC数据的R语言包，但是该包发表后并没有引起多大水花，最直观的就是，Github上该包的Issues只有10条，其中两条还是本人最近贡献的。这里附上Socrates的Github链接：https://github.com/plantformatics/Socrates。
虽然这个包没啥人用，但是由于本人爱折腾，最终还是尝试了一把，接下来讲一些使用该包的技巧，绝对是全网首发！！！
第一步：读取输入文件
输入文件总共需要3个，分别是Tn5结合位点文件，基因组的gff注释文件，基因组染色体长度文件。
但是通常情况下，我们上游通过cellranger-atac得到的是fragment文件，这是需要在读取时设置is.fragment=T。
library(Socrates)
#读取文件
bed="SampleName.fragments.tsv.gz"
ann="genome.gff3"
chr="genome.chr.bed"
#建立obj，注意Socrates包默认的bed为包含正负链信息的单碱基Tn5插入位点，但我们的fragment文件是从cellranger-atac中得到的，因此在读取时需要设置is.fragment=T
obj=loadBEDandGenomeData(bed,ann,chr,is.fragment=T)
这里仍然存在一个小问题就是如果fragment文件过大的话，是无法读取的，这里需要自行将fragment文件转换成单碱基的Tn5插入位点文件，直接进行读取。至于转换的方式，则是根据cellranger-atac产生fragment文件的方式来逆向还原Tn5插入位点的位置。
cellranger-atac count的结果中会生成fragments.gz文件，其中每一个fragment是由两个独立的转座事件创造的，也就是说fragment的两端正好是两次Tn5酶插入的位点。每一个特异的fragment会产生多条冗余的reads，这些冗余的reads被折叠成单个fragment记录。
Fragment的位置信息是从bam文件中获得的。fragment的起点是将bam中的reads比对的位置最左侧向前移动4bp，终点是将reads比对位置的最右侧向后移动5bp。
因此将fragment文件中每条fragment的起点和终点都分别看做一个Tn5结合位点即可。
将fragment文件转换成Tn5结合位点文件
unpigz -p 24 -k SampleName.fragments.tsv.gz &&
echo "unpigz done" &&
awk '{print $1 "\t" $2 "\t" $2+1 "\t" $4 "\t" "+" }' SampleName.fragments.tsv > SampleName.start &&
echo "start done" &&
awk '{print $1 "\t" $3 "\t" $3-1 "\t" $4 "\t" "-" }' SampleName.fragments.tsv > SampleName.end &&
echo "end done" &&
cat SampleName.start SampleName.end | sort -k1,1 -k2,2n | uniq > SampleName.Tn5.sites.txt &&
echo "cat done" &&
pigz -p 24 SampleName.Tn5.sites.txt &&
echo "pigz done" &&
第二步：进行peak calling并过滤低质量的细胞
#call peak，鉴定ACRs
##output指定输出peak文件和peak的前缀，tempdir指定peak的输出目录
obj=callACRs(obj,genomesize=1.46e10,shift=-75,extsize=100,fdr=0.05,output="SampleName",tempdir="./SampleName")
#生成metadata
##注意这里生成metadata的时候需要提前call peak，否则不会计算frip，在后续的过滤过程中会出问题
obj=buildMetaData(obj,tss.window=3000)
head(obj$meta)
#过滤低质量的细胞
#findCells后会生成多个meta，分别为meta.v1，meta.v2，meta.v3，他们分别对应pdf图中的三次过滤，其中meta.v3为最终的过滤结果。
##min.cells和max.cells代表最低和最高的细胞数目
##min.tn5代表Tn5位点处最低的read depth
##filt.tss=T代表根据tss进行过滤
##tss.min.freq=0.2和tss.z.thresh代表reads比对到TSS位点处的比例超过0.2，同时z-score的阈值为3
##filt.frip=T代表根据frip进行过滤
##frip.min.freq=0.2代表细胞的frip需高于0.2并且z-score的阈值为3
obj=findCells(obj,doplot=T,prefix="S4-1",min.cells=1000,max.cells=15000,min.tn5=1000,filt.tss=T,tss.min.freq=0.2,tss.z.thresh=3,filt.frip=T,frip.min.freq=0.2,frip.z.thresh=3)
这里在经过质控过滤后，Socrates并不会丢弃那些被过滤掉的细胞，而是增加多个meta来进行区分。
第三步：建立Socrates对象并进行下游分析
#生成Socrates对象
##filtered=T，指定用过滤后的细胞建立matrix，windows=1000指定基因组划窗的大小为500，peaks=F不适用peak文件进行定量，如果设置peaks=T，则对peak进行定量，并覆盖windows=500的参数。
obj=generateMatrix(obj,filtered=T,windows=500,peaks=F,verbose=T)
obj=generateMatrix(obj,filtered=T,peaks=T,verbose=T)
#将该对象转换为Socrates格式
soc.obj=convertSparseData(obj)
saveRDS(obj,file="QC_object.rds")
saveRDS(soc.obj,file="Socrates_object.rds")
#对Socrates对象进行过滤
cell.counts=Matrix::colSums(soc.obj$counts)
site.freq=Matrix::rowMeans(soc.obj$counts)
layout(matrix(c(1:2),ncol=2))
par(mar=c(3,3,1,1))
plot(density(cell.counts),main="Log10 cell counts",log="x")
abline(v=1000,col="red")
plot(density(site.freq),main="peak accessibility frequency",log="x")
#对Socrates对象进行过滤
##min.c设置每个细胞中accessible ACRs的最小值
##min.t过滤掉在细胞中出现频率低的ACRs
##max.t过滤掉在细胞中出现频率高的ACRs，默认为前0.005
soc.obj=cleanData(soc.obj,min.c=1000,min.t=0.005,max.t=0.005)
#Normalization
##
soc.obj=regModel(soc.obj,nthreads=16)
tfidf.obj=tfidf(soc.obj)
lr.obj=logisticModel(soc.obj,nthreads = 16)
regMod2.obj=refModel2(obj)
#SVD降维
##method="SVD"指定使用SVD降维，也可选择method="NMF"
soc.obj=reduceDims(soc.obj, n.pcs=50, cor.max=0.7, verbose=T)
soc.obj=reduceDims(soc.obj, n.pcs=50, cor.max=0.7, method="NMF")
#UMAP
soc.obj=projectUMAP(soc.obj, verbose=T)
plot(soc.obj$UMAP, pch=16, cex=0.2)
#Graph-based clustering
soc.obj=callClusters(soc.obj, res=0.6, verbose=T)
plotUMAP(soc.obj)
saveRDS(soc.obj, file="Socrates_object.rds")
