starskyzheng / panpop

Application of pan-genome for population
MIT License
98 stars 9 forks source link

我想知道Fill_DP.pl的用处 #48

Closed hanshenglei closed 6 months ago

hanshenglei commented 6 months ago

你好,郑博士! 我有一个群体的sv文件,我想知道Fill_DP.pl在PART_run.pl之前运行还是在它之后运行,以及RADEME中你提到Fill_DP.pl仅添加深度信息,那该如何做后续处理,比如文中提到的软屏蔽;还是PART_run.pl运行完了之后再运行Fill_DP.pl就过滤完了 期待您的回复 祝好!

hanshenglei commented 6 months ago

看了您的文章之后,我想知道过滤深度信息和过滤缺失率的脚本分别是什么 是panpop/scripts/flt_raw_vcfs.pl和panpop/scripts/flt_vcf_by_miss.pl吗

starskyzheng commented 6 months ago
  1. 先Fill_DP.pl,再PART_run.pl
  2. 对的。具体参数可以看snakefile里写的,比如: https://github.com/starskyzheng/panpop/blob/df6d0e72151d0456bf1bcccff6ad00ad27169f6f/subworkflows/panpopbase.py#L124
hanshenglei commented 6 months ago

感谢您的回复!

hanshenglei commented 6 months ago

1) 2.callSV/s1/s1-Ref.gaffe.aug.q5.call.ext.vcf.gz 2.callSV/s2/s2-Ref.gaffe.aug.q5.call.ext.vcf.gz perl /home/hsl/workdir/panpop/scripts/merge_vcf.pl --inlist 3.merge_rawvcf/1.inputvcfs.list --out 3.merge_rawvcf/2.merge_rawvcf.1.vcf.gz --tmp_dir tmp --vcfs_per_run 10000 --threads 11 --bcftools_threads 4 -r 1 --m_none 0 >> logs/3.2.merge_rawvcf.1.log 2>&1

6.aug_dp/3.vcf_with_dp.1/s1.sorted.vcf.gz 6.aug_dp/3.vcf_with_dp.1/s2.sorted.vcf.gz perl /home/hsl/workdir/panpop/scripts/merge_vcf.pl --inlist 7.aug_merge_rawvcf/1.inputvcfs.1.list --out 7.aug_merge_rawvcf/2.merge_rawvcf.1.vcf.gz --tmp_dir tmp --vcfs_per_run 10000 --threads 11 --bcftools_threads 4 --m_none 1 >> logs/7.2.merge_rawvcf.1.log 2>&1

perl /home/hsl/workdir/panpop/scripts/merge_vcf.pl -h usage: perl xx.pl --inlist FILE.vcfs.list --out FILE.vcf.gz --tmp_dir PATH [--chrome STRING] Chromosome name [--threads INT (default: 1)] Number of bcftools running together [--vcfs_per_run INT (default: 999999)] [--bcftools_threads INT (default: 4)] Number of threads of 'bcftools merge' [--noclean] Do not clean tmp files [--invcf_nocheck] Do not check input vcf.gz files wather they are exists and has .tbi files. [--help] Print this help message -r 和 -m_none 您好像没有解释清楚? 2)merge_vcf_same_pos.pl这个脚本是用来解决以下情况吗:如果两个VCF文件中存在相同位置的变异记录,需要进行合并。这通常涉及到解决冲突,例如选择保留哪一个变异记录或者合并它们。 3)为什么flt_raw_vcfs.pl在example_short_reads/3.merge_rawvcf中4.filter_raw_vcf.1.vcf.gz使用了一次,此时还没有添加深度信息好像;6.aug_dp/之后才带有深度信息,并在7.aug_merge_rawvcf使用了flt_raw_vcfs.pl

hanshenglei commented 6 months ago

我观察到 GetOptions ( 'help|h!' => \$opt_help, 'inlist|i=s' => \$inlist, 'out|o=s' => \$out, 'threads|t=i' => \$threads, 'tmp_dir|T=s' => \$tmp_dir_top, 'vcfs_per_run|V=i' => \$vcfs_per_run, 'invcf_nocheck!' => \$invcf_nocheck, 'noclean!' => \$noclean, 'bcftools_threads=i' => \$bcftools_threads, 'chrome|r=s' => \$chrome, 'bcftools_bin=s' => \$bcftools, 'm_none=i' => \$m_none, 'verb|v!' => \$verb, ); -r 应该是chrome 默认应该是所有 -m_none=1 为bcftools -m none 是这样吗

hanshenglei commented 6 months ago

以上都不是最重要的,我现在遇到的主要问题为怎么更改REF ALT GT;如您所说https://github.com/starskyzheng/panpop/issues/18%EF%BC%9B%E4%BB%8E%E4%B8%8A%E4%B8%8B%E6%96%87%E4%B8%AD%E6%88%91%E5%BE%97%E7%9F%A5%E5%AF%B9%E4%BA%8Esv 的处理会比cnv更合理,请问怎么将不同软件的sv结果转变成PART需要的形式。 我的流程应该是得到manta delly二代重测序call出来的sv文件,对文件进行更改(REF ALT),进行bcftools merge -m none,然后去除一些same_position,进行FILL_DP;然后在进行PAPT,最后进行split sv得到最终的结果,请问我从您的文章和所有的问答当中得到的理解对吗

hanshenglei commented 6 months ago

我已经脚本改好,可以转换delly和manta的结果了,请回答除此之外的问题, 4)后续分析的时候,比如跟表型关联,也可以 沿用这种格式吗? 5)vcf本来就包含深度信息,那么FILL_DP是填充的跟什么有关的深度信息,比如vcf是跟bam和基因组相关 6)merge_vcf_same_pos.pl 中--skip_mut_at_same_pos 您使用的值是2 请问这个是什么意思 是采用这个参数吗 7)long_inv_bed2vcf.pl得到的00.inv.vcf.gz之后再与thin之后的vcf进行concat是为什么

starskyzheng commented 6 months ago

-r 和 -m_none 您好像没有解释清楚?

我加了一些说明 -r就是染色体 --m_none 就是跑bcftools merge时候的参数 不建议单独调用这个脚本 07ab3f15c44f70bc9e0b077f711229a9557f4cdc

2)merge_vcf_same_pos.pl这个脚本是用来解决以下情况吗:如果两个VCF文件中存在相同位置的变异记录,需要进行合并。这通常涉及到解决冲突,例如选择保留哪一个变异记录或者合并它们。

不存在这种情况。这一步只是粗合并,只有当CHR、POS、ALT、REF都一样时候才会合并。

3)为什么flt_raw_vcfs.pl在example_short_reads/3.merge_rawvcf中4.filter_raw_vcf.1.vcf.gz使用了一次,此时还没有添加深度信息好像;6.aug_dp/之后才带有深度信息,并在7.aug_merge_rawvcf使用了flt_raw_vcfs.pl

将不准的位点删掉,可以减轻补全深度时候的计算量

4) 这是标准vcf格式,后续分析都可用。但是建议只保留二等位基因的SV

5)原本的VCF里只有突变个体的深度信息,补全的是参考基因型的个体的深度信息

6)见: https://github.com/starskyzheng/panpop/blob/07ab3f15c44f70bc9e0b077f711229a9557f4cdc/scripts/realign.pl#L637

7)panpop的这种方法处理inv容易不准确,所以inv是单独进行计算的。最后再合并

hanshenglei commented 6 months ago

感谢郑博士的回复!

1)您好,我在使用merge_vcf_same_pos.pl之后,进行bcftools sort的时候报错:[W::vcf_parse_format] FORMAT 'DP' at LG1:1746 is not defined in the header, assuming Type=String Error encountered while parsing the input at LG1:1746 而且是生成merge_vcf_same_pos.pl的结果文件中,GT:DP才出现,之前并不是DP;运行时的log为[W::bcf_hdr_merge] Trying to combine "SVLEN" tag definitions of different lengths [W::bcf_hdr_merge] Trying to combine "HOMLEN" tag definitions of different lengths 2)我发现将两个软件call出来的vcf中的INV进行合并的时候,运行scripts/long_inv_bed2vcf.pl,得到的vcf为空,也就是说INV全被过滤了吗

starskyzheng commented 6 months ago

1) DP信息在vcf头里,没有的话手动补上吧 ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">

2) inv的分离是prase那步进行的。不用管了,inv反正也很少。

hanshenglei commented 6 months ago

您好,郑博士 1)我看到您的example_long_reads中不同样本中的10.concat_inv.vcf.gz,也就是INV与PART之后的vcf进行concat之后的vcf文件--包含INV的vcf,进行过split之后就没有再使用了,在下一步整合不同样本的sv然后填充深度信息PART的时候,使用的vcf文件都为之前每个样本剔除INV的vcf文件,所有说INV是不需要分析吗还是panpop不适合分析INV

hanshenglei commented 6 months ago

2)我还发现您在example_short_reads中使用flt_raw_vcfs.pl过滤添加完深度信息的vcf时没有加 --flt_append_only参数,但是在example_short_reads中过滤所有样本的添加完深度信息的vcf时使用了这个参数,请问这是为什么?是因为之前已经过滤了一遍,所以这一次为了加快计算只过滤添加的吗?

starskyzheng commented 6 months ago

1) INV建议作bed单独合并,用PART处理很可能不准确 2) 不需要了,因为所有位点都有深度信息了。

hanshenglei commented 6 months ago

感谢您的回复!

starskyzheng commented 6 months ago

不客气!没有更多问题我就关闭这个issue了。如果有更多问题请再新开一个issue