Querying positions (bed file) towards BAM files using samtools mpileup
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
0 answers
No answers yet.
Log in to answer this question.
More posts like this
-
Filtering for homozygous reference alleles
written by AGE •
-
How to set up my vcf with a minimum depth?
written by apl00028
-
metagenomic assembly coverage, for multiple samples
written by willnotburn •
-
How can I add PV4 or BQ or NS to my VCF file
written by sahar.voghoei •
-
Samtools mpileup: AF1 and DP4 not provided
written by filippo.martignano •
-
Samtools pileup problem.
written by salvatore.digiorgio •
-
Can I run nucleoatac with a consensus peak file?
written by BioinfGuru
-
How to annotate your vcf files with the GT (genotype) field
written by Michi
-
Cuffdiff Command Line For Labels And Replicates
written by Johan
-
Samtools Mpileup Vs Pileup: Multiple Samples And Gt=0/0
written by Pierre Lindenbaum
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
chrprefix before chromosome numbers? Does your bed file also contain thatchrprefix? 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