Tweet
filtering out reads from abundant transcripts before using velevt
Hi everybody,
My situation is as following: We did a whole Solexa flow cell RNA-seq from a non model organism (so no refrence genome) and got 225 Million paired end reads (so a 112 Million reads for each direction). Additionally we did a 454 sequencing (also mRNA) of the same organism and got ~8000 contigs from that (many of these are very highly abundant in the tissue I am looking at). Additionally I have to say that a substantial amount of reads matches to ribosomal RNA (you can not get rid of these in the experimental protocol is a problem of this organism).
From this I want to assemble the transcriptome using velvet (or abyss).
As discussed here this many reads are difficult to handle in velvet (alone hashing the reads in velvet with 16 CPU's takes forever). The other problem is the very non-random coverage of the transcriptome (as discussed here as well) which is in my case pretty extreme (the tissue we look at is highly specialized). So I was thinking of the following strategy: Aligning all the short reads to the ribosomal RNA and the abundant transcripts from the 454 sequences with bowtie and then filtering out these reads before I start using velvet. In this way I would reduce complexity for velvet and at the same time get a more even coverage (yes I do not assemble or correct the abundant transcripts I have but this is not my focus right now). So far so good but my problem is that I don't know exactly how to do that. Bowtie gives me a SAM file which I would now have to compare with the Solexa reads and take out the ones that gave a match in bowtie. Since I want to use paired ends to resolve errors in velvet I need to keep the right order in the file for velvet (I used the shufflesequence algorythm coming with the velvet package) and also if a read aligns in bowtie I have to filter out this one and the corresponding mate pair. Does anybody know of script of some kind or can give me any advice on this. I found in this forum a post describing the filtering of ribosomal reads with the grep -v command but I am not sure this will work with amount of data.
I would appreciate any help you could give me with this.
Thanks for reading this long post.
Marco
My situation is as following: We did a whole Solexa flow cell RNA-seq from a non model organism (so no refrence genome) and got 225 Million paired end reads (so a 112 Million reads for each direction). Additionally we did a 454 sequencing (also mRNA) of the same organism and got ~8000 contigs from that (many of these are very highly abundant in the tissue I am looking at). Additionally I have to say that a substantial amount of reads matches to ribosomal RNA (you can not get rid of these in the experimental protocol is a problem of this organism).
From this I want to assemble the transcriptome using velvet (or abyss).
As discussed here this many reads are difficult to handle in velvet (alone hashing the reads in velvet with 16 CPU's takes forever). The other problem is the very non-random coverage of the transcriptome (as discussed here as well) which is in my case pretty extreme (the tissue we look at is highly specialized). So I was thinking of the following strategy: Aligning all the short reads to the ribosomal RNA and the abundant transcripts from the 454 sequences with bowtie and then filtering out these reads before I start using velvet. In this way I would reduce complexity for velvet and at the same time get a more even coverage (yes I do not assemble or correct the abundant transcripts I have but this is not my focus right now). So far so good but my problem is that I don't know exactly how to do that. Bowtie gives me a SAM file which I would now have to compare with the Solexa reads and take out the ones that gave a match in bowtie. Since I want to use paired ends to resolve errors in velvet I need to keep the right order in the file for velvet (I used the shufflesequence algorythm coming with the velvet package) and also if a read aligns in bowtie I have to filter out this one and the corresponding mate pair. Does anybody know of script of some kind or can give me any advice on this. I found in this forum a post describing the filtering of ribosomal reads with the grep -v command but I am not sure this will work with amount of data.
I would appreciate any help you could give me with this.
Thanks for reading this long post.
Marco
Comment