5
$\begingroup$

I would like to get some summary statistics on a vcf file from one individual, which has over a million variant calls. I've tried to make a histogram of indel sizes with this command,

vcftools --vcf sample.vcf --out sample --hist-indel-len

but it returns an empty file. I see thousands of deletions in the vcf, although strangely no insertions. A typical line in the vcf file I have looks like

NC_025968.3 168363 191 N <DEL> . . SVTYPE=DEL;STRANDS=+-:19;SVLEN=-1026;END=169389;CIPOS=-9,330;CIEND=-152,9;CIPOS95=-3,61;CIEND95=-32,3;IMPRECISE;SU=19;PE=19;SR=0 GT:SU:PE:SR ./.:19:19:0

Update: I aligned reads to my reference with BWA MEM, then used samblaster to extract discordant pairs and split-reads. I sorted and indexed the three bam files, then used lumpyexpress to generate the vcf.

samblaster -i sample.sam --excludeDups --addMateTags --maxSplitCount 2 --minNonOverlap 20 | samtools view -S -b - > sample.bam

samtools view -h sample.bam | samblaster -a -e -d sample.disc.sam -s sample.split.sam -o /dev/null

lumpyexpress -B ../sample.sorted.bam -S ../sample.split.sorted.bam -D ../sample.disc.sorted.bam -o sample.vcf 

I would like to see the size distribution of deletions, and if possible, figure out why there are no insertions.

$\endgroup$
5
  • 1
    $\begingroup$ Thanks for your question, Joanne, and welcome to Bioinformatics SE. Can you post some lines from your VCF file so that we can get a better idea of what is going on? It is hard to tell without looking at the data directly. $\endgroup$
    – conchoecia
    Oct 16 '18 at 18:29
  • $\begingroup$ Thanks for the welcome! Here is a typical line: NC_025968.3 168363 191 N <DEL> . . SVTYPE=DEL;STRANDS=+-:19;SVLEN=-1026;END=169389;CIPOS=-9,330;CIEND=-152,9;CIPOS95=-3,61;CIEND95=-32,3;IMPRECISE;SU=19;PE=19;SR=0 GT:SU:PE:SR ./.:19:19:0 $\endgroup$
    – Joanne
    Oct 18 '18 at 0:51
  • 1
    $\begingroup$ Hello @Joanne just a note. You can edit your question to add details (like the line I have already added for you) or any other progress or idea you will have regarding the question. Usually more details you give, better chances that there will be somebody who will know the answer. $\endgroup$
    – Kamil S Jaron
    Oct 18 '18 at 9:37
  • $\begingroup$ Just curious, how are you calling this VCF file? If you are using GATK I think there are some options to output only SNPS, or SNPs + indels, or just indels. Also, can you clarify exactly what you mean by summary statistics? Right now it is hard to tell what your question is exactly: IE do you want a size distribution of the deletions that are in there, or do you also want to figure out why insertions are missing from your VCF? If you would like the latter, we'll need info on the commands you used to generate the VCF file. $\endgroup$
    – conchoecia
    Oct 18 '18 at 14:42
  • $\begingroup$ Thanks for the forum protocol tips - I have edited my question to include more details. $\endgroup$
    – Joanne
    Oct 19 '18 at 17:24
1
$\begingroup$

It looks like lumpy looks for breakpoints in your reference. Presumably places where the is a haplotype switch detectable with the data on-hand?

My only thought is that since lumpy is being used to produce the VCF file that it is inherently not looking for insertions. I would suggest running the same BAM file through the GATK best practices pipeline (easier said than done) and vcftools pipeline. If those produce insertions but lumpy does not, it must be due to how lumpy produces the calls. Good luck!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.