18

It's not really possible to convert bam to vcf. bam is a mapping file, it does not contain the information about variants, this information needs to be inferred in process called variant calling. I find important to mention that it's not just a different format of the same thing. How to call variants (a vcf file) from mapped reads (a bam file) is very broad ...


10

I'll follow up to the great answer from Kamil S Jaron: Regarding predicting what the variant ("mutation" is a very loaded term) will do, there are a variety of tools. Chief among these are annovar and VEP. The general idea behind these is to classify the variants according to their overlap with genes, which codons they change (if any), how big that change ...


10

There are two types of INDELs: short indels and long indels. Some put the threshold at 50bp; others choose 1000bp. Short and long indels are called differently. <50bp short indels are called from read alignment directly. Modern indel callers essentially build a multi-alignment of reads and call an indel with enough read supports. Short indels may break ...


10

Bcftools has sample/individual filtering as an option for most of the commands. You can subset individuals by using the -s or -S option: -s, --samples [^]LIST Comma-separated list of samples to include or exclude if prefixed with "^". Note that in general tags such as INFO/AC, INFO/AN, etc are not updated to correspond to the subset samples. bcftools ...


10

Another solution for repeating one command with different parameters is to use parallel. Create a tab-delimited file samples.txt, which contains the sample names and the information for your read groups, e.g.: Sample1 Lib1 Sample2 Lib2 Sample3 Lib3 Then you can use parallel like this: parallel --colsep "\t" 'bwa mem -M -R "@RG\tID:{1}\tLB:{2}\tPL:...


10

TLDR: yes! be careful with someone's genomic data! There are two aspects to this question: can I find a match for a random VCF in a database of genomes (YES) or can I identify a subject who is not in a (public) database. But even when it may not possible to identify the subject itself, using DTC genetic testing sites may enable you to identify (close or ...


9

Insertions and deletions (indels) are one type among many different types of genetic variation, such as single nucleotide variants (SNVs), copy number variants (CNVs), and structural variants (SVs). I'll assume here that you know how indels are defined, but are simple trying to understand the importance of discovering and analyzing them. The goal of indel ...


9

Part 1 : how to detect mutations The keywords you are searching for are "variant calling". Basically you have to map sequencing reads to a reference genome (or gene) and then estimate for each position of the genome if the observed difference of mapped reads and the reference is more likely a sequencing error or a mutation (in genomic glossary - variant). ...


8

For qualitative analysis, you're probably better off using something less granular like IGV or IGB. However, if you really want to look at a couple of reads: If you're willing to ignore sequencing errors, you can inspect the CIGAR string or the MD tag, both of which give information on the alignment of a single read. The CIGAR string gives details on ...


8

If you don't mind a bit of manual counting, then samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam will produce output where you can count the bases for that position. You could, of course, use the command line to do most of that automatically: samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam | cut -f 5 | tr '[...


8

Something like this might do it: #!/bin/bash # make an array to hold your list of samples/base filenames samples=( S1 S2 S3 S4 S5 ) # make a loop to go through the sample IDs for sample in "${samples[@]}" do : # do something with each sample/file, e.g. make a directory for it mkdir "${sample}"_home # do something else with the sample, e.g. unzip ...


8

Yes, definitely identifiable. The combination of ~80 unlinked common SNPs can be fairly unique in the entire human population, let alone the whole VCF file. EDIT: 30 in the original answer is an underestimate. We need ~80 unlinked common SNPs to uniquely identify every individual to a low false positive rate. Here is the derivation based on the approximate ...


7

1. Adapter Trimming One of the first things I do after encountering a set of reads is to remove the adapter sequences from the start and end of reads. Most basecalling software includes some amount of built-in adapter trimming, but it is almost always the case that some adapter sequence will remain. Removing adapters is helpful for mapping because it ...


7

I used this in the past for ChIP-seq data and it generated SNVs: samtools mpileup \ --uncompressed --max-depth 10000 --min-MQ 20 --ignore-RG --skip-indels \ --fasta-ref ref.fa file.bam \ | bcftools call --consensus-caller \ > out.vcf This was samtools 1.3 in case that makes a difference.


7

I wrote a program, ASCIIGenome, that I find handy in cases where you want to have a quick look at genomic data. It's a genome browser for the command line. To view only reads containing mismatches you can use the internal function awk. To filter for reads where the NM tag (number of mismatches) is >0: ASCIIGenome -fa genome.fa aln.bam ... [h] for help: ...


7

I more like to use "ts/tv" for transition-to-transversion ratio. This abbreviation had been used in phylogenetics. When NGS came along, some important developers started to use "ti/tv", but I am still used to the old convention. Why is the expected value for random substitutions for the Ti/Tv ratio 0.5? There are six types of base changes. Two of them ...


6

I think the best method or combination of methods will depend on aspects of the data that might vary from one dataset to another. E.g. the type, size, and frequency of structural variants, the number SNVs, the quality of the reference, contaminants or other issues (e.g. read quality, sequencing errors) etc. For that reason, I'd take two approaches: Try a ...


6

By default, a CRAM you create with samtools is lossless. It typically halves the input BAM in terms of file size. If you want to compress more, you can let samtools convert most read names to integers. You won't be able to tell optical duplicates from read names, but this is a minor concern. You can also drop useless tags depending on your mapper and the ...


5

using vcfilterjs and the following script: function accept(vc) { var i,j; for(i=0;i< vc.getNSamples();++i) { var genotype = vc.getGenotype(i); if(!genotype.hasAD()) continue; var ad = genotype.getAD(); /* loop over AD starting from '1' = first ALT */ for(j=1 ;j< ad.length ;++j) ...


5

I personally don't think BQSR has a huge impact on variant calling, but you don't really need to guess. If you run GATK BQSR, it outputs a table and charts of exactly how much quality scores are adjusted. The adjustment will vary depending on the position in the read and genomic context (previous and following base). In my experience, the difference is a few ...


5

I've been using the LoFreq* caller for exactly this. It is designed to find variants with very low frequency, so is well suited for this type of analysis. LoFreq* (i.e. LoFreq version 2) is a fast and sensitive variant-caller for inferring SNVs and indels from next-generation sequencing data. It makes full use of base-call qualities and other sources of ...


5

There is a recent paper that attempts to do this: ISOWN: accurate somatic mutation identification in the absence of normal tissue controls. In this work, we describe the development, implementation, and validation of ISOWN, an accurate algorithm for predicting somatic mutations in cancer tissues in the absence of matching normal tissues. ...


5

(edit)you can filter the VCF annotations with snpsift, I've also written a VcfFilterSequenceOntology http://lindenb.github.io/jvarkit/VcfFilterSequenceOntology.html I've written vcf2table: http://lindenb.github.io/jvarkit/VcfToTable.html It decodes VEP and SNPeff annotations: >>chr1/10001/T (n 1) Variant +--------+--------------------+ | Key | ...


5

Another approach is htsbox. You can get a candidate list with: htsbox pileup -Cvcf ref.fa -q20 -Q20 -s5 file.bam > out.vcf Here, -q sets min mapping quality, -Q sets min base quality, -v outputs variants only -c outputs VCF, -C gives you base counts on both strands and finally -s5 requires at least 5 high-quality bases to call out an allele. It is ...


5

On the GATK forum they've recommended the population stratified VCF file for this purpose.


5

Check this python script vcf2maf.py


5

You can use annovar to annotate the vcf, then convert it to maf using the function annovarToMaf of maftools bioconductor package.


5

if you have a repeating workflow, I strongly recommend to have a look at workflow management systems like snakemake. I've also wrote a little tutorial on biostars about this topic, which might be useful for you. fin swimmer


4

What's your reference for the definition of STRs? I think it is still ambiguous among the community. Wikipedia states motifs of >=2 bp. However other references include homopolymers as well: Source 1, source 2 and source 3.


4

You can do this in Hail: from hail import * hc = HailContext() (hc.import_vcf('test.vcf') .filter_variants_expr('gs.exists(g => g.ad[1:].exists(d => d > 10))') .export_vcf('filtered.vcf')) This works with any number of samples and will keep the variants where at least one sample has a genotype with an alternate allele support by more than 10 ...


Only top voted, non community-wiki answers of a minimum length are eligible