Calculate number of multi-mapped reads?

[KAP]

Hi, does anyone know how to calculate the number of multi-mapped reads from your bam file?

This seems like an obvious/silly question but the latest versions of Cufflinks no longer provides this statistic in the basic output and I have been looking everywhere (log files, samtools stats, this website, etc) but can't figure out how to find this out!

Thanks!

[chadn737]

Although primarily to tally up read counts for genes, HT-seq count also adds up the multimap reads "alignment not unique."

[KAP]

Thanks! I'll look into it.

[polyatail]

Something like this might work.

samtools view -F 4 file.bam | awk '{printf $1"\n"}' | sort | uniq -d | wc -l

Explanation
-----------
samtools view -F 4 file.bam
    print all lines in SAM format except those flagged as unmapped

awk '{printf $1"\n"}'
    print the first field in the line (the read name)

sort | uniq -d | wc -l
    sort the read names
    print only the duplicates (i.e. those appearing more than once)
    count the lines

[KAP]

Thanks, great idea and thanks for the detailed, explanation! I'll give it a try

[Johnwang]

Retrieve mapped reads to certain regions (eg. between 120-140 in a gene)?

I have my own unannotated reference data base. With it, the reads were mapped with Bowtie, now I need to retrieve all reads which mapped to certain region between 120-140 nt/each gene, how can I do it with samtools ? Thanks in advance.

[dpryan]

Make a BED file with the gene locations, flank that with bedtools (or awk), and then use it with the "-L" option to "samtools view".

[Johnwang]

more specifics

Would you please give me more details ? Thanks.

[dpryan]

Which part in particular requires more detail?

[Johnwang]

example

>Sg01_138835_A_G
TGCGAAGTATGCCACCAACATGTTTCACCGCATAACCTTGATCACCTTGGACATTCTGCCACAAGCAGACGACGAGTTCCTCAGAAAATTCATGGACAAATTTGGGCGCACGACCATCTA [GATCGAAAATCTTAATCCCG]AGCTAGTGTTGCATTCCATGCTCCTCGTCTTACACTCCAACAAGTTTGCAGAAAGTTTGTCAAAAATCCACTCGGTAGCATGGACGAGCTGCACGAGTGAGCCAAGGGTTACATATAGATGGAAGAAATGTCCAAATTCTGGAAAGAAGTCCGTCAAGCCGAACACAAGCACGAAAAGTGCAAAGCAAACTCCAAGGCCGACTTGCACAAGTCAAACAAGAGGCCCAATTCAGACAAGCGCTAACCTCTCCCCAAAGAGC

I want to retrieve reads aligned to [GATCGAAAATCTTAATCCCG] within gene above.

[dpryan]

If you aligned to the gene itself, then "samtools view alignments.bam Sg01_138835_A_G:120-140". If you aligned to the genome then find out where in the genome that is.

[Johnwang]

how about multiple genes in reference

>Sg01_138835_A_G
TGCGAAGTATGCCACCAACATGTTTCACCGCATAACCTTGATCACCTTGGACATTCTGCCACAAGCAGACGACGAGTTCCTCAGAAAATTCATGGACAAATTTGGGCGCACGACCATCTAGATCGAAAATCTTAATCCCGAGCTAGTGTTGCATTCCATGCTCCTCGTCTTACACTCCAACAAGTTTGCAGAAAGTTTGTCAAAAATCCACTCGGTAGCATGGACGAGCTGCACGAGTGAGCCAAGGGTTACATATAGATGGAAGAAATGTCCAAATTCTGGAAAGAAGTCCGTCAAGCCGAACACAAGCACGAAAAGTGCAAAGCAAACTCCAAGGCCGACTTGCACAAGTCAAACAAGAGGCCCAATTCAGACAAGCGCTAACCTCTCCCCAAAGAGC
>Sg01_281790_T_C
CATCTATCTATTGCTTCAACCATGAAAATAATTATTTTCAAAGGCTTGAGTAATAATGTTTACCAATATAATGCCTATAATAATGTTACAAAATTCTAGTTTATATACCAAGTTTGATCGATGTTTGTAAAATTTAGCATATTCTTGCCTTGAATATCCTTTTGACACACACAAATACCATATAACTAACCCACAATATTCTGACTTATATGGGTGCAGCGGATCATTTGAGCGAAGCTGTGGGCGACAACGAATCCAAGGAACTTTTTGAATCATTTCGAGTCTTCTGCAATAGAAAAAGTGTTAGTACTAATGCATTTTCATCATTCAACAATTAATATAAATATATATATATATAGAGAGAGAGAGAGAGAGAGAAGTGCACTCATGAACAAAATCC
>Sg01_334575_A_G
ATTTTCTTAGTCCGCAAAAGTTAAACAGTATGTTAAGAATATAGGTTAAGAAACTTAAAAGTAAAAATATTTTTATTAAGAGGTGAAAAATTTTGTACAATAACTTGACTTGGAAAATTCTGCAGTTCCTGTTCTTCCAAATGTAATCTATGGACTATGGCAACAACAAAACTCAGAGACGCTACGTCCATTACACTTCACTTGGACAGTAAGTTAGACATCAACCAGGGGCCAAATTCTTGATGAAAAAAGCCTGCACGAACAGCATCTGGGATTAAAACATTCCAGACCTCCTTAGCCGCATAACATTTGGATGGATCATATAATTATATGGATAGATAACTATATGGATCATAATTAATTCAGCATAATTTTTTTTAATACTCTCTAAGCGCCTCTC
>Sg01_402061_A_C
CATAACCTGATAATATATTTTATATATCTACTATGAAAATAAATAGAATATCATATATACTAAGTAAGAAACGATCTGGTGTGTATCAAACAATGAATAATAAAATGTTGTTATCCAATTCTTATATAATATCTTATTTTAAGCAACAATCAAACAATGAATATTTTAAGAACTTATTTATCATGTTCTTATGGAGGCCTCTAGGGTAAGTCGTGACTGTATAGCATAGAAAATAAAGCTGCAGGAGCGTGGTATATCTTGTCAGAGCATGTATGTGTCTTGAGTGTGAAAGAGGGTTGGAGAAAGCAAAGGAATGTTGGAGGTCATTATTGGCTTCAAGCAGAGAGTGCTTACGAACCTATTAAAGTTGGAGCAGTGAGATATATAGTATAGTTGAGCT

If I want retrieve reads mapped to same region on multiple (several thousands) genes, do I need to list all, like this, Sg01_138835_A_G:120-140, Sg01_281790_T_C:120-140 ... ?

[dpryan]

Make a BED file of the regions and use that with the "-L" option.

[Johnwang]

thanks.

I will do that.