SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to extract mapped and unmapped raw reads from bwa's sam file ? vaibhavvsk Bioinformatics 11 02-07-2013 10:01 AM
BWA rescue of multi-mapping or unmapped reads kjlee Bioinformatics 6 07-13-2012 11:29 PM
BWA output bitwise flag for mapped/unmapped reads wenhuang Bioinformatics 1 08-29-2011 04:54 PM
BWA, mostly unmapped reads Michael.James.Clark Bioinformatics 4 03-03-2011 03:20 PM
What are the unmapped reads beelu Illumina/Solexa 1 09-09-2010 06:18 AM

Reply
 
Thread Tools
Old 06-17-2010, 01:05 AM   #1
Adamo
Member
 
Location: Paris

Join Date: Jun 2010
Posts: 28
Default BWA - unmapped reads

Hello everyone,

I'm trying to combine the use of bwa and blat but I have a few difficulties to extract the unmapped reads from bwa.
It appears they are in the first .sam file created using bwa bwasw.
I used the command below to filter the SAM file and extract unaligned reads into the file unalntest.bam:

Quote:
$ samtools view -f 4 /home/biopuce11/Documents/Adam/bwa-0.5.8/aln/aln.sam -t -b -S -o /home/biopuce11/Documents/Adam/bwa-0.5.8/aln/unalntest.bam
The problem is that when I flagstat it to see how many reads it contains, the soft tells me that there is no read in it.

By the way, I have this message using "samtools view..." (above):
Quote:
[bam_header_read] EOF marker is absent.
[bam_flagstat_core] Truncated file? Continue anyway.
Could someone help me, please? I'm turning mad!

Thanks.
Adamo is offline   Reply With Quote
Old 06-17-2010, 07:42 AM   #2
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

Quote:
Originally Posted by Adamo View Post
Hello everyone,

I'm trying to combine the use of bwa and blat but I have a few difficulties to extract the unmapped reads from bwa.
It appears they are in the first .sam file created using bwa bwasw.
I used the command below to filter the SAM file and extract unaligned reads into the file unalntest.bam:

The problem is that when I flagstat it to see how many reads it contains, the soft tells me that there is no read in it.
What about:

Code:
$ samtools view -b -S -f 4 input.sam > output.bam
Quote:
By the way, I have this message using "samtools view..." (above):
Could someone help me, please? I'm turning mad!
Read this:

http://seqanswers.com/forums/showpos...72&postcount=2
__________________
-drd
drio is offline   Reply With Quote
Old 06-21-2010, 02:46 AM   #3
Adamo
Member
 
Location: Paris

Join Date: Jun 2010
Posts: 28
Default

I have resolved the problem concerning the error message, thanks.
But for the other one, I'm still try to figure out where the problem is.
I've tried your command (btw, I didn't paste perfectly mine ("-f 4" is missing...) in my first post), but the problem remains...
Any idea?

Last edited by Adamo; 06-21-2010 at 07:49 AM.
Adamo is offline   Reply With Quote
Old 06-21-2010, 07:50 AM   #4
tcezard
Member
 
Location: Edinburgh

Join Date: Dec 2008
Posts: 13
Default

if you type
Code:
samtools view -?
you'll see that to convert sam to bam you'll need the reference fasta file

Code:
2. SAM->BAM conversion: `samtools view -bT ref.fa in.sam.gz'.
your command should probabaly be
Quote:
samtools view -bT ref.fa -f 4 input.sam > output.bam
or
Quote:
samtools view -bt ref.fa.fai -f 4 input.sam > output.bam
tcezard is offline   Reply With Quote
Old 06-22-2010, 01:09 AM   #5
Adamo
Member
 
Location: Paris

Join Date: Jun 2010
Posts: 28
Default

I think the command's ok but still, I have no read in the output file. I don't get it!
I have only about 2% of mapped reads in my .sam file according to flagstat...
I've tried both commands you suggest.
Is there something I have forgotten?
Adamo is offline   Reply With Quote
Old 06-22-2010, 01:45 AM   #6
tcezard
Member
 
Location: Edinburgh

Join Date: Dec 2008
Posts: 13
Default

It should work ...
Can you paste your commands and the first few lines of your sam file
tcezard is offline   Reply With Quote
Old 06-22-2010, 01:58 AM   #7
Adamo
Member
 
Location: Paris

Join Date: Jun 2010
Posts: 28
Default

The first lines of my sam file are:

Quote:
@SQ SN:Clostridium_botulinum_chromosome LN:3886916
@SQ SN:Clostridium_botulinum_plasmid LN:16344
SRR037154.1271 0 Clostridium_botulinum_chromosome 9824 0 35S229M3S * 0 0 CAGCTTAGGCTGCTAAAGCGCACGCAGGCGGTCTGTTAAGTCAGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTTGAAACTGGCAGGCTTGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGATAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACGCCGTAAACGATGTCG FFFFFFFFFFFFFI@@@GIIIIIIIIIIIIIIIIIIIIIIIIIIIFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFG;55555CGEGGFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFBBBBBGGFFFGGGGFFFFFFFFFFFFFFFFFFFFFFFFGA@@EE@7777@@@EEEEEEEEBBBGGIIIAAAAAAAAAAAAAA@@@@@@@@ AS:i:73 XS:i:73 XF:i:1 XE:i:0 XN:i:0
My command is:

Quote:
$ samtools view -bT /home/db.fa -f 4 /home/alncomplet.sam > /home/unalncomplet.bam
I also tried -bt instead of -bT (replacing db.fa by db.fa.fai (I indexed my database with samtools faidx)) but the result is just the same...
Adamo is offline   Reply With Quote
Old 06-22-2010, 02:12 AM   #8
tcezard
Member
 
Location: Edinburgh

Join Date: Dec 2008
Posts: 13
Default

I'm wondering if you have unmapped read in your sam file to beggin with.

try
Code:
samtools view -bT ref.fa sam_file | samtools sort - bam_file 
(with no .bam extension)
then
 samtools flagstat bam_file
This will tell you what you have in your bam file.
Then you can extract what you want from the bam file with
Code:
samtools view bam_file -f 4 | less
to see the unmapped reads
or which ever flag you want

but if your original file does not contain any unmapped read then this is not going to help
tcezard is offline   Reply With Quote
Old 06-22-2010, 02:31 AM   #9
Adamo
Member
 
Location: Paris

Join Date: Jun 2010
Posts: 28
Default

Well, I'm a little confused.
My sam file contains only the reads mapped. The thing is that I read on a topic a post from lh3 (who know what he speaks about, I think!) telling someone that all the informations concerning unmapped reads where in the sam file. I assumed he talked about the one created using:
Quote:
$bwa bwasw ...
That's why I'm trying to extract them from this sam file.
Otherwise, which one was lh3 speaking about? Or which one are you speaking about?
Adamo is offline   Reply With Quote
Old 06-22-2010, 02:48 AM   #10
tcezard
Member
 
Location: Edinburgh

Join Date: Dec 2008
Posts: 13
Default

I'm sure Heng li knows what he's talking about so you should check your commands and log messages, parameters used, maybe rerun the alignment command (and keep the command you're running).
If it really does not output the unmapped read then send a email to the bwa mailing list.
and try to be a bit more details on what you've done and what you have otherwise you won't get any help.
tcezard is offline   Reply With Quote
Old 06-22-2010, 02:54 AM   #11
Adamo
Member
 
Location: Paris

Join Date: Jun 2010
Posts: 28
Default

I'll try the mailing list, thanks for the tip.
I've already run bwa several times, with default parameters and followed the steps suggested by the readme. Nothing more, nothing less...

I'm quite a newbie here, excuse me for not being precise enough.

Last edited by Adamo; 06-22-2010 at 02:59 AM.
Adamo is offline   Reply With Quote
Old 06-22-2010, 08:49 AM   #12
Adamo
Member
 
Location: Paris

Join Date: Jun 2010
Posts: 28
Default

After further researches, it seems like bwasw doesn't include any unmapped read in the output file, unlike bwa aln.
I think I might have to write a script which could handle such a job.

Thanks for your help anyway!
Adamo is offline   Reply With Quote
Old 10-05-2010, 12:01 PM   #13
freeseek
Junior Member
 
Location: Cambridge MA

Join Date: Feb 2010
Posts: 8
Default

I confirm that I find that 'bwa samse/sampe' does output the unaligned reads in the sam file while 'bwa bwasw' does not. This difference in behavior is quite odd to me. Is there a reason for it? This is not even mentioned in the bwa manual. But then, bwa is still at version 0.5.8a, so I should not expect too much. I would be really interested in seeing what is being left after aligning all the reads with samse and bwasw, using blat or blast.
freeseek is offline   Reply With Quote
Old 10-07-2010, 03:48 PM   #14
freeseek
Junior Member
 
Location: Cambridge MA

Join Date: Feb 2010
Posts: 8
Default bwasw patch

This simple patch:
Code:
--- bwtsw2_aux.c
+++ bwtsw2_aux.c
@@ -520,8 +520,10 @@
 			bsw2_resolve_query_overlaps(b[0], opt.mask_level);
 		} else b[1] = 0;
 		// generate CIGAR and print SAM
-		gen_cigar(&opt, l, seq, pac, b[0]);
-		print_hits(bns, &opt, p, b[0]);
+		if(b[0]->n > 0) {
+		    gen_cigar(&opt, l, seq, pac, b[0]);
+		    print_hits(bns, &opt, p, b[0]);
+		}
 		// free
 		free(seq[0]);
 		bsw2_destroy(b[0]);
@@ -583,6 +585,7 @@
 	for (i = 0; i < _seq->n; ++i) {
 		bsw2seq1_t *p = _seq->seq + i;
 		if (p->sam) printf("%s", p->sam);
+		else printf("%s\t4\t*\t0\t0\t*\t*\t0\t0\t%s\t%s\n", p->name, p->seq, p->qual);
 		free(p->name); free(p->seq); free(p->qual); free(p->sam);
 		p->tid = -1; p->l = 0;
 		p->name = p->seq = p->qual = p->sam = 0;
will make bwasw spit out also the unmapped reads, though it only works with .fastq files, so use it carefully. To use it, create a file contained the above code, like "bwtsw2_aux.diff", enter the bwa source directory, and issue the command "patch < bwtsw2_aux.diff". (PS You need to be careful to not make the tabs become spaces, or the patch will not work)

Last edited by freeseek; 10-07-2010 at 04:09 PM.
freeseek is offline   Reply With Quote
Old 10-08-2010, 04:41 PM   #15
Pepe
Member
 
Location: Germany

Join Date: Mar 2009
Posts: 28
Default

This thread contains a perl script to get reads that are in a .fq file but not in a sam file. It was written with TopHat in mind but it may work for you too.
I'm not sure if it helps.
http://seqanswers.com/forums/showthr...eferrerid=2547
Pepe is offline   Reply With Quote
Old 10-08-2010, 07:46 PM   #16
freeseek
Junior Member
 
Location: Cambridge MA

Join Date: Feb 2010
Posts: 8
Default

@Pepe I am sure you can make a script to perform the same task, though it is conceptually easier to avoid dropping the unmapped reads while you are processing them with 'bwa bwasw' one by one rather than figuring out later what you dropped. I also wanted to mention that occasionally 'bwa bwasw' reports multiple alignments in the sam file for the same read. The only way to modify that behavior at the moment is to get into the source code.
freeseek is offline   Reply With Quote
Old 02-21-2012, 04:12 PM   #17
Bgansw
Member
 
Location: Sydney, Australia

Join Date: Nov 2011
Posts: 24
Default Unmapped reads bwasw + patch

Quote:
Originally Posted by freeseek View Post
This simple patch:
Code:
--- bwtsw2_aux.c
+++ bwtsw2_aux.c
@@ -520,8 +520,10 @@
 			bsw2_resolve_query_overlaps(b[0], opt.mask_level);
 		} else b[1] = 0;
 		// generate CIGAR and print SAM
-		gen_cigar(&opt, l, seq, pac, b[0]);
-		print_hits(bns, &opt, p, b[0]);
+		if(b[0]->n > 0) {
+		    gen_cigar(&opt, l, seq, pac, b[0]);
+		    print_hits(bns, &opt, p, b[0]);
+		}
 		// free
 		free(seq[0]);
 		bsw2_destroy(b[0]);
@@ -583,6 +585,7 @@
 	for (i = 0; i < _seq->n; ++i) {
 		bsw2seq1_t *p = _seq->seq + i;
 		if (p->sam) printf("%s", p->sam);
+		else printf("%s\t4\t*\t0\t0\t*\t*\t0\t0\t%s\t%s\n", p->name, p->seq, p->qual);
 		free(p->name); free(p->seq); free(p->qual); free(p->sam);
 		p->tid = -1; p->l = 0;
 		p->name = p->seq = p->qual = p->sam = 0;
will make bwasw spit out also the unmapped reads, though it only works with .fastq files, so use it carefully. To use it, create a file contained the above code, like "bwtsw2_aux.diff", enter the bwa source directory, and issue the command "patch < bwtsw2_aux.diff". (PS You need to be careful to not make the tabs become spaces, or the patch will not work)
Hi Freeseek,



I'm trying to align my sample illumina bacterial genome to an already published reference bacterial genome. The aim is to find out what genes are present in my sample genome, but missing in the reference. The plan of action was to run bwasw, and then pull out the sequence in the sample genome that do not align to the reference, then look at what these sequences are.

I aligned my query genome to the ref using bwasw.

The output was in sam format. I converted that to bam format using

samtools view -bS aln.sam > aln.bam

Then I ran the command lines quoted above.

When I headed my mapped.sam file, it appears to have worked, however, my unmapped.sam file is empty,the size is zero.
This is consistent for two of my sample genomes. Reading this thread, it seems that bwasw does not give out unmapped reads, thus I was trying your patch.

This is the message I'm getting:
bgansw:~/bwa$ patch < bwtsw2_aux.diff
can't find file to patch at input line 3
Perhaps you should have used the -p or --strip option?
The text leading up to this was:
--------------------------
|--- bwtsw2_aux.c
|+++ bwtsw2_aux.c
--------------------------
File to patch:

I'm probably going to ask a very silly question. I'm a newbie, so am not good with scripts. What should I enter in file to patch? Do I have to make changes or input a file name in the script? I've got paired end illumina hiseq data.

Any help would be greatly appreciated.

bgansw
Bgansw is offline   Reply With Quote
Old 02-21-2012, 04:21 PM   #18
Bgansw
Member
 
Location: Sydney, Australia

Join Date: Nov 2011
Posts: 24
Default

Hi Freeseek,

I inputted my faq file in the file to patch. That gave me .fastq.orig file and .fastq.rej

Copy pasting the .orig file

bgansw:~/bwa$ head 6_combined_paired.faq.orig
@GAPC_0027_FC:7:1:4477:1032#TGACCA/1
AGTATGGGGGAAATGCAAGCAAGCATTNATGCTAGTGGTTGTCAAATAGTTACCGTAGCTGTGCGCA
+
hhhhhhhhhhdhhghhhhhhhheddeeBededaad^dfdfffhhhhhhfhhfhhfceffdfbhhhad
@GAPC_0027_FC:7:1:4477:1032#TGACCA/2
TGCTTCTAAGGTACCAATGGNNNCAGGTAACAGGTATTTAGGATCGGGTATGACCTNNNNTTTAACA
+
hhhhhhhhhhhefhfhf_]_BBB`]aa^acdcff`fffffffhghhgfafccY][ZBBBBW\\[[[[
@GAPC_0027_FC:7:1:4986:1032#TGACCA/1
GAAAATCCGAAAGAACCGGCAATGATTNTTTCTACGTTCCGTAAAGAAAGCAATGGTGATAAGCCTG


Am copy pasting the .fastq.rej file:


bgansw:~/bwa$ cat 6_combined_paired.faq.rej
***************
*** 520,527 ****
bsw2_resolve_query_overlaps(b[0], opt.mask_level);
} else b[1] = 0;
// generate CIGAR and print SAM
- gen_cigar(&opt, l, seq, pac, b[0]);
- print_hits(bns, &opt, p, b[0]);
// free
free(seq[0]);
bsw2_destroy(b[0]);
--- 520,529 ----
bsw2_resolve_query_overlaps(b[0], opt.mask_level);
} else b[1] = 0;
// generate CIGAR and print SAM
+ if(b[0]->n > 0) {
+ gen_cigar(&opt, l, seq, pac, b[0]);
+ print_hits(bns, &opt, p, b[0]);
+ }
// free
free(seq[0]);
bsw2_destroy(b[0]);
***************
*** 583,588 ****
for (i = 0; i < _seq->n; ++i) {
bsw2seq1_t *p = _seq->seq + i;
if (p->sam) printf("%s", p->sam);
free(p->name); free(p->seq); free(p->qual); free(p->sam);
p->tid = -1; p->l = 0;
p->name = p->seq = p->qual = p->sam = 0;
--- 585,591 ----
for (i = 0; i < _seq->n; ++i) {
bsw2seq1_t *p = _seq->seq + i;
if (p->sam) printf("%s", p->sam);
+ else printf("%s\t4\t*\t0\t0\t*\t*\t0\t0\t%s\t%s\n", p->name, p->seq, p->qual);
free(p->name); free(p->seq); free(p->qual); free(p->sam);
p->tid = -1; p->l = 0;
p->name = p->seq = p->qual = p->sam = 0;

My apologies, since I dont know if I'm on the right track or not. Any help would be greatly appreciated.

Regards

bgansw
Bgansw is offline   Reply With Quote
Old 02-21-2012, 04:37 PM   #19
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

If you have Illumina reads, just use bwa aln. It doesn't much matter how long your reads are, you are only aligning to a bacterial genome, how long can it take?

The post you are citing is a year old, do you know that it's suppsoed to work for whatever version of bwa you have?
swbarnes2 is offline   Reply With Quote
Old 02-21-2012, 05:43 PM   #20
Bgansw
Member
 
Location: Sydney, Australia

Join Date: Nov 2011
Posts: 24
Default

Hi swbarnes2,

I ran bwa using the aln parameter, instead of the bwasw, but am not sure that its worked. The output .sam file is only 4K, while the bwasw output file was 796K. Additionally, the aln.sam file looks nothing like the bwasw.sam file. It doesnt contain reads, looks more like a .bam file. The aln.sam file also doesnt contain @SQ headers, while the bwasw file contained them, so I'm unable to convert the sam file to .bam.

I tried

samtools faidx ref.fa (but I get an error - [fai_build_core] different line length in sequence 'NZ_A01000034'.
Segmentation fault

So am stuck...

I would really appreciate any help.
Bgansw is offline   Reply With Quote
Reply

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 09:00 PM.


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