This is a test version of Biostars. For the public version, visit https://www.biostars.org.
finding reads from fastq files corresponds to particular transcript from fasta file

Hello all,

I have performed mRNA seq for few samples and analysed them by using edgeR pipeline. Here I am looking for one particular trancript's expression profile. After analysis I find out that it is not expressing at all in any sample. This transcript does not have single read matching from the data set after featureCount analysis.

But when I did tblastx with the trinity denovo fasta file, this transcript is present there. which means there is a transcript present matching exactly with the query transcript that I am looking for. This means there are some reads present corresponds to my query transcript. but this is not being detected by the edgeR pipeline.

Can somebody please tell me how can I retrieve all the fastq reads corresponding to the particular transcript from fastq file?

Thank you in advance.

rna-seq
• 4,479 views • link
updated by Ram Ram • written by amoltej amoltej

3 answers

You could try a tool like RapMap

You should be able to use your assembled transcripts and map your reads to your assembly. This would answer the question of which reads were used in assembly that were not counted for differential expression analysis.

• 0 views • link
updated by Ram Ram • written by Matt Shirley Matt Shirley

Yup; RapMap should handle this use case well. You could simply map all of your reads to your assembled transcriptome and then look for ones that map to your transcript of interest, or, if you're only interested in knowing what maps to this transcript, simply build the index on that transcript alone.

• 0 views • link
updated by Ram Ram • written by Rob Rob

You can use the bbmap toolkit as such:

bbduk.sh -Xmx128544m in=YOUR_FASTQ_file.fastq ref=Your_transcript.fasta k=13 outm=OUTPUT_FASTQ_file_with_transcript_reads_only.fastq

This will generated a fastq file with reads matching your transcript as you requested. The fasta file you supply should have only your transcript of interest for this to work as you want it.

• 0 views • link
updated by Ram Ram • written by Adrian Pelin Adrian Pelin

Additionally, you can separate the reads into one file per fasta sequence in a single command with Seal (also part of BBMap):

seal.sh ref=transcripts.fasta in=reads.fastq pattern=out_%.fastq
• 0 views • link
updated by Ram Ram • written by Brian Bushnell Brian Bushnell

You could make a fasta file with your transcripts and then align the reads to these transcripts. Or, if you have a genome and transcript annotations, you could use tophat for aligning reads to exons.

If you are only interested in one particular transcript, you can rely on blast as you already have done.

HTH
Federico

• 0 views • link
updated by Ram Ram • written by abascalfederico abascalfederico

Log in to answer this question.

More posts like this