Question: Querying positions (bed file) towards BAM files using samtools mpileup
0
gravatar for a.james
2.3 years ago by
a.james190
Germany
a.james190 wrote:

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

rna-seq • 1.2k views
ADD COMMENTlink written 2.3 years ago by a.james190

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) or bedtools coverage (documentation).

ADD REPLYlink written 2.3 years ago by genomax67k

I have also one more suggestion (hopefully understood the question properly),

http://deeptools.readthedocs.io/en/latest/content/tools/bamCoverage.html

ADD REPLYlink modified 2.3 years ago • written 2.3 years ago by EagleEye6.3k

Can you post the error message you get? Does your bed file contains "chr" prefix?

ADD REPLYlink written 2.3 years ago by reza.jabal320

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.

ADD REPLYlink modified 2.3 years ago • written 2.3 years ago by a.james190

Does your Reference fasta file contain chr prefix before chromosome numbers? Does your bed file also contain that chr prefix? You get empty VCF probably because your bed files and reference have different designation for chromosome!

ADD REPLYlink written 2.3 years ago by reza.jabal320

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

ADD REPLYlink modified 2.3 years ago • written 2.3 years ago by a.james190
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 713 users visited in the last hour