search

twolinin / longphase

GNU General Public License v3.0
102 stars 9 forks source link

modcall resulting vcf not phased #89

Closed sahuno closed 1 week ago

sahuno commented 3 weeks ago

Hi, i run modcall with modified.bams that have MN & ML tags however the resulting don't have | at the genotype position meaning the results were unphased. here's command that i used for modcall and to check phasing. Que: is the output not supposed to be phased? pls let me know if i'm missing any step

longphase modcall -b /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/results/mark_duplicates/D-0-1_4000/D-0-1_4000_modBaseCalls_sorted_dup.bam -t 24 -o modcall_D-0-1_4000.vcf -r /data1/greenbab/database/mm10/mm10.fa
# grep -v "^#" modcall_D-0-1_4000.vcf.vcf | awk '{if ($10 ~ /\|/) print "Phased genotype found at position: " $2}'

pls see snapshot of vcf 

##longphaseVersion=1.7.3
##commandline="modcall -b /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/results/mark_duplicates/D-0-1_4000/D-0-1_4000_modBaseCalls_sorted_dup.bam -t 24 -o modcall_D-0-1_4000.vcf -r /data1/greenbab/datab
ase/mm10/mm10.fa "
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  SAMPLE
chr1    3348011 .       c       N       .       PASS    RS=P;MR=acb15622-fbae-4dac-b3ef-840801870adc,4cfa7809-6b0c-41e1-85ef-4783c4d90aa4,762389c0-9de2-4f87-b401-0e1a43ea3a79,08b39281-c4c9-42bd-86bf-257620cc8baf,001dcd7d-f2c0-4e
7a-8780-0e45bf654d98,454403f1-a91d-47b5-a1a2-4847f87aead0;NR=05dfa142-5a31-4e2f-a403-b7bb63017758,5e70609a-8118-40e3-a797-1edfb0946f03,3b32029f-ab8f-495f-8da6-40409609a159,601c67dd-b463-42dd-ba9b-01205ea933fc,113f4d42-4a76-44e2-
8b26-1d3c0a72632c;      GT:MD:UD:DP     0/1:6:5:13
chr1    3348012 .       g       N       .       PASS    RS=N;MR=0f4d1ff8-0676-42e6-9262-b7f5a06a3f1b,cbb6aa20-df72-434f-9a4f-fc39740cb467,2407bdfb-cac6-4f15-b884-e37e2cb63921,60a1c1a5-3c8e-48e3-b8c7-48151cbefacf;NR=83143f28-df43
-42d0-9e20-aa25f4a88e82,a6bf78fa-cf63-44b1-a980-f471ae3f6501,1f0caf97-8992-47c0-9e9a-0c1f7167b8b9,c828e195-19c8-4eb8-a43e-7066d621504c,14ca456e-d0cd-492b-8ebd-b3abd84df1de;        GT:MD:UD:DP     0/1:4:5:11
chr1    3565870 .       c       N       .       PASS    RS=P;MR=a08605cf-0ff5-4566-8c95-ee00a57ad49e,1f5fb8e4-1bbb-4130-b630-9027ad3f716b,710e99da-731d-4a4a-ba77-f6f364145d02,e7261e07-2389-48c7-9a8c-13ba258df036,9e517f93-0597-47
3a-a576-ad2de1bfe16c;NR=359f0394-1ca1-4d4c-a594-09c5d026164d,0c836564-de22-46b4-894d-80223ead1725,c2b22a84-d8e3-491f-ae2b-478554c02674,a998278c-25a9-4bb1-9753-bcbd42c8b44d,e713cccd-3c1f-4959-8fbd-170314b4be79,9b2bf6df-7b41-47e8-
a485-5c0f6c57e534,f6f2f70b-13e1-4ec1-b395-b2b4867cc53c; GT:MD:UD:DP     0/1:5:7:15
chr1    3565871 .       g       N       .       PASS    RS=N;MR=590c358f-c285-48c1-99c0-cbdc9c02a5ca,a95e8580-f9db-448e-b19c-f0cbb7f3fd35,089f7968-5b6f-49d6-ac27-85fdfc5cbbaa,075dd019-8555-4bca-b1af-d00ce2fb8379;NR=8e01b2b8-4489
-45ca-b332-b484d426a25a,60beabd9-e8bf-4aa7-9c1c-1e415ffd2f8e,633cf7ef-200d-4bfc-8cc7-b608a8f5e569,99ea111d-4c59-4f5e-ba4e-b7e3f0fdcea4,bfca7804-8700-4532-80a9-99fec633d4f3,38308d07-8491-4738-b216-e99170e8ae4d;   GT:MD:UD:DP     
0/1:4:6:10
chr1    3565873 .       c       N       .       PASS    RS=P;MR=a08605cf-0ff5-4566-8c95-ee00a57ad49e,1f5fb8e4-1bbb-4130-b630-9027ad3f716b,710e99da-731d-4a4a-ba77-f6f364145d02,e7261e07-2389-48c7-9a8c-13ba258df036,9e517f93-0597-47
3a-a576-ad2de1bfe16c;NR=359f0394-1ca1-4d4c-a594-09c5d026164d,0c836564-de22-46b4-894d-80223ead1725,c2b22a84-d8e3-491f-ae2b-478554c02674,4a7c25d4-9517-48a5-a29e-db498c1264f9,a998278c-25a9-4bb1-9753-bcbd42c8b44d,e713cccd-3c1f-4959-
8fbd-170314b4be79,9b2bf6df-7b41-47e8-a485-5c0f6c57e534; GT:MD:UD:DP     0/1:5:7:15
chr1    3565874 .       g       N       .       PASS    RS=N;MR=590c358f-c285-48c1-99c0-cbdc9c02a5ca,a95e8580-f9db-448e-b19c-f0cbb7f3fd35,089f7968-5b6f-49d6-ac27-85fdfc5cbbaa,075dd019-8555-4bca-b1af-d00ce2fb8379;NR=8e01b2b8-4489
-45ca-b332-b484d426a25a,60beabd9-e8bf-4aa7-9c1c-1e415ffd2f8e,633cf7ef-200d-4bfc-8cc7-b608a8f5e569,99ea111d-4c59-4f5e-ba4e-b7e3f0fdcea4,bfca7804-8700-4532-80a9-99fec633d4f3,38308d07-8491-4738-b216-e99170e8ae4d;   GT:MD:UD:DP     
0/1:4:6:10

showing a single read from the modified bam file

88ada23edc17842dccb104ed02d30eb1f_dna_r10.4.1_e8.2_400bps_sup@v4.1.0        NM:i:20 AS:i:4124       sd:f:26.9664    zd:i:1  de:f:0.00754717 ch:i:341        rl:i:14005      cm:i:128        sm:f:102.699    fn:Z:sample_rate-4000.po
d5      nn:i:0  rn:i:48266      tp:A:P  ms:i:4128       ns:i:434025     qs:i:20 ts:i:485        st:Z:2023-06-09T00:14:41.762+00:00      du:f:108.506    sv:Z:quantile   dx:i:0  mx:i:4
bb914429-7312-4522-a8ac-438002d4ea93    2048    chr1    3000001 60      12077H19M1I78M2I161M1I287M1I9M1D373M1D18M1I229M1D59M2D142M2I5M1I21M1I14M1I74M2I1M2I26M1I2M2I605M26H     *       0       0       TTCTGTTTCTATTTTGTGGAATAATTTG
AGGAGAGTTGGAATTAGGTCTTCTTTGAAGGTCTGGTAGAACTCTGCATTAAACCCATCTGGTCCTGGGCTTTTTTTTTTTTTTTTTTTTTTTTTTTTGGTTGGGAGACTATTGATGACTGCCTCTATTTCTTTAGGGGAAATGGGACTTTTAGTCCATGAATCTGATCCTGATTTAGCTTTGGTACCTGGTATCTGTCTAGGAAGTTGTCCATTTCATCCAGGTTTT
CCTGGTTTTTTTTTTAGTATAGCCTTTCATAGTAGAATCTGATGATGTTTTTGATATCCTCATGTTCTGTTGTTATGTCTCCTTTTTCATTTCTGATTTTGTTAATTATAGTACAGTCCCTATGCCCTCTAGTTAGTCTGGCTAAGGGTTTATCTATCTTGTTGACTTTCTCAAAGAACCAGCTACTAGTTTGGTTGATTCTTTGAATATTTCTTTTTGTTTCCACTT
GGTTGATTTCAGCTCTGAGTTTGATTATTTCCTGCTGTCTACTCATCTTGGGTGAATTTGCTTCCTTTTTGTTCTGAGCTTCTAGATTTGCTGTCAGGCTGCTAGTGTATACTCTAGTTTCCTTTTGGAGGCACACAGGCCTGTGAGTTTTACTCTTAGGACTGCCTCATTGTGCCCCATATGTTTGGCTATGTTGTGGATTTATTTTCATTAAACTTTAAAACATCT
TTAATTTTTTTCTTTATTTCATCATTGACCAAGCTATCATTAAGTAGAGTATTGTTCAGTTTCCAGGTGAATGTTGGCTTTCTATTATTTATGCTGTTATTGAAGATCAGCCTTAGTCCGTAGTTATCTGAAAGGATGCATGGGAAAATTTCAATATTTTTGTATCTGTTGAGGACTTTTTGTGAGTGACTATATGGTCAATTTTGGAGGATTTGGTACTAGAAGAAG
GTATATATCCTTTTTGTCTTATGATAAAATGTTCTGTAGATATCTATTAAATTCATTTGTTTCATAACTTCGGTTAGTGTCTGTGTGTCTCTGTTCAGTTTCTGCTTCCAGGATCTGTCTCTTGGTGAGAGTGTGGTCTTGAAGTCTCCCAGTATTATTTTATGAGGTGCAATGTGTGCTTTGATCTTTAGCAAAGTGTATTTAATGAATGTGGCTGCTCTTGCATTT
AGAGCATAGACATCAGAATTGAGAGGTCATCTTGGTAGATTTTGCCTTTGATGAGTATGAAGTGCCCCTCATTTTTTTTGATAACTTTGAGCTGAAAGTAAATTTGGTTCGATATTAGAATTGCTACTCCAGCTTAATTCTTCATCCCATTTGCTTGGAAATTTGTTTTCCAACCTTTAACTCTGAGGTAGTGTCTGTCTTTTTCCCTGAGGTCAGAGTTATCCTGTA
AGCAACAAAATGTTCAGGTCCTGTTTGTGTTAGCCAGTCTATTAGTCTATGTCTTTTTACTGGGGAATTGAGTCCATTGATGTTAAGAGAAATTAAAGACAAGTGGATTATTGTTGCTTCCTGTTATTTTTGTTGCTTGGAAACCTGGGATTCTGTTCTTGAGGCTGTCTTCTTTTAGGTTTGTTAAAGGATTGCATTCTTGCTTTTTCTAAGGTGTAGTTTCTGTCC
TTGTGTTTGTGTTATCCCGTTATTATCCTTTGAAAGGCTAGATTCATGGAAAGATATTGTGTGAATTTTGTTTGGTCATGGAATACTTTGGTTTCTCCATCTATGGTAATTGAGAGTTTGGCTGGGTATAGCAACCTGGGCTGGCACTTCTGCTTTCTTAGGGTCTGTATAACATCTGTCCAGGATCTTCTGGCTTTCATAGTCTCTGGTGAGAAGTCTGGTGTAATT
CTAATAGGCCTGCATTTATATGTTACTTGACCTTTTTCCCTTACTGCTAAAGATATTCTATCTTTATTTAGTGAATTTGTTGTTCTGATTATTATGTGTCAGGAGGAATTTCTTTTCTGGTCCAATCTATTTGGAATTCTGTAGGCTTCTTTTATGTTCATGGATATCTCTTTCTTTAAGTTTGGGAAGTTTTCTTCTCTAATTTTGTTAAAGATATTTGCTGGTCCT
TTAAGTTGAAAATCTTCATTCTCACCTACTCTTGTTGTCCATATGTTTGGGCTTTTCATTG   8879;:<;*****C?9B><<>>=AA@ACA><<::;<<><===>?===9891/.-(((&'$$$$%%())))*11111BAAAA:9989>?@@@><===A@==:6,./147;>BDFGIGFFFGGHHDAA855568734;;:=>?????>=<<:89889;;A-==3..
8899:=?@>;;889:D8766667:=>?98888=::::<AC>=:9::;=;;::;8888899::;>??>>?>=>:;;:45<=:<<13++888>BD;;:98;9>=EDCEB<9978555677<9@??<;<<<@@AABBB@????ABBK@@@?@=98778>>?>>?>??>?>>>>;:88876SEF;;;;??>>@CCFBAABBGIFCA><;:::;987779:974434477;<=
@>=;;;<;@<77777:;33---.../75579<60155599?@>=?=:::8874221-111233336666@DABDCD>>>>EG@??@=?=::<;89989<=@@?=<;<<=?>@@@>?>?>@BD==:98:==:;977779<=A?>>>>>@EC?;9998786321133..+**+/41,,,5:::88;::;<<<;;:;<=>DFJHCF?>=<544555755884;;97888;?
@?;:8889<76(((((9;222//09::876767797766551/3454100-035789>>>>;9899;@@@B@C>>=?CB>===>FFE>=>;;==:CEDDGJHGS@??<@FF=;;99;<<<98999;;;<<>>?@G???==<<<==BBC@@>>?@>>:9::::=?<==<;::;<=A?@@ABBDDG@>>=<===???@<==<<<:99:;:;8100015:;;<===>A=<9
989;;7778=?EDEC@@@??CEBAC>==>;;;;;==;9889;>DE@@=<:::::;;<<@>:966666BC?>=;9:::<=<<877+)('--076999:=<>:4223=;;@32223<<<9:=;<<<<<=CA@;::::=;;33333=>>===<?>??:;?>>???611124;<<<;????A?@<;33322,,677?>C>:::99<:89979997:++++057;<=<===::
<9:;;<A?@<=77766;4=???><8-99;9==;;<==>?>=<:9:::?@==>=;:;;;<;;:::>?ABFDBA@@?@><::98999;::999992221././////)'(34552..011;::;;<>;:><8999:BB=;::==DF33333;99:;;<9765667:>9559:C=<<=>888623349;;=@>=<<;=>>,+**876444DA@@A@A@=<;88888<;=>>
>?A@>=88987878=@A;;;;9<<A<EE@@@@@H><<;;<<?=C?>==:;87777888:<;=;<;<A@@=;499998-*('&'((()))7:=>;::::<??CE665.....077650/.//78::;?>;;:::==@<<;::=<><<<=DE=;8//24:;<><;998889:=??>=>==?@CDC?>98746965556*(&&&''(()*))+//0::=>==?@?=744..
1*%%$%&%%%%%).4522233336<;=<;;;<89888<<@>=;;<<;9988;:;889;<;;;:<=;8789?@>>:::<<<;;==>?@AA?:7788:?>>?>?>=>>@<910/004@@?;999;=BAA>=999:;=>?A@>==>A<;;;;C@@>@>AAAACDAABA=;;;;<>;;;;<>>?>=;;;;8;:99=:<<<;;;;<=@@?==;;<=<>=<<;;;;;:877889
;;:;=<<<;;:66667<9;;376)(((*+++*+.))+**(('&&)*...-+,,+,,/0::;=@B@??>==>>=>>>>?@::::::::9;??C????>?=;876778;;>BCDC@??>:::99:877:GS<<98<:;@?==?>A???>>D<;;;9;;<<:<7222229<<DBAB@?A?>>>?AAEDA@??7765578:;<>AFHA@@@FCE@>>=;;;;;>==>D@=<<
<=AAA@??;87788@AFDE@@A?@A=::::<>?<;:;:AB@<BCCFF=:8669=@AACBDA@=<<==DDICAA@@C;;;;:AB@=<;99766764695655679?86666;==><<:9999;<<===@==:::;;@A@AAA>=8889HC@??>   s1:i:1371       s2:i:222        SA:Z:chr8,29938254,+,1S3053M9D11190S,1,6
3;chr9,102364003,+,3640S1110M1I9493S,17,222;chr18,73591923,-,7155S647M6I6436S,36,77;    MD:Z:19T2C102G160A37G115T113^A373^G88C158^T59^TT143G24G119T0T600        PG:Z:MarkDuplicates     RG:Z:8f728d488ada23edc17842dccb104ed02d30eb1
f_dna_r10.4.1_e8.2_400bps_sup@v4.1.0        NM:i:34 AS:i:4066       sd:f:25.8816    de:f:0.0130841  ch:i:1645       rl:i:4319       cm:i:117        sm:f:105.673    fn:Z:sample_rate-4000.pod5      nn:i:0  rn:i:3456       tp:A:P  
ms:i:4071       ns:i:149980     qs:i:20 ts:i:420        st:Z:2023-06-06T22:36:04.013+00:00      du:f:37.495     sv:Z:quantile   dx:i:0  mx:i:4
twolinin commented 3 weeks ago

Hi @sahuno ,

If you want to obtain a phased modcall VCF, could you please confirm whether you ran longphase phase after completing longphase modcall?

longphase phase \
-s SNP.vcf \
--mod-file modcall.vcf \
-b alignment.bam \
-r reference.fasta \
-t 8 \
-o phased_prefix \
--ont

Thanks

sahuno commented 2 weeks ago

Thank you @twolinin for the heads up! the output from the command below is phased D-0-1_4000.vcf

longphase phase -s results/call_snps_indels/D-0-1_4000/snv.vcf.gz \
--mod-file results/longphase_modcall/D-0-1_4000/modcall_D-0-1_4000.vcf \
-b /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/results/mark_duplicates/D-0-1_4000/D-0-1_4000_modBaseCalls_sorted_dup.bam \
-r /data1/greenbab/database/mm10/mm10.fa \
-t 12 \
-o results/longphase_phase/D-0-1_4000/D-0-1_4000 \
--ont

My goal is to do 5mC + SNP co -phasing. However my final haplotagged bam file doesn't have phasing information. i've checked manually and by command

# samtools view /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/D-0-1_4000_haplotagged_new.bam | awk '{for(i=12;i<=NF;i++) if($i ~ /^HP:|^PS:/) print $i}'

this is command i used

longphase haplotag \
-s s4000/results/call_snps_indels/D-0-1_4000/snv.vcf.gz \
--mod-file s4000/results/longphase_phase/D-0-1_4000/D-0-1_4000.vcf \
-r mm10.fa -b /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/results/mark_duplicates/D-0-1_4000/D-0-1_4000_modBaseCalls_sorted_dup.bam \
-t 12 -o D-0-1_4000_haplotagged_new

please see snapshot of phased modcall file

cat /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/snps_longphase_modcalls/s4000/results/longphase_phase/D-0-1_4000/D-0-1_4000.vcf | grep -v "^##" | head
chr19   3079763 .       A       C       24.5814 PASS    H;FAU=6;FCU=4;FGU=0;FTU=0;RAU=5;RCU=1;RGU=0;RTU=0;SB=0.60268    GT:GQ:DP:AF:AD:AU:CU:GU:TU:PS       0|1:24:16:0.3125:11,5:11:5:0:0:3079248
chr19   3086424 .       G       C       19.2765 PASS    FAU=0;FCU=3;FGU=15;FTU=0;RAU=0;RCU=2;RGU=7;RTU=0;SB=0.39155     GT:GQ:DP:AF:AD:AU:CU:GU:TU:PS       0/1:19:27:0.1852:22,5:0:5:22:0:.
chr19   3093019 .       C       G       18.3241 PASS    H;FAU=0;FCU=7;FGU=6;FTU=0;RAU=0;RCU=6;RGU=4;RTU=0;SB=1.0        GT:GQ:DP:AF:AD:AU:CU:GU:TU:PS       1|0:18:23:0.4348:13,10:0:13:10:0:3079248
chr19   3093020 .       T       C       24.2024 PASS    H;FAU=0;FCU=7;FGU=0;FTU=11;RAU=0;RCU=4;RGU=0;RTU=7;SB=1.0       GT:GQ:DP:AF:AD:AU:CU:GU:TU:PS       1|0:24:29:0.3793:18,11:0:11:0:18:3079248
chr19   3103636 .       T       A       32.2231 PASS    H;FAU=5;FCU=0;FGU=0;FTU=5;RAU=3;RCU=0;RGU=0;RTU=8;SB=0.65944    GT:GQ:DP:AF:AD:AU:CU:GU:TU:PS       1|0:32:21:0.3810:13,8:8:0:0:13:3079248
chr19   3109477 .       A       T       9.7446  LowQual H;FAU=4;FCU=0;FGU=0;FTU=2;RAU=2;RCU=0;RGU=0;RTU=2;SB=0.6044     GT:GQ:DP:AF:AD:AU:CU:GU:TU:PS       0|1:9:14:0.2857:6,4:6:0:0:4:3079248
chr19   3109520 .       G       A       15.9325 PASS    FAU=1;FCU=0;FGU=3;FTU=0;RAU=3;RCU=0;RGU=2;RTU=0;SB=0.03571      GT:GQ:DP:AF:AD:AU:CU:GU:TU:PS       0|1:15:15:0.2667:5,4:4:0:5:0:3079248
chr19   3129778 .       A       T       31.7292 PASS    H;FAU=3;FCU=0;FGU=0;FTU=3;RAU=5;RCU=0;RGU=0;RTU=5;SB=1.0        GT:GQ:DP:AF:AD:AU:CU:GU:TU:PS       0|1:31:16:0.5000:8,8:8:0:0:8:3079248
bash:islogin01:/data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed 1025 $

i called somatic variants with

##cmdline=/opt/bin/run_clairs_to --tumor_bam_fn /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/results/mark_duplicates/D-0-1_4000/D-0-1_4000_modBaseCalls_sorted_dup.bam --ref_fn /data1/greenbab/database/mm10/mm10.fa --threads 12 --platform ont_r10_dorado_sup_4khz --output_dir /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/snps_longphase_modcalls/s4000/results/call_snps_indels/D-0-1_4000 --ctg_name chr19 --conda_prefix /opt/micromamba/envs/clairs-to

i'm i supposed to phase D-0-1_4000/snv.vcf.gz separately before using for haplotagging? did i miss something?

benson10120 commented 2 weeks ago

Hi @sahuno , If you want to do 5mC + SNP co-phasing, could you please try 

longphase phase -s results/call_snps_indels/D-0-1_4000/snv.vcf.gz \
--mod-file results/longphase_modcall/D-0-1_4000/modcall_D-0-1_4000.vcf \
-b /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/results/mark_duplicates/D-0-1_4000/D-0-1_4000_modBaseCalls_sorted_dup.bam \
-r /data1/greenbab/database/mm10/mm10.fa \
-t 12 \
-o results/longphase_phase/D-0-1_4000/D-0-1_4000 \
--ont

After phasing , you will get two output files, D-0-1_4000.vcf and D-0-1_4000_mod.vcf Haplotag using the two phased output files.

longphase haplotag \
-s s4000/results/longphase_phase/D-0-1_4000/D-0-1_4000.vcf \
--mod-file s4000/results/longphase_phase/D-0-1_4000/D-0-1_4000_mod.vcf \
-r mm10.fa -b /data1/greenbab/projects/triplicates_epigenetics_diyva/DNA/preprocessed/results/mark_duplicates/D-0-1_4000/D-0-1_4000_modBaseCalls_sorted_dup.bam \
-t 12 -o D-0-1_4000_haplotagged_new

Thanks

sahuno commented 1 week ago

your approach worked! thanks!