Extracting unmapped reads to draft genome

Question

[Bioinfonext]

I mapped de novo assembled transcript to genome using GMAP but still large no. of transcripts showing similarity with draft genome sequences. My aim is it extract unique sequence which do not mapping to draft genome?

I have 32 paired end RNA seq libraries. Can I take all raw reads R1 and R2 reads in two separated files and is it possible to map those raw reads to Draft genome and to extract unmapped R1 and R2 reads in two separete files and followed by de novo assembly?

Please suggest which software should I use to complete this analysis?

Answer 1

[Brian Bushnell]

There are many ways to do this. I suspect you would want to assemble reads in which neither read in a pair maps to the reference, which you can do like this with BBMap:

bbmap.sh in=r1.fq in2=r2.fq ref=reference.fasta outm=mapped.sam outu=unmapped1.fq outu2=unmapped2.fq

Now unmapped reads are in unmapped1.fq and unmapped2.fq. Alternately, you could use reformat.sh and repair.sh on an existing sam file:

reformat.sh in=all.sam out=unmapped.fq unmappedonly
repair.sh in=unmapped.fq  out=r1.fq out2=r2.fq outs=singleton.fq

This will give you both unmapped pairs and unmapped singletons.

[Bioinfonext]

Hi Brian,

Thanks for helping me.

I have downloaded BBmap. so before running cammand Do I need to do indexing of reference genome.

[Brian Bushnell]

The command I gave will do indexing first, then map. Alternately you could do it in two steps, like this:

bbmap.sh ref=reference.fasta
bbmap.sh in=r1.fq in2=r2.fq outm=mapped.sam outu=unmapped1.fq outu2=unmapped2.fq

Either way gives the same result.

[Bioinfonext]

Hi Brian

Thanks

I got result as below, I have some queries:

1) ) During mapping how many mismatch it allowed, is there any option by which we can adjust the mismatch. 2) How it map splicing variant?

Pairing data:           pct reads       num reads       pct bases          num bases

mated pairs:             79.8502%         8232674        79.8502%         2058168500
bad pairs:                1.2823%          132203         1.2823%           33050750
insert size avg:          264.21


Read 1 data:            pct reads       num reads       pct bases          num bases

mapped:                  83.9105%         8651295        83.9105%         1081411875
unambiguous:             71.2570%         7346701        71.2570%          918337625
ambiguous:               12.6535%         1304594        12.6535%          163074250
low-Q discards:           0.0165%            1702         0.0165%             212750

perfect best site:       27.3360%         2818387        27.3360%          352298375
semiperfect site:        27.4375%         2828843        27.4375%          353605375
rescued:                  8.0776%          832816

Match Rate:                   NA               NA        66.6516%         1057771997
Error Rate:              67.3220%         5824231        33.3197%          528789715
Sub Rate:                54.7312%         4734961         1.2152%           19285478
Del Rate:                30.2153%         2614018        31.8588%          505604634
Ins Rate:                 8.2683%          715313         0.2457%            3899603
N Rate:                   0.3170%           27421         0.0287%             454797

[GenoMax]

BBMap does not have a specific mismatch number. To quote @Brian from a recent answer:

The default is roughly 76% identity. You can adjust this with the "minid" flag (e.g. "minid=0.80" for 80% identity.) If you want to restrict alignments to a maximum number of substitutions, you can use "subfilter"; e.g., "subfilter=5" will discard alignments with more than 5 substitutions.

Splice variants would be mappable based on the setting used for (maxindel and intronlen).

[Brian Bushnell]

Splice variants would be mappable based on the setting used for (maxindel and intronlen).

Hmmm, that was my mistake, for some reason I neglected to mention maxindel. When mapping RNA-seq data to a genome, maxindel is a useful flag to adjust; the default (16000) is fine for fungi and many plants, which have short introns, but for things like mammals which have long introns, I suggest setting adding the flags "maxindel=400000 intronlen=10". That allows mapping across introns of up to around 400 kbp or so. "intronlen" is normally unnecessary but may affect some downstream programs.

[Bioinfonext]

Thanks,

I want to map raw reads to CDS instead of genome than also should I run this with default setting?

[Brian Bushnell]

Yes, for mapping RNA-seq reads to transcripts default settings are fine. You may want to add "ambig=all" because some transcripts have multiple isoforms.