I want to count the number of reads falling into specified regions. Why is there a difference when I used -F 4 -L and just -L.
I used
samtools view -c -L ../../../gene_matrix/matrix_gene.bed ENCFF001KYF.bam
31168853
samtools view -c -F 4 -L ../../../gene_matrix/matrix_gene.bed ENCFF001KYF.bam
14868564
samtools view -c -F 4 ENCFF001KYF.bam
23240270
So, the number of mapped reads is 23240270. And when I search for the number of reads in a specified regions not taking into account -F option, I am getting more reads (31168853) then I have the mapped reads. How is it possible? Why do I get some strange reads that fall into my region which flag is not set to 4? I though if the read is in the region then it will be for sure mapped.
2 answers
If bit 4 of the flag is set, the chromosome and position information can be ignored. Many aligners will give a position of * and chromosome of * for unmapped alignments, but others won't. It is rather annoying when aligners give unmapped alignments position information though, since it screws up people's pipelines.
FYI, picard will complain about stuff like this (i.e., stuff that follows the spec. but doesn't follow what really should ideally be done).
I'm not sure for your exact case, but some aligners break the "rules" (or at least your reasonable interpretation of them) and give mapped locations for reads that are flagged as unmapped. This could be causing what you're seeing.
See this FAQ from Picard about bwa, for instance: "BWA can produce SAM records that are marked as unmapped but have non-zero MAPQ and/or non-"*" CIGAR. Typically this is because BWA found an alignment for the read that hangs off the end of the reference sequence."
You can examine some of your reads manually to see if this is the case, or look at the output from the complementary -f 4 query.
Log in to answer this question.
-
bam file: unmapped but matching reads (samtools view -f 4)
written by Tommy •
-
Differing number of reads for read 1 and 2 in fastq's from a subset bam file
written by michael •
-
extracting mapped reads from the bam file
written by KVC_bioinfo
-
Larger number of reads after filtering 'chrM'??
written by cpc5 •
-
Extracting mapped paired end FastQ reads from BAM file using Samtools - trying to understand -f 0x0…
written by acxcrew007 •
-
samtools counting unique mapped reads
written by bitpir
-
Getting read counts per loci
written by JameB11 •
-
samtools view, overlap with multiple regions
written by tonja.r
-
featurecounts vs samtools view
written by tonja.r
-
coverage per base using samtools depth and mpileup
written by tonja.r
