Dear all,
I used Bowtie to map RNA-seq reads on the genome, and got a sam file. I also have a LIST file with many ranges, below shows an example. I would like to exact the reads that fall into the ranges in the LIST file from the sam file.
LIST file
Chr1 + 2-999
Chr1 + 2999-5000
Chr3 + 90-2000
......
I am new user of Bowtie,and can't figure out the format of sam file. If I can get the read start and stop position, I can write a simple script to check if the reads locate at the ranges in the LIST file. Could anyone help? Certainly if there exist available tools, that will simplify the job. Thanks in advance!
2 answers
http://samtools.github.io/hts-specs/SAMv1.pdf
To be succinct - look at column 3 and 4 for the position information.
•
written
by
Brian Bushnell
It's not in the required columns. You have to calculate it by parsing the cigar string for insertion and deletion symbols.
If you map using BBMap, you can use the stoptag flag to add the stop position to each read (prefixed by YS:i:).
•
written
by
Brian Bushnell
Hi Brian, thanks for helpful answer. Can I ask one more question: I noticed Cigar strings have: M, I, D, P, N. To my understanding M means match, all the others are mismatches, INDELs. So I just need to make sum of the digits in front of the letters other than M, then plus the length of read and start position, this gives me the end position. Am I right? Thanks!
no. you need no skip somes bases with the clipping operator 'S' and 'H'. 'M' means that bases are aligned but there could be some mismatches (A vs G, G vs G...).
Just change the format of your list file to match BED format (this is untested):
sed -e 's/\:/\t/' list.txt | awk 'BEGIN{OFS="\t"}{print $1, $3-1, $4,".",".",$2}' > list.bed
Convert your SAM file to BAM, sort and index it with samtools and then just:
samtools view -b -L list.bed sorted.bam > filtered.bam
filtered.bam now contains all of the alignments overlapping the regions in your list.
BTW, don't use bowtie to map possibly spliced reads to the genome, the results won't be very good.
•
written
by
Devon Ryan
Log in to answer this question.
-
Summing Chromosome Sizes
written by selplat21 •
-
Calculate gene % covered at depth X
written by ellieuk •
-
How Many Reads In A Bam File Are Aligned To a Specific Region
written by h.rokny •
-
How do I replace a value used in script with a range of values in a file?
written by josephmccarter1
-
How to convert genomic position to transcript position?
written by biolab
-
bowtie problem: cannot locate bowtie index
written by biolab
-
how to produce mapping graph based on blat result?
written by biolab
-
Question about bowtie
written by biolab
-
Setting elementMetadata in genomicranges list objects
written by reviewer3 •
-
How Tophat Treats Additional Contigs
written by GR
