1、计算1)全基因组平均测序深度 2）不同深度（1X 4X 10X 20X）下的测序覆盖度

for i in `ls -d SRR*`
do
        {
                bedtools genomecov -ibam $i/$i'.sorted.bam'|grep genome>$i/$i'.bedtools.genomecov.res'
        }&
done
#接R的统计代码
info.list<-data.frame(matrix(0,nrow =length(grep("SRR*",list.files())) ,ncol =6 ))
n=0
for (i in list.files()[grep("SRR*",list.files())]) {
  n=n+1
  dir<-paste(i,"/",i,".bedtools.genomecov.res",sep = "")
  x<-read.table(dir,sep = "\t")
  x<-x[grep("genome",x$V1),]
  info.v<-vector()
  info.v[1]<-i
  info.v[2]<-sum(x$V3[-1]*x$V2[-1])/x$V4[1]
  info.v[3]<-sum(x$V3[-1])/x$V4[1]
  info.v[4]<-sum(x$V3[-c(1:4)])/x$V4[1]
  info.v[5]<-sum(x$V3[-c(1:10)])/x$V4[1]
  info.v[6]<-sum(x$V3[-c(1:20)])/x$V4[1]
  info.list[n,]<-info.v
}
colnames(info.list)<-c("ID","mean_depth","1X","4X","10X","20X")

image.png

2、卡合适的阈值(质量值至少为30，深度至少为4，以及其它的硬阈值)的比对结果纳入高质量SNP范围（这里使用的gatk工具）

#!/bin/bash
for id in `ls -d SRR*`
do
        {
                input=$id/$id'_raw.vcf'
                output=$id/$id'_snp.vcf'
                gatk SelectVariants -select-type SNP -V $input -O $output
        }&
done

for id in `ls -d SRR*`
do
        {
                input=$id/$id'_raw.vcf'
                output=$id/$id'_indel.vcf';
                gatk SelectVariants -select-type INDEL -V $input -O $output
        }&
done
wait



for id in `ls -d SRR*`
do
        {
        input=$id/$id'_snp.vcf'
        output=$id/$id'_snp.filtered.vcf'
        gatk VariantFiltration -V $input --filter-expression "DP<4 ||  QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "Filter" -O $output
        }&
done 1>hard_filter_snp.log 2>hard_filter_snp.errlog


for id in `ls -d SRR*`
do
        {
                input=$id/$id'_indel.vcf'
                output=$id/$id'_indel.filtered.vcf'
                gatk VariantFiltration -V $input --filter-expression "DP<4 ||  QD < 2.0 ||  FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "Filter" -O $output
        }&

done 1>hard_filter_indel.log 2>hard_filter_indel.errlog

wait

for id in `ls -d SRR*`
do
        {
                input1=$id/$id'_snp.filtered.vcf'
                input2=$id/$id'_indel.filter.vcf'
                output=$id/$id'.filtered.vcf'
                gatk MergeVcfs -I $input1 -I $input2 -O $output
        }&
done 1>merge_snp_indel.log 2>merge_snp_indel.errlog

#这里暂未考虑质量值的问题！！！

3、四组数据的比对结果中，公共区域的SNP的交并补情况

因为四组数据是单细胞数据，本身测序深度为1X左右，gap很多，因此希望将四个样本共同测到的区域用于做后续的统计

#1、找到四个样本的公共区域
#方案：
#1）将参考基因组用小滑窗划分为小区间
#2）使用bedtools计算四个样本在step1文件中的平均覆盖度
#3）挑选出覆盖度>0的区间认为是单个样本在全基因组上的覆盖区域！
#4）将四个样本的step3的结果文件取交集，认为是四个样本的公共区域！
#缺点：存在误差，只是以平均覆盖度为参考依据。
#bamtobed
for id in `ls -d SRR*`
do
        {
                bedtools bamtobed -i $id/$id'.sorted.bam' >$id/$id'.sorted.bam.to.bed'
        }&
done
#做-g文件
samtools dict index/hg38.fa >hg38.dict
cat hg38.dict|grep -v '^@HD'|sed 's/:/\t/g'|cut -f 3,5 >genome.chr.ln
less genome.chr.ln |sed -n '/chr[0-9XYM]\{1,2\}\t/p' >genome.chr.ln.new#只保留chr1-22 X Y M
bedtools makewindows -g genome.chr.ln.new -w 2000 >2000.genome#两千的窗口
bedtools makewindows -g genome.chr.ln.new -w 50 >50.genome#五十的窗口
#将四个样本分别于-g文件做比较，统计overlap的频率，我们把overlap>0的窗口认为该样本在此窗口整体有read覆盖！
ls -d fastq/SRR*|while read id ;do  input=$id/${id##fastq/}'.sorted.bam.to.bed';bedtools intersect -a 2000.genome -b $input -c >$id/${id##fastq/}'.2000windows.hg38.intersect.bed';done &
ls -d fastq/SRR*|while read id ;do  input=$id/${id##fastq/}'.2000windows.hg38.intersect.bed';less $input|awk '$4>0{print $0}'  >$input'.filter.bed';done
#计算四个样本的公共区域
bedtools intersect -a fastq/SRR630578/SRR630578.2000windows.hg38.intersect.bed.filter.bed -b fastq/SRR630579/SRR630579.2000windows.hg38.intersect.bed.filter.bed >sub1.bed&
 bedtools intersect -a fastq/SRR630580/SRR630580.2000windows.hg38.intersect.bed.filter.bed -b fastq/SRR630581/SRR630581.2000windows.hg38.intersect.bed.filter.bed >sub2.bed&
bedtools intersect -a sub1.bed -b sub2.bed >sub3.bed &
#则sub3.bed则为四个样本的公共区域

#2、将四个样本的SNP和indel与公共区域取交集
ls -d fastq/SRR*|while read id ;do input2=sub3.bed;input1=$id/${id##fastq/}'_snp.filtered.vcf';output=$id/${id##fastq/}'_snp.filtered.intersect.vcf';bedtools intersect -a $input1 -b $input2 -wa >$output;done

ls -d fastq/SRR*|while read id ;do input2=sub3.bed;input1=$id/${id##fastq/}'_indel.filtered.vcf';output=$id/${id##fastq/}'_indel.filtered.intersect.vcf';bedtools intersect -a $input1 -b $input2 -wa >$output;done


#3、计算step2输出的vcf的交并补
#接R语言
library(stringr)
library(UpSetR)
fq_dir<-"../../../fastq/"
id<-c("SRR630578","SRR630579","SRR630580","SRR630581")
vcf_list<-lapply(id,function(x){
  file_dir<-paste0(fq_dir,"/",x,"/",x,"_snp.filtered.intersect.vcf")
  vcf<-read.table(file_dir)
  vcf<-vcf[,c(1,2,3,4,5,7)]
  vcf<-vcf[vcf$V7=="PASS",]
  vcf$pos<-paste0(vcf$V1,vcf$V2)
  vcf$alter<-paste0(vcf$pos,vcf$V5)
  return(vcf)
  
})

#做突变位置的统计
upset(fromList(list(a1=vcf_list[[1]]$pos,
                    a2=vcf_list[[2]]$pos,
                    a3=vcf_list[[3]]$pos,
                    a4=vcf_list[[4]]$pos)))
#做突变位置+突变碱基的交并补统计
upset(fromList(list(a1=vcf_list[[1]]$alter,
                    a2=vcf_list[[2]]$alter,
                    a3=vcf_list[[3]]$alter,
                    a4=vcf_list[[4]]$alter)))