SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to extract multi-mapped reads by samtools? mavishou RNA Sequencing 5 12-05-2016 05:27 AM
Number of Mapped reads ninad Genomic Resequencing 2 08-29-2013 04:15 AM
RNAseq and number of mapped reads seqfan RNA Sequencing 7 06-30-2011 02:01 PM
GAII low number of mapped reads aligenie Bioinformatics 3 06-21-2011 11:55 PM
different number of reads mapped to plus strand and minus strand gfmgfm Bioinformatics 2 02-03-2011 10:26 AM

Reply
 
Thread Tools
Old 12-22-2011, 11:39 AM   #1
KAP
Member
 
Location: Adelaide, Australia

Join Date: Mar 2010
Posts: 12
Default Calculate number of multi-mapped reads?

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!
KAP is offline   Reply With Quote
Old 12-22-2011, 12:38 PM   #2
chadn737
Senior Member
 
Location: US

Join Date: Jan 2009
Posts: 392
Default

Although primarily to tally up read counts for genes, HT-seq count also adds up the multimap reads "alignment not unique."
chadn737 is offline   Reply With Quote
Old 12-23-2011, 09:19 AM   #3
KAP
Member
 
Location: Adelaide, Australia

Join Date: Mar 2010
Posts: 12
Default

Thanks! I'll look into it.
KAP is offline   Reply With Quote
Old 12-27-2011, 10:31 AM   #4
polyatail
Member
 
Location: New York, NY

Join Date: Dec 2010
Posts: 25
Default

Something like this might work.

Code:
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
polyatail is offline   Reply With Quote
Old 12-27-2011, 02:38 PM   #5
KAP
Member
 
Location: Adelaide, Australia

Join Date: Mar 2010
Posts: 12
Default

Thanks, great idea and thanks for the detailed, explanation! I'll give it a try
KAP is offline   Reply With Quote
Old 02-16-2017, 11:04 AM   #6
Johnwang
Member
 
Location: NE, USA

Join Date: Jun 2016
Posts: 14
Default 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.
Johnwang is offline   Reply With Quote
Old 02-16-2017, 12:30 PM   #7
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,466
Default

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".
dpryan is offline   Reply With Quote
Old 02-16-2017, 12:33 PM   #8
Johnwang
Member
 
Location: NE, USA

Join Date: Jun 2016
Posts: 14
Default more specifics

Would you please give me more details ? Thanks.
Johnwang is offline   Reply With Quote
Old 02-16-2017, 12:39 PM   #9
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,466
Default

Which part in particular requires more detail?
dpryan is offline   Reply With Quote
Old 02-16-2017, 12:46 PM   #10
Johnwang
Member
 
Location: NE, USA

Join Date: Jun 2016
Posts: 14
Default example

>Sg01_138835_A_G
TGCGAAGTATGCCACCAACATGTTTCACCGCATAACCTTGATCACCTTGGACATTCTGCCACAAGCAGACGACGAGTTCCTCAGAAAATTCATGGACAAATTTGGGCGCACGACCATCTA [GATCGAAAATCTTAATCCCG]AGCTAGTGTTGCATTCCATGCTCCTCGTCTTACACTCCAACAAGTTTGCAGAAAGTTTGTCAAAAATCCACTCGGTAGCATGGACGAGCTGCACGAGTGAGCCAAGGGTTACATATAGATGGAAGAAATGTCCAAATTCTGGAAAGAAGTCCGTCAAGCCGAACACAAGCACGAAAAGTGCAAAGCAAACTCCAAGGCCGACTTGCACAAGTCAAACAAGAGGCCCAATTCAGACAAGCGCTAACCTCTCCCCAAAGAGC

I want to retrieve reads aligned to [GATCGAAAATCTTAATCCCG] within gene above.
Johnwang is offline   Reply With Quote
Old 02-16-2017, 01:21 PM   #11
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,466
Default

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.
dpryan is offline   Reply With Quote
Old 02-16-2017, 01:46 PM   #12
Johnwang
Member
 
Location: NE, USA

Join Date: Jun 2016
Posts: 14
Smile 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 ... ?
Johnwang is offline   Reply With Quote
Old 02-16-2017, 11:01 PM   #13
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,466
Default

Make a BED file of the regions and use that with the "-L" option.
dpryan is offline   Reply With Quote
Old 02-17-2017, 06:07 AM   #14
Johnwang
Member
 
Location: NE, USA

Join Date: Jun 2016
Posts: 14
Default thanks.

I will do that.
Johnwang is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:50 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO