Closed ke-shi closed 6 years ago
Hi, Your samtools/VarScan2 pipeline produces AD, so you should be okay with that. For a samtools/bcftools pipeline we found the following to work for us (at least for v1.1):
samtools mpileup -f genomic.fa -I -u -g -B -t DP,DPR,DV,DP4,INFO/DPR,SP result.bam | bcftools call -cv - | bcftools view -M2 - > result.vcf
The key piece is samtools mpileup -t DP4
I hope this helps. Cheers, Rudiger
Hi Rudiger,
Thanks for the comments. I retried with samtools v1.3.1 to generate a vcf file including DP,DPR,DV,DP4,INFO/DPR,SP, as you suggested:
samtools mpileup -t DP,DPR,DV,DP4,INFO/DPR,SP
but GUSMap stopped again after the VCFtoRA command with the message like this:
Error in VCFtoRA("out.vcf") : object 'ad_pos' not found Execution halted
Next, I retried again to make a new vcf with the samtools command
samtools mpileup -t DP,DPR,DV,DP4,INFO/DPR,SP,AD,AO,RO
but samtools said and died:
[warning] tag DPR functional, but deprecated. Please switch to
AD
in future. [warning] tag DV functional, but deprecated. Please switch toAD
in future. [warning] tag DP4 functional, but deprecated. Please switch toADF
andADR
in future. [warning] tag INFO/DPR functional, but deprecated. Please switch toINFO/AD
in future. Could not parse tag "AO" in "DP,DPR,DV,DP4,INFO/DPR,SP,AD,AO,RO"
What should I modify?
Thanks, Kenta
Hi Kenta,
There was a bug in the code which is what gave you the first error, e.g.,
Error in VCFtoRA("out.vcf") :
object 'ad_pos' not found
Execution halted
I have made some changes which hopefully has fixed things. Download GUSMap again and try running the VCF created via
samtools mpileup -t DP,DPR,DV,DP4,INFO/DPR,SP
again and let me know if it still breaks.
Thanks, Timothy
Hi Timothy,
I downloaded the new one and tried it with a vcf from
samtools mpileup -t DP,DPR,DV,DP4,INFO/DPR,SP
, but still have an error:
% /app/gen/R-3.5.0/bin/R
library(GUSMap)
Genotyping Uncertainty with Sequencing data for linkage MAPping (GUSMap) Copyright (C) 2017-2018 Timothy P. Bilton This program comes with ABSOLUTELY NO WARRANTY. This is free software, and you are welcome to redistribute it under certain conditions.
Welcome to GUSMap v0.1.1
RAfile <- VCFtoRA('test.s1.vcf')
Processing VCF file: Converting to RA format.
Found 184 samples Error in VCFtoRA("test.s1.vcf") : (list) object cannot be coerced to type 'integer'
Thanks, Kenta
Hmm, there is obviously another bug somewhere. I'll need to find some data generated in the same format to run some tests on.
Timothy
Hi Kenta,
I found the bug causing the latest issue. I've fixed it in commit 642d069. Download GUSMap again and try again.
Sorry about the bugs in the code. I hope there are no more.
Timothy
The latest version worked well!! Thanks!
That's great to hear. Thanks for highlighting the bugs in the code.
What SNP callers do you recommend for making vcf?
library(GUSMap)
RAfile <- VCFtoRA('my.vcf')
readRA(RAfile, pedfile='my_ped.csv')
Then, I got an error message: Error in if (sum(d_p) > filter$DEPTH) { : missing value where TRUE/FALSE neededWe usually use two SNP calling methods to make vcf files, samtools/bcftools and samtools/VarScan2. I assume the vcf formats do not include sufficient information for GUSMap saying "In order to convert the VCF file to RA format, some information regarding allelic depth is required. Specifically, one of the fields 'AD', 'RO' and 'AO', or 'DP4' needs to be supplied in the VCF file."
VCF files from samtools/bcftools include only GT:PL:DP:GQ.
FORMAT=
FORMAT=
FORMAT=
FORMAT=
VCF files from samtools/VarScan2 include GT:GQ:SDP:DP:RD:AD:FREQ:PVAL:RBQ:ABQ:RDF:RDR:ADF:ADR.
FORMAT=
FORMAT=
FORMAT=
FORMAT== 15">
FORMAT=
FORMAT=
FORMAT=
FORMAT=
FORMAT=
FORMAT=
FORMAT=
FORMAT=
FORMAT=
FORMAT=