问题描述
我之前曾使用 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 (将#修改为@)