read counts without gff3 file

[cyendrek]

I am working with a species that does not have a genome sequenced. There is a partial unigene/EST database that I have used to align some Illumina RNA-Seq reads with tophat2/bowtie2. My problem is with obtaining read counts from the .bam files. Normally, I use HT-Seq but in this case I don't have a gff3 file. Our core bioinformatics director said to use custom perl scripts but that is beyond my expertise. Can anyone help?
Thanks,
Craig

[chadn737]

Why not do a de novo assembly?

[jgibbons1]

@cyendrek I've had the same issue before. Rather than using bowtie2 (which I typically use too) try using seqmap and rseq. Both programs are in the rseq package http://www-personal.umich.edu/~jianghui/rseq/

You first need to map the reads to the reference (or in your case unigene/EST) using seqmap then you can generate RPKM and read count calculations per gene/EST using rseq. I've used this pipeline quite a bit in the past so let me know if you have any problems.

Here's an example command line of seqmap allowing 2 mismatches:

[user]$ seqmap 2 ReadFile.fasta Reference.fasta Output.seqmap /eland:3

Here's an example command line of rseq assuming read length is 50 bp:

[user]$ rseq comp_exp -r 50 Reference.fasta Output.seqmap

This will create a file with the "comp_exp" extension that has the number of mapped reads, number of uniquely mapped reads and rpkm values (among other stats).

A few words of wisdom, your reads must be in fasta format, so convert fastq to fasta (I use the fastxtoolkit for this). Also, seqmap uses ALOT of memory, so I usually break my read file up into batches of 5-10 million reads. I then mapped these independently against the reference, merge the output files then run rseq on the merged output.

Good luck!