Entering edit mode
7.8 years ago
a.james
▴
240
Dear All,
I need to query a set of chromosomal positions towards a couple of BAM files, and I am using samtools mpileup for this. Following is my command line,
samtools mpileup -L 1000000 -u -A -B -d 100000 -Q 0 -t DP -t AD -t ADF -t ADR -f GRCh37.fa --positions mybed.bed sample1.bam sample2.bam sample3.bam sample3.bam sample4.bam ... sampleN.bam >test.vcf
I could find the positions in my bed file within that sample.bam files when I look into IGV for certain samples but the above command line is not able to identify those positions.
Any help would be great . My aim is to get the read depth of the chromosomal position defined in my bed file from BAM files
You should be able to use
bedtools genomecov
(documentation) orbedtools coverage
(documentation).I have also one more suggestion (hopefully understood the question properly),
http://deeptools.readthedocs.io/en/latest/content/tools/bamCoverage.html
Can you post the error message you get? Does your bed file contains "chr" prefix?
I didn't get any error message. Whereas I got an empty vcf file (I mean there is header but no results). I am supposed to get a vcf file with read depth and all the information I asked in the parameter of command line.
Does your Reference fasta file contain
chr
prefix before chromosome numbers? Does your bed file also contain thatchr
prefix? You get empty VCF probably because your bed files and reference have different designation for chromosome!The prefix for both the query bed file and genome file for chromosomes has no "chr". And so I don't think that the potential reason