低深度测序下统计各个细胞的共同测到的SNV情况

前言：在低深度（小于1X，覆盖度低于50%）下进行单个细胞的全基因组测序，与高深度（30X或60X）的近乎100%覆盖率相比，如果想得到每个细胞的SNV的交并补情况是比较困难的，因为在vcf中，如果某个位点没有被记录，可能有两种情况：1、这个位点没有发生突变 2、这个位点未被测到。因此我们需要：1、统计出每个细胞的数据中共同测到的部分，2、该部分中进行交并补统计

注意：这里仅限于统计细胞间的交并补情况，如果是对单个细胞的统计，不应该做第一步！

一、统计每个细胞的数据中共同测到的部分

#!/bin/bash
#得到不同样本的碱基覆盖情况，在这里我将最低深度调整为4
nohup time ls -d FDSW20259055*|while read id ;do cd $id;samtools depth  $id.sorted.bam|awk '$3>4{print $0}' >../all_chr_stat/$id.all_chr.depth ;cd ../;done &
#求样本覆盖的交集  此处为R环境
library(stringr)
library(UpSetR)
file_pattern="depth$"
sample.name<-list.files()[grep(".*depth$",list.files())]
sample.name<-str_split(sample.name,"\\.",simplify = T)[,1]
mat_list<-lapply(grep(file_pattern,list.files()),function(x){
                           a=read.table(list.files()[x],sep = "\t")
                             a[,4]=paste(a[,1],a[,2],sep = "_")
                             return(a)
})
mat_pos_list=lapply(mat_list,function(x){x[,4]})
all_intersect=Reduce(intersect,mat_pos_list)
all_intersect_pos=all_intersect
all_intersect=data.frame(chr=str_split(all_intersect,"_",simplify = T)[,1],
                                                  pos=str_split(all_intersect,"_",simplify = T)[,2])

write.table(all_intersect,file = "intersect.res.txt",quote = F,row.names = F)

二、求每个vcf文件在区域的交集

#接第一步的R环境
vcf_list<-lapply(grep("vcf$",list.files()),function(x){
                           a=read.table(list.files()[x],sep = "\t")
                             a$pos<-paste(a$V1,a$V2,sep = "_")
                             a=a[a$pos%in%all_intersect_pos,]
                                                  })

vcf_pos_list<-lapply(vcf_list,function(x){
                               x$pos
                                                  })
names(vcf_pos_list)<-sample.name
pdf("interset_plot.pdf")
upset(fromList(vcf_pos_list))
dev.off()

缺点：R读大文件，尤其是samtools计算出来的depth文件，太吃内存了，有没有好的解决办法？？？