Question: Extracting nucleotide immediately prior to mapped read
1
gravatar for graeme.thorn
4.6 years ago by
graeme.thorn50
London, United Kingdom
graeme.thorn50 wrote:

Hi,

I have mapped (paired-end) RNA-seq reads (so BAM file information) and I want to extract the nucleotide in the genome immediately prior to the (forward) mapping (essentially to test whether a RNAse digest has worked correctly). 

As it is in BAM format, then extracting the position of the read (chromosome and location) is straightforward (just samtools view it, and process the line output for the relevant field), but I was wondering if there was a quicker way than brute-force (i.e. extracting, sorting, then reading from the genome FASTA file).

rna-seq sequence • 1.4k views
ADD COMMENTlink modified 4.6 years ago by Brian Bushnell17k • written 4.6 years ago by graeme.thorn50
1

Your time spent properly processing the SAM format will probably outweigh indexed FASTA access overhead.

ADD REPLYlink modified 7 months ago by RamRS28k • written 4.6 years ago by Matt Shirley9.4k
4
gravatar for Pierre Lindenbaum
4.6 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum129k wrote:

Using my tool samslop and bioalcidae:

$ java -jar dist/samslop.jar -c -m 1 -M 1 -r ref.fa  S1.bam |\
java -jar dist/bioalcidae.jar -F SAM -e  'while(iter.hasNext()) {var read=iter.next(); if(read.getReadUnmappedFlag()) continue; out.println(read.getReferenceName()+":"+read.getAlignmentStart()+"-"+read.getAlignmentEnd()+" "+read.getReadString().substr(0,1)+"/"+read.getReadString().substr(read.getReadLength()-1));}' | uniq | head

rotavirus:1-71 G/A
rotavirus:1-72 G/A
rotavirus:1-53 G/A
rotavirus:1-72 G/A
rotavirus:1-65 G/A
rotavirus:1-68 G/A
rotavirus:2-73 G/T
rotavirus:3-74 C/A
rotavirus:3-69 C/T
rotavirus:3-74 C/A
ADD COMMENTlink modified 7 months ago by RamRS28k • written 4.6 years ago by Pierre Lindenbaum129k
2
gravatar for Brian Bushnell
4.6 years ago by
Walnut Creek, USA
Brian Bushnell17k wrote:

I think the easiest approach would be to just add an "N" to both ends of the reads before you map them. Then you can see the reference base in the MD tag.

ADD COMMENTlink modified 7 months ago by RamRS28k • written 4.6 years ago by Brian Bushnell17k
2
gravatar for Matt Shirley
4.6 years ago by
Matt Shirley9.4k
Cambridge, MA
Matt Shirley9.4k wrote:

Here is some skeleton code that can help you.

Depending on the question you are asking, you might find the Counter in the collections module useful for summarizing the results.

ADD COMMENTlink modified 7 months ago by RamRS28k • written 4.6 years ago by Matt Shirley9.4k
0
gravatar for Sean Davis
4.6 years ago by
Sean Davis26k
National Institutes of Health, Bethesda, MD
Sean Davis26k wrote:

At least R, python, and perl have indexed fasta capabilities that allow you to read a "slice" from a fasta file or files. You could also look at using the Biostrings package in Bioconductor. Finally, if you have the UCSC kent tools, you can convert your fasta to .2bit format and then use the kent tools to very quickly get bases from the genome.

ADD COMMENTlink modified 7 months ago by RamRS28k • written 4.6 years ago by Sean Davis26k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 725 users visited in the last hour