Hot answers tagged

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.


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 ...


4

That documentation has been expanded recently (see PR #1055) and now describes the characters seen in this column in more detail: Forward Reverse Meaning --------------------------------------------------------------- . dot , comma Base matches the reference base ACGTN acgtn Base is a mismatch to the reference base > ...


3

Examining the mpileup output for a single read on its own is a good way to figure out why bases are not appearing as expected. So if you suspect that the cigar=18S58M read is the missing one (you are correct), prepare a SAM file containing just that read: @SQ SN:chr1 LN:248956422 NS500355:NS500355:HHYN5BGXB:1:13209:10692:11045 83 chr1 569888 53 ...


3

Why do you believe it's the read with cigar=18S58M? The mapping qualities given in sam files are different to those in the mpileup output. I guess samtools is doing some recalculation. You are using the parameters -q20 -Q20. So you skip all reads with mapping quality <20. This will result in 7 reads got to mpileup. Furthermore you are skipping bases with ...


2

By using -h in the samtools view command, you're including all the header lines in your word count. If you happen to have about 89500 reference sequences, then the lengths of those would all appear in the header and inflate the -h word count, but not the mpileup count. Try piping it through an additional samtools view (i.e. without -h) and see if the counts ...


2

Why not use tee? Please see this post. You can write mpileup to file while also piping it to another process. For example, samtools mpileup ... | tee file.mpileup | bcftools call ... the bcftools will run while generating the file.mpileup. After it's done, you can use file.mpileup for your VarScan2. You can test it with the following: touch helloworld.txt ls ...


2

I broke down and finally scratched this itch. I have implemented a tool that performs read-group-aware text-pileup from a single SAM/BAM file. The tool is called streaming_pileupy and it's available from pypi for installation with pip. After installation the command would be: spileup input.bam sample_names.txt The tool is pre-alpha and missing some features....


2

The main issue here is that your bam files have different chromosome labels. (i.e. 1,2,3 vs. chr1,chr2,chr3) as you mentioned. This suggests that the data was mapped to a slightly different version of the reference genome. As a first layer, you might want to double check that these reference chromosomes are cross compatible (e.g. both correspond to hg19 for ...


1

Bcftools and samtools both follow the very common Unix filter pattern: they read from all input files specified (as non-option plain arguments), or if there are none they read from standard input. bcftools call is no exception: $ bcftools call About: SNP/indel variant calling from VCF/BCF. To be used in conjunction with bcftools mpileup. ...


1

Skipped references are similar to deletions, but the different symbols used are indicating that it's an expected deletion. The most common use for a skipped reference is when excluding intronic sequence from cDNA reads that are being matched to a genome reference. Note that the skipped reference bases are added by the mapping program, which generates the ...


1

I am not familiar with getting variants from bam files, so I will just answer the second part. From the error message, it seems that it is necessary to have consistent chromosome labels across datasets. It is also a good idea to do so because this keeps your datasets consistent and makes things a lot easier. To change the chromosomes labels, you can write ...


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