Socrates是一款用于分析snATAC數(shù)據(jù)的R語言包嘿悬，但是該包發(fā)表后并沒有引起多大水花实柠，最直觀的就是，Github上該包的Issues只有10條善涨，其中兩條還是本人最近貢獻(xiàn)的窒盐。這里附上Socrates的Github鏈接：https://github.com/plantformatics/Socrates茶行。
雖然這個包沒啥人用，但是由于本人愛折騰登钥，最終還是嘗試了一把，接下來講一些使用該包的技巧娶靡，絕對是全網(wǎng)首發(fā)Ｄ晾巍！姿锭！
第一步：讀取輸入文件
輸入文件總共需要3個塔鳍，分別是Tn5結(jié)合位點(diǎn)文件，基因組的gff注釋文件呻此，基因組染色體長度文件轮纫。
但是通常情況下，我們上游通過cellranger-atac得到的是fragment文件焚鲜，這是需要在讀取時設(shè)置is.fragment=T掌唾。
library(Socrates)
#讀取文件
bed="SampleName.fragments.tsv.gz"
ann="genome.gff3"
chr="genome.chr.bed"
#建立obj，注意Socrates包默認(rèn)的bed為包含正負(fù)鏈信息的單堿基Tn5插入位點(diǎn)忿磅，但我們的fragment文件是從cellranger-atac中得到的糯彬，因此在讀取時需要設(shè)置is.fragment=T
obj=loadBEDandGenomeData(bed,ann,chr,is.fragment=T)
這里仍然存在一個小問題就是如果fragment文件過大的話，是無法讀取的葱她，這里需要自行將fragment文件轉(zhuǎn)換成單堿基的Tn5插入位點(diǎn)文件撩扒，直接進(jìn)行讀取。至于轉(zhuǎn)換的方式吨些，則是根據(jù)cellranger-atac產(chǎn)生fragment文件的方式來逆向還原Tn5插入位點(diǎn)的位置搓谆。
cellranger-atac count的結(jié)果中會生成fragments.gz文件，其中每一個fragment是由兩個獨(dú)立的轉(zhuǎn)座事件創(chuàng)造的豪墅，也就是說fragment的兩端正好是兩次Tn5酶插入的位點(diǎn)泉手。每一個特異的fragment會產(chǎn)生多條冗余的reads，這些冗余的reads被折疊成單個fragment記錄偶器。
Fragment的位置信息是從bam文件中獲得的螃诅。fragment的起點(diǎn)是將bam中的reads比對的位置最左側(cè)向前移動4bp，終點(diǎn)是將reads比對位置的最右側(cè)向后移動5bp状囱。
因此將fragment文件中每條fragment的起點(diǎn)和終點(diǎn)都分別看做一個Tn5結(jié)合位點(diǎn)即可术裸。
將fragment文件轉(zhuǎn)換成Tn5結(jié)合位點(diǎn)文件
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" &&
第二步：進(jìn)行peak calling并過濾低質(zhì)量的細(xì)胞
#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袭艺，否則不會計(jì)算frip，在后續(xù)的過濾過程中會出問題
obj=buildMetaData(obj,tss.window=3000)
head(obj$meta)
#過濾低質(zhì)量的細(xì)胞
#findCells后會生成多個meta叨粘，分別為meta.v1猾编，meta.v2瘤睹，meta.v3，他們分別對應(yīng)pdf圖中的三次過濾答倡，其中meta.v3為最終的過濾結(jié)果轰传。
##min.cells和max.cells代表最低和最高的細(xì)胞數(shù)目
##min.tn5代表Tn5位點(diǎn)處最低的read depth
##filt.tss=T代表根據(jù)tss進(jìn)行過濾
##tss.min.freq=0.2和tss.z.thresh代表reads比對到TSS位點(diǎn)處的比例超過0.2，同時z-score的閾值為3
##filt.frip=T代表根據(jù)frip進(jìn)行過濾
##frip.min.freq=0.2代表細(xì)胞的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)
這里在經(jīng)過質(zhì)控過濾后瘪撇，Socrates并不會丟棄那些被過濾掉的細(xì)胞获茬，而是增加多個meta來進(jìn)行區(qū)分。
第三步：建立Socrates對象并進(jìn)行下游分析
#生成Socrates對象
##filtered=T倔既，指定用過濾后的細(xì)胞建立matrix恕曲，windows=1000指定基因組劃窗的大小為500，peaks=F不適用peak文件進(jìn)行定量渤涌，如果設(shè)置peaks=T佩谣，則對peak進(jìn)行定量，并覆蓋windows=500的參數(shù)实蓬。
obj=generateMatrix(obj,filtered=T,windows=500,peaks=F,verbose=T)
obj=generateMatrix(obj,filtered=T,peaks=T,verbose=T)
#將該對象轉(zhuǎn)換為Socrates格式
soc.obj=convertSparseData(obj)
saveRDS(obj,file="QC_object.rds")
saveRDS(soc.obj,file="Socrates_object.rds")
#對Socrates對象進(jìn)行過濾
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對象進(jìn)行過濾
##min.c設(shè)置每個細(xì)胞中accessible ACRs的最小值
##min.t過濾掉在細(xì)胞中出現(xiàn)頻率低的ACRs
##max.t過濾掉在細(xì)胞中出現(xiàn)頻率高的ACRs茸俭，默認(rèn)為前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")
