search

quanc1989 / SV-ONT-Tibetan

Characterization of Structural Variation in Chinese samples
MIT License
15 stars 7 forks source link

使用paragraph对sniffles的结果进行genotype #3

Closed Yechao94 closed 2 years ago

Yechao94 commented 2 years ago

您好， 我看到您再文章里面使用paragraph对sniffles的结果进行了genotyping，我也尝试使用使用这个piepline进行计算，但是发现在这个过程中，总是会出现报错，想请教一下，我需要对vcf文件进行怎么样的预先处理吗？

感谢~

quanc1989 commented 2 years ago

您好， 我看了一下，我做的预处理的步骤最主要的就是删除TRA，然后将SV在vcf中的indel序列转成了Symbol，增加INS的seq序列，确定DEL的起始点、长度、和原来indel格式中的seq是否一致，不一致就按照参考基因组的序列。

你可以把错误发到paragraph的github中的issues去，也可以发给我看看。

希望对你有帮助。

Yechao94 commented 2 years ago

感谢您的回答。 对于TRA和INV我进行了删除操作。 这是我的部分log文件以及部分变异位点~ log：

2021-12-20 16:42:19,122 ERROR      File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/lib/python3/grm/vcf2paragraph/__init__.py", line 286, in run_vcf2paragraph    alt_paths=params["alt_paths"])
2021-12-20 16:42:19,122 ERROR      File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/lib/python3/grm/vcf2paragraph/__init__.py", line 86, in convert_vcf    ref, indexed_vcf.name, ins_info_key, chrom, start, end, ref_node_padding, allele_graph)
2021-12-20 16:42:19,122 ERROR      File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/lib/python3/grm/vcfgraph/vcfgraph.py", line 128, in create_from_vcf    graph.add_record(record, allele_graph, varId, ins_info_key)
2021-12-20 16:42:19,122 ERROR      File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/lib/python3/grm/vcfgraph/vcfgraph.py", line 178, in add_record    f"Missing key {ins_info_key} for <INS> at {self.chrom}:{vcf.pos}; ")
2021-12-20 16:42:19,122 ERROR    Exception: Missing key SEQ for <INS> at 12:86259792; 
2021-12-20 16:42:19,475 WARNING  1:219904230 Padding base in genome is different from VCF. Use the one from genome.
The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/bin/multigrmpy.py", line 353, in <module>
    main()
  File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/bin/multigrmpy.py", line 349, in main
    run(args)
  File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/bin/multigrmpy.py", line 261, in run
    graph_files = load_graph_description(args)
  File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/bin/multigrmpy.py", line 52, in load_graph_description
    header, records, event_list = convert_vcf_to_json(args, alt_paths=True)
  File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/lib/python3/grm/vcf2paragraph/__init__.py", line 156, in convert_vcf_to_json
    variants = pool.map(run_vcf2paragraph, zip(to_process, itertools.repeat(params)))
  File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/lib/python3.7/multiprocessing/pool.py", line 268, in map
    return self._map_async(func, iterable, mapstar, chunksize).get()
  File "/public/home/fan_lab/ychhuang/miniconda3/envs/genotype/lib/python3.7/multiprocessing/pool.py", line 657, in get
    raise self._value
Exception: Missing key SEQ for <INS> at 12:86259792;

SVs

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  final.sorted.bam
1       66224   0       TATTATATAATATATATTATATTATATAATATATAATATAAATATAATATAAATT T       .       PASS    PRECISE;SVMETHOD=Snifflesv1.0.12;CHR2=1;END=66278;ZMW=16;STD_quant_start=1.264911;STD_quant_stop=0.707107;Kurtosis_quant_start=-0.766439;Kurtosis_quant_stop=-0.670589;SVTYPE=DEL;RNAMES=m161114_050501_42177R_c101134222550000001823248606131710_s1_p0/20519/8698_36789,m161226_095604_42156_c101148962550000001823276307071774_s1_p0/9578/0_12422,m161229_035312_42177R_c101137192550000001823250106131745_s1_p0/161232/0_16589,m170102_141109_42177R_c101134132550000001823248606131736_s1_p0/101325/0_16030,m170104_142136_42177R_c101158452550000001823249707191772_s1_p0/7884/0_17669,m170116_022514_42177R_c101159392550000001823249707191756_s1_p0/141636/0_22229,m170118_060525_42156_c101134752550000001823248606131730_s1_p0/157517/19758_34195,m170118_060525_42156_c101134752550000001823248606131730_s1_p0/157517/34246_38062,m170119_131124_42163R_c101160152550000001823271007071782_s1_p0/141707/362_19489,m170120_105240_42156_c101160782550000001823271007071771_s1_p0/3964/0_21927,m170123_133418_42177R_c101159262550000001823249707191716_s1_p0/97583/27559_41997,m170315_224025_42156_c101178502550000001823273908211743_s1_p0/134752/0_5762,m170315_224025_42156_c101178502550000001823273908211743_s1_p0/44779/0_21459,m170315_224025_42156_c101178502550000001823273908211743_s1_p0/44779/21502_32355,m170320_060736_42177R_c101178492550000001823273908211781_s1_p0/73809/0_14712,m170322_024920_42163R_c101187242550000001823278408081785_s1_p0/15797/0_22657,m170330_065326_42163R_c101170682550000001823273408151751_s1_p0/52704/0_24907,m170331_150548_42156_c101201392550000001823280410241774_s1_p0/163291/1855_6764;SUPTYPE=AL;SVLEN=-54;STRANDS=+-;STRANDS2=10,8,10,8;RE=18;REF_strand=35,42;Strandbias_pval=0.600978;AF=0.233766;IRIS_PROCESSED=1;IRIS_REFINED=0    GT:DR:DV        0/0:59:18

1       90367   2       C       CTGGAGGAAGACAGTTCCCTCTTCCTCTGTTCTGCCAACCAGTTACTGCTGCTCTGGAGAAGACAGTCCATGTCCCTCTGTCTCTGCAACCAGCTTTTCTGCTCTTT     .       PASS    IMPRECISE
;SVMETHOD=Snifflesv1.0.12;CHR2=1;END=90367;ZMW=20;STD_quant_start=29.637200;STD_quant_stop=29.232609;Kurtosis_quant_start=-0.796489;Kurtosis_quant_stop=-0.720155;SVTYPE=
INS;RNAMES=m161106_080829_42163R_c101115672550000001823260305221756_s1_p0/79599/0_15594,m161107_093215_42163R_c101038532550000001823258912031682_s1_p0/150591/0_7822,m161
107_093215_42163R_c101038532550000001823258912031682_s1_p0/4696/0_22262,m161115_170044_42163R_c101134722550000001823248606131764_s1_p0/155803/0_11365,m161226_133708_4216
3R_c101134022550000001823248606131774_s1_p0/248/0_20636,m161227_010503_42156_c101148962550000001823276307071777_s1_p0/50560/5711_22100,m161227_072656_42156_c101134802550
000001823248606131753_s1_p0/159092/0_13309,m161230_045914_42156_c101134792550000001823248606131791_s1_p0/84535/0_8651,m170106_075143_42156_c10113468255000000182324860613
1735_s1_p0/139044/23504_38538,m170106_075143_42156_c101134682550000001823248606131735_s1_p0/139044/38585_52021,m170106_075143_42156_c101134682550000001823248606131735_s1
_p0/139044/52067_67023,m170113_123611_42177R_c101134762550000001823248606131723_s1_p0/42393/4120_18619,m170113_215244_42156_c101159192550000001823249707191717_s1_p0/9644
2/1562_16851,m170115_132920_42156_c101160652550000001823271007071736_s1_p0/102011/1808_15785,m170115_132920_42156_c101160652550000001823271007071736_s1_p0/93899/0_14041,
m170116_112830_42156_c101158082550000001823249707191762_s1_p0/55390/2497_24863,m170117_153826_42163R_c101159202550000001823249707191772_s1_p0/36774/0_19243,m170120_06231
8_42163R_c101160152550000001823271007071785_s1_p0/154245/0_15990,m170320_060736_42177R_c101178492550000001823273908211781_s1_p0/119076/4738_15973,m170321_202016_42163R_c
101187242550000001823278408081784_s1_p0/53905/0_13195,m170321_202016_42163R_c101187242550000001823278408081784_s1_p0/61358/272_16503,m170331_084902_42163R_c1011706825500
00001823273408151755_s1_p0/146191/12433_26593;SUPTYPE=AL;SVLEN=106;STRANDS=+-;STRANDS2=14,8,14,8;RE=22;REF_strand=39,37;Strandbias_pval=0.341023;AF=0.289474;IRIS_PROCESS
ED=1;IRIS_REFINED=0        GT:DR:DV        0/0:54:22

再次感谢

quanc1989 commented 2 years ago

虽然paragraph写的是在symbol格式才需要在INFO中加入SEQ属性，但是你这个错误不是明显写了需要一个SEQ的键值对么？ 你可以把序列加入info中试试，应该就可以了

Yechao94 commented 2 years ago

好的 我试试。 万分感谢