SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Extracting unmapped pairs of mapped reads from a particular chromosomal position Xespok Bioinformatics 2 01-07-2014 10:48 PM
extracting unmapped reads from BAM using Samtools? Lspoor Bioinformatics 17 08-25-2013 01:22 AM
Picard, extracting unmapped reads and converting it to FASTQ gabrieledcjr Bioinformatics 0 05-15-2013 04:57 PM
454: Unmapped contigs after Reference Assembly Alex Clop 454 Pyrosequencing 5 03-28-2011 03:57 AM
mate unmapped and read unmapped rururara Bioinformatics 1 02-25-2011 02:31 AM

Reply
 
Thread Tools
Old 10-13-2015, 12:16 PM   #1
antifolate
Member
 
Location: Chesterfield, MO

Join Date: Aug 2015
Posts: 52
Default Extracting Unmapped Contigs

Hello,

I've assembled some contigs from a sequenced strain and aligned them back to the reference with Mauve. I'm interested in studying those contigs that did not align. How can I extract those contigs?

Mauve can export a file with the positions of all those contigs in the input file so I can probably write a script to get the sequences, but I thought there's probably an easier way. Also, if you know a better software than Mauve, please tell me. My purpose is comparing several de-novo-assembled strains with each other and with the reference genome (~20MB) to see sequence variations.

Thanks!
antifolate is offline   Reply With Quote
Old 10-13-2015, 01:28 PM   #2
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Well, Mauve is usually used for genome-to-genome alignments. Although it can be used to map a draft genome with many contigs to a more finished genome. Depending on the size of your contigs and the reference perhaps a mapper such as BBMap might work better -- you would at least get a BAM file that you could work with.

But since you are already using Mauve, well, it seems simple to me to write a short script via 'cut', 'grep', 'paste', 'diff' etc. to pull out the non-mapping contigs.

I know my response isn't much help. Sorry about that. Just random thoughts at the end of a long day.
westerman is offline   Reply With Quote
Old 10-13-2015, 02:54 PM   #3
antifolate
Member
 
Location: Chesterfield, MO

Join Date: Aug 2015
Posts: 52
Default

I don't know why I went and installed Mauve while I already had the BB package. But now that I tried BBMap, I don't see how a BAM file is better (I'm not even sure this output is BAM). Is there an option in BBMap to retain the contigs that don't get mapped?

I'll be doing this for several assemblies so it'd be easier and faster if there was a simple flag for it, rather than having to massage the data every time with unix tools.
antifolate is offline   Reply With Quote
Old 10-13-2015, 03:07 PM   #4
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

BBMap can separate mapped and unmapped sequences, and keep them in fasta format, like this:

mapPacBio.sh in=contigs.fasta ref=reference.fasta outm=mapped.fasta outu=unmapped.fasta

That has a sequence length limit of 6000bp; longer sequences will be broken into fragments (bbmap.sh has a length limit of 600bp while mapPacBio is 6000bp). On the other hand, BBDuk is alignment-free and can do the same thing using kmer matching, without breaking up the input contigs:

bbduk.sh in=contigs.fa ref=reference.fa outm=matching.fa outu=unmatching.fa k=31 mcf=0.7

"mcf" means "min covered fraction", so this will make input sequences in which at least 70% of their bases are covered by reference kmers go to outm, and the others will go to outu.

Another way of doing it is to shred the input (or reference) contigs and map them to the other, get a bed file of the resulting sam, and extract the covered regions. There are a lot of ways, I guess, depending on your specific goal.

Last edited by Brian Bushnell; 10-13-2015 at 03:09 PM.
Brian Bushnell is offline   Reply With Quote
Old 10-13-2015, 03:52 PM   #5
antifolate
Member
 
Location: Chesterfield, MO

Join Date: Aug 2015
Posts: 52
Default

Thanks Brian! About BBMap, does that mean that I shouldn't be using it (or even mapPacBio) to map contigs with n50=91k into the reference (s. pombe)?

Also could you please please go into the details of your last method?
antifolate is offline   Reply With Quote
Old 10-13-2015, 06:09 PM   #6
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Quote:
Originally Posted by antifolate View Post
Thanks Brian! About BBMap, does that mean that I shouldn't be using it (or even mapPacBio) to map contigs with n50=91k into the reference (s. pombe)?
It really depends on your purpose, but for anything that requires the contigs to remain intact, it's not ideal; Seal (or an aligner for unlimited-length reads) would be better.

Quote:
Also could you please please go into the details of your last method?
For example -

bbduk.sh in=contigs.fa ref=reference.fa out=masked_contigs.fa kmask=N k=31

That would produce a file like the contigs.fa, but with all bases covered by reference kmers masked (converted to N). This is handy because it still works even if both fastas contain a complete, single-contig genome, without shredding anything or changing any coordinates. Filtering methods would not really work well for single-contig genomes.
Brian Bushnell is offline   Reply With Quote
Old 10-14-2015, 04:28 AM   #7
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,016
Default

@antifolate: You should post this question on Mauve list. Aaron will generally address mauve questions in a day or two. It may be possible to extract the info you need from one of the mauve files.
GenoMax is offline   Reply With Quote
Old 10-27-2015, 01:59 PM   #8
antifolate
Member
 
Location: Chesterfield, MO

Join Date: Aug 2015
Posts: 52
Default

@Brian Bushnell Is there a way to make bbduk.sh not count mismatches resulting from Ns towards the 0.7?

Also, what is k?

@GenoMax Can you link me to the Mauve list? I can't find it.
antifolate is offline   Reply With Quote
Old 10-27-2015, 02:08 PM   #9
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Quote:
Originally Posted by antifolate View Post
@Brian Bushnell Is there a way to make bbduk.sh not count mismatches resulting from Ns towards the 0.7?
You cannot do that with the "mcf" (min covered fraction) flag, but you can with the "mkf" (min kmer fraction) flag.

mcf=0.7: reference kmers must cover 70% of the sequence's bases.

mkf=0.7: sequence must share 70% of its valid kmers with the reference.

The distinction is subtle, but mkf takes into account Ns - sequence with Ns do not have valid kmers, so it is ignored with respect to the 70%. mcf does not take Ns into account. Not currently, at least. I might change that in the future.

Quote:
Also, what is k?
"k" is the kmer length. 31 is BBDuk's maximum; you can make it shorter if you want, but specificity starts declining rapidly once you drop below 20.
Brian Bushnell is offline   Reply With Quote
Old 10-27-2015, 02:48 PM   #10
antifolate
Member
 
Location: Chesterfield, MO

Join Date: Aug 2015
Posts: 52
Default

God bless Brian Bushnell! Thanks for the quick reply!
antifolate is offline   Reply With Quote
Old 10-27-2015, 03:03 PM   #11
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,016
Default

Quote:
Originally Posted by antifolate View Post
@GenoMax Can you link me to the Mauve list? I can't find it.
mauve-users at lists.sourceforge.net

You may need to become a member first before you can post.
GenoMax is offline   Reply With Quote
Reply

Tags
alignment, assembly, contigs, mauve

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 07:11 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO