Draft genome scaffolding with RNAseq paired-end reads

[JPZ]

Hello all,

I used Tophat to map 100bp, PE Illumina transcriptome reads to a draft genome (133062 contigs).
Our main goal was SNP mining, but I have been suggested the reads could also be used for scaffolding.

I have no experience in genome assembly and scaffolding, but I assume that if I can find read pairs where the 2 reads are mapped to different genomic contigs, the 2 genomic contigs could then be connected.

How can I search the BAM alignments for such read pairs?

Alternatively I could use an assembler that can combine different types of reads such as Mira, but I thought it would take longer, and the genomic reads are not available anyway.

Thank you!

[Brian Bushnell]

Look for reads where the "rname" field and "rnext" field are different (and rnext is not "=" or "*"); those have the reads mapped to different contigs.

[JPZ]

Thank you for the info Brian,
good starting point, saved me a lot of reading and guesswork.

I also thought of filtering for MAPQ = 50 (should be uniquely mapped reads)
and properly paired reads (FLAG = 83|99|147|163)

The following command should then extract the alignments of interest:

samtools view -q 50 accepted_hits.bam |gawk '($2 == 83 || $2 == 99 || $2 == 147 || $2 == 163) && $7 !~/[*=]/ {print $3, $7}' > output

And thus obtain a list of joined contigs.
However, while it is possible to determine which contigs are joined, I assume the lenght of N bases padding between them cannot.
Not only the region may not be transcribed, but the insert size for paired reads that have a mate in a different contig appears to be always 0 (at least that is what Tablet shows).

Or do I have other options I'm unaware of?

[Brian Bushnell]

The insert size of reads mapped to different contigs is unknown. Scaffolding tools can use the distribution of insert sizes of pairs on the same contig, or user-supplied insert size numbers, to determine how many Ns to pad.

This might be easier if you just use a standalone scaffolding tool. There are various out there, but I don't have a recommendation. Here's a paper comparing some of them:

http://genomebiology.com/2014/15/3/R42

[JPZ]

Thank you again for the input,
I'll check the paper and see if using the above filtered alignments can work.

[boetsie]

Maybe this program will help:

https://github.com/svm-zhang/AGOUTI

I have not used it myself though.

[bioBob]

L-rna-scaffolder may also help.

I have used it with varying success.

[JPZ]

Thank you all for your answers,
I'll try some of the suggested tools, more likely those that do not have too many dependencies..