This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Is this range of VAF is ordinary?

Hi

I have calculated Variant Allele Frequency (VAF) for called SNVs and INDELs called by Strelka separately. For getting VAF, I done

VAF = Tumour Variant Allele Count / Tumour Read Count

For some position of the genome I have VAF > 30 , so DOES these big VAFs are normal or I am doing something wrong?

I was supposing VAFs should be in the range of 0 < VAF <1

Can you help in getting some idea?

Thanks

wgs vcf vaf
• 2,306 views • link
written by zizigolu zizigolu

I was supposing VAFs should be in the range of 0 < VAF <1

Yes. I don't see how VAF can be >1. There's something wrong in your calculations.

• 0 views • link
written by Ram Ram

Did you follow the instructions to calculate VAFs as suggested in one of your previous posts (=as in the Strelka manual)? I guess not given this result.

• 0 views • link
written by ATpoint ATpoint

Actually somebody wrote a script for me; Assuming a Strelka .vcf for SNV

#Basic information
chrom=$(echo "$line" | sed 's/ /\t/g' | cut -f 1) #&& echo $chrom;
Pos=$(echo "$line" | sed 's/ /\t/g' |  cut -f 2) #&& echo $Pos;
Ref=$(echo "$line" | sed 's/ /\t/g' | cut -f 4)
Alt=$(echo "$line" | sed 's/ /\t/g' | cut -f 5)

#Tumor sample read, variant and reference information
ReadCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 13 | sed 's/ReadCount=//' )
VariantAlleleCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 26| sed 's/VariantAlleleCount=//')
ReferenceAlleleCount=$(echo $ line | awk -v rc="$ReadCount" -v vac="$VariantAlleleCount" '{print rc-vac}')


#Control or Normal read, variant, reference information
ReadCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 14 | sed 's/ReadCountControl=//')
VariantAlleleCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 27 | sed 's/VariantAlleleCountControl=//')   

ReferenceAlleleCountControl=$(echo "$line" | awk -v rcc="$ReadCountControl" -v vacc="$VariantAlleleCountControl" '{print rcc-vacc}')


VAF=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 28 | sed 's/VariantAlleleFrequency=//')
• 0 views • link
written by zizigolu zizigolu

Somebody wrote a script for me

No wonder you're running into problems you can't explain. You can execute each statement line by line and see where the logic goes awry, or you can contact the author and hope they have the time to explain what could be going wrong. I'd recommend the former approach.

• 0 views • link
written by Ram Ram

Please read the Strelka manual towards calculating AFs. The way one calculates this is different for Indels and SNPs, so running one script (the one below) on both is not going to work. For Indels you need the TIR and TAR values while for SNPs you will have to extract something different, I do not remember. I fugired it out back in the day entirely by reading the manual, I am sure you can do that as well.

• 0 views • link
written by ATpoint ATpoint

For INDELs he has written 

#Basic information
chrom=$(echo $line | sed 's/ /\t/g' | cut -f 1) #&& echo $chrom;
Pos=$(echo $line | sed 's/ /\t/g' |  cut -f 2) #&& echo $Pos;
Ref=$(echo $line | sed 's/ /\t/g' | cut -f 4)
Alt=$(echo $line | sed 's/ /\t/g' | cut -f 5)

#Tumor sample read, variant and reference information (I used TIER 1)
ReadCountTumor=$(echo $line | sed 's/ /\t/g' | cut -f 11 | sed 's/:/\t/g' | cut -f 1)
SupportAltAlleleTumor=$(echo $line | sed 's/ /\t/g' | cut -f 11 | sed 's/:/\t/g' |  cut -f 3 | sed 's/,/\t/g' | cut -f 1)
SupportIndelTumor=$(echo $line | sed 's/ /\t/g' | cut -f 11 | sed 's/:/\t/g' | cut -f 4 | sed 's/,/\t/g' | cut -f 1)

#OMMITTEDSupportOtherTumor=$(echo $line | sed 's/ /\t/g' | cut -f 10 | sed 's/:/\t/g' | cut -f 5 | sed 's/,/\t/g' | cut -f 1)


AltAlleleFrequencyTumor=$(echo $line | sed 's/ /\t/g' | awk -v RCT=$ReadCountTumor -v AltAlleleTumor=$SupportAltAlleleTumor '{print AltAlleleTumor/RCT}')                                                                         
IndelFrequencyTumor=$(echo $line | sed 's/ /\t/g' | awk -v RCT=$ReadCountTumor -v AltINDELTumor=$SupportIndelTumor '{print AltINDELTumor/RCT}')  



#Control or Normal read, variant, reference informationi (I used TIER 1)
ReadCountControl=$(echo $line | sed 's/ /\t/g' | cut -f 10 | sed 's/:/\t/g' | cut -f 1)
SupportAltAlleleControl=$(echo $line | sed 's/ /\t/g' | cut -f 10 | sed 's/:/\t/g' | cut -f 3 | sed 's/,/\t/g' | cut -f 1) 
SupportIndelControl=$(echo $line | sed 's/ /\t/g' | cut -f 10 | sed 's/:/\t/g' | cut -f 4 | sed 's/,/\t/g' | cut -f 1)

#OMMITTED##SupportOtherControl=$(echo $line | sed 's/ /\t/g' | cut -f 10 | sed 's/:/\t/g' | cut -f 5 | sed 's/,/\t/g' | cut -f 1)

AltAlleleFrequencyControl=$(echo $line | sed 's/ /\t/g' | awk -v RCC=$ReadCountControl -v AltAlleleNormal=$SupportAltAlleleControl '{print AltAlleleNormal/RCC}')

IndelFrequnecyControl=$(echo $line | sed 's/ /\t/g' | awk -v RCC=$ReadCountControl -v AltINDELNormal=$SupportIndelControl '{print AltINDELNormal/RCC}')

He believes that

In cases where the frequency is above 100% or 1, this is likely an error where there is more information in support of the variant than there is read depth???. In these cases, you could consider the frequency to be around 1.

• 0 views • link
written by zizigolu zizigolu

Is the “he” you mention the author of the Strelka manual? If not, your comment does not add to the discussion, it’s just asking ATPoint to debug a random person’s script. That’s your job, not ATPoint’s.

• 0 views • link
written by Ram Ram

Sorry in my 5 years membership in biostars I realised that most of times you used to commenting on OP's posts in a harsher way than the other mediators used to. I admit that I am too bad in deep and logical thinking but I don't think this would justify the way you are commenting. By mentioning "HE", I just meant somebody who has written this script. However, by mentioning this I solely did not mean that I expect @ATPoint to debug this script for me rather firstly I wanted to share this script for people who may face the same problem with INDELs in future. I then posted if @ATPoint wants, of course he can argue this script which would be a lot of help.

• 0 views • link
written by zizigolu zizigolu

I apologize for the harsh comment. I take it all back. Good luck with the project.

• 0 views • link
written by Ram Ram

Not at all. You are very welcome to comment

• 0 views • link
written by zizigolu zizigolu

Sorry @ATpoint, I googled but I failed to find a full documentation explaining what each part of a vcf from Strelka means especially INFO column. Please can you share if you found such documentation?

• 0 views • link
written by zizigolu zizigolu

0 answers

No answers yet.

Log in to answer this question.

More posts like this