RNA-SEQ SNP calling

Question

[babji2318]

I am comparing the effect of knocking down a gene in Aedes aegypti, I did a RNA-SEQ to understand expression profiles, I am also using this RNA-SEQ against the genome to check for variants I did a variant analysis using the GATK variant calling program, since GATK 4 doesnt support RNA-seq variant calling I used GATK 3. Used STAR for alignment I also did variant calling using samtools and bcftools pipeline

Is there a way to count the number of SNP's and INDEL's in my samples from the VCF files generated

I am looking for something that is going to say for example: Total SNP's detected: 3450, Total INDEL's detected: 460 Something along those lines so I can have a grasp of what is happening in the treatment samples

[ATpoint]

Please use the search function:

Answer 1

[JC]

A simple way is to parse and count the elements in the file, for example, to get how many variants are in your VCF, you can use grep command to count all lines which are not a header (marked with #):

grep -vc "#" file.vcf

Then you can filter all variant which looks like an SNP, I prefer to use Perl but you can use whatever you want:

perl -lane 'print if ($F[3] =~ m/^[ACGT]$/ and $F[4] =~ m/^[ACGT]$/)' | wc -l

In this case, we filter lines which have the fields #4 (reference call) and #5 (sample call) as a single DNA letter (A, C, G, or T), then we just count them with wc -l

Indels will be the subtraction of SNP from the total variants unless you are calling another type of variants, but GATK reports SNPs and short-Indels by default.