Multimapped and unique mapped for BAM file of RRBS data
Hi all,
I have a BAM files from RRBS data, and, I want to have number of all reads, number of multimapped and number of unique mapped. I think it can be done with samtools, and, with this we can have all reads.
Thanks in advance.
samtools view -c sample.bam
1 answer
This command is not enough to do all you want. As far as I know 0x100 is that reads that are marked as multimapped. So the following commands will give you everything you need:
All reads :
samtools view -c input.bam
Multimapped:
samtools view -F 0x4 -f 0x100 input.bam | wc -l
unique map reads;
bedtools bamtobed -i input.bam > output.bed
sort -k 4,4 output.bed | uniq -f 3 -u | wc -l
Then you can merge all the information into one.
• 0 views
•
link
written
by
young_bioinformatician
Log in to answer this question.
More posts like this
-
Why Is My Uniquely Mapped Read Count Higher in Samtools Compared to Bowtie2?
written by daffodil •
-
Extract information of mapped reads bam files
written by sata72 •
-
Problem with ONT sequencing pipeline. Probably in the Samtools part
written by Emilio •
-
How to retrieve the number of reads mapped for each chromosome in the BAM file of methylation data
written by sata72 •
-
Total number of trimmed reads less than number of mapped reads
written by kspata
-
Counting uniquely mapped reads - ambiguous results
written by mickey_95 •
-
Questions about STAR and Bowtie (not Bowtie2) mapping
written by debitboro
-
Counting and extracting in a new file the unique mapped reads from BWA 0.7.15-r1140 mapping
written by valopes •
-
Larger number of reads after filtering 'chrM'??
written by cpc5 •
-
samtools counting unique mapped reads
written by bitpir
To get you started. Multiple other threads exist.
Multi-mapped reads: How to filter/view reads mapped to multiple locations using samtools or pysam
Unique reads: Help with Samtools flags to get only primary, unique, non duplicate aligned reads
The exact definition of what is "unique" depends on the aligner. Typically, one sets a MAPQ threshold, for example 20, and everything below is ignored during analysis because mapping quality is too low. You could simply count reads with MAPQ > 20 (samtools view -q 20) and compare to all reads.