我以前使用 GATK 和 vcftools 来调用我的 scRNA-seq 数据中的变体。现在我正在尝试使用 bcftools (v.1.9) 来查看是否得到相同的变体。
到目前为止，我已经完成了以下工作：
for fn in $(find /folder/ -name "*.bam")
do
bcftools mpileup -O u -o /folder/"$(basename $fn .bam)".bcf -f /home/cluster/colpe/scratch/ref_sequences/Mus_musculus.GRCm38.dna.primary_assembly.fa "$fn"
done
这创建了一个我用 htsfile 检查的 bcf 文件，它看起来很好。然后我运行 bcftools 调用：
for fn in $(find /folder/ -name "*.bam")
do
bcftools call -m -v -o /folder/"$(basename $fn .bcf)".vcf "$fn"
done
这给了我一个vcf文件。再次使用 htsfile 检查它给了我我所期望的。简要检查它具有所有预期列的文件：（以下是一个示例行）#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 20190219_Gli1_5d_A03 1 6226714。TT 37.4152。INDEL;IDV=2;IMF=1;DP=2;VDB=0.1;SGB=-0.453602;MQ0F=0;AC=2;AN=2;DP4=0,0,2,0;MQ=60 GT： PL 1/1:67,6,0
现在我想根据读取深度进行过滤。我只想要 DP=10 或更多的职位。我已经尝试过 vcftools (v0.1.16)，但它给了我一个空文件，即使我知道有 DP>10 的位置。这是我运行的 vcftools 代码：
for fn in $(find /folder/ -name *.vcf)
do
vcftools --vcf "$fn" --out "$fn"_dp10 --min-meanDP 10 --recode --recode-INFO-all
done
然后我用这段代码尝试了 bcftools 视图：
for fn in $(find /folder/ -name *.vcf)
do
bcftools view -O v -i INFO/DP>9 & TYPE="snp" -o /folder/"$(basename $fn .vcf)".vcf "$fn"
done
但是，这给了我错误：
/var/lib/slurm/slurmd/job3764707/slurm_script: line 11: -o: command not found
Failed to open -: unknown file type
我已经尝试了几个小时，但在任何论坛中都找不到任何答案。非常感谢一些指导！
提前感谢科拉