如何使用 bcftools 视图根据深度过滤我的 vcf 文件?

问题描述

我之前曾使用 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

这创建了一个 bcf 文件,我用 htsfile 检查了它,看起来不错。 然后我运行 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 过滤器信息格式 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,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

我已经尝试了几个小时,但在任何论坛中都找不到任何答案。非常感谢您的指导!

提前致谢 科拉

解决方法

暂无找到可以解决该程序问题的有效方法,小编努力寻找整理中!

如果你已经找到好的解决方法,欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@)