SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to batch download thousands of FASTA files? J:Mo Bioinformatics 6 12-18-2014 02:00 PM
Various methods of creating .sam/.bam files wmyashar Bioinformatics 3 03-14-2013 12:49 PM
Grouping sam or bam files ahmadsam Bioinformatics 2 06-13-2012 04:59 AM
[help] Cuffdiff with BAM or SAM files stavi Illumina/Solexa 2 03-30-2012 05:52 AM
NEw to Chip-seq and have .bam/.sam/.bam.bai files... then what? NGS newbie Bioinformatics 11 05-25-2011 07:48 AM

Reply
 
Thread Tools
Old 06-05-2013, 11:16 AM   #1
jazz
Member
 
Location: us

Join Date: Nov 2008
Posts: 28
Default Split SAM/BAM files into thousands of contigs

Hi,

I am new to bioinformatics and programming, so this might be too easy to tackle for some of you. I am working on a newly sequenced genome with thousands of contigs. I have Illumina sequence data that I have mapped using bowtie and BWA. When I try to import that file into my genome viewer, Geneious, it crashes. So, I want to split my BAM files into many, so the program can handle them easily. My question is if there is an easy way to split bam files by each of the contigs. So, in effect, it should split one BAM file into thousands of mini-files, each corresponding to a particular contig. I am familiar with samtools view (e.g. samtools view in.bam chr1 > chr1.bam) and have used it to make it for individual cases. Is there a better way to automate that?
Any help much appreciated.
jazz is offline   Reply With Quote
Old 06-05-2013, 12:16 PM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

No; I'd just write my own script that will read the line of the .bam, and then write it to the proper file for its contig.
swbarnes2 is offline   Reply With Quote
Old 06-05-2013, 01:13 PM   #3
vivek_
PhD Student
 
Location: Denmark

Join Date: Jul 2012
Posts: 164
Default

A unix oneliner should work right?

Code:
for i in {1..22};do samtools view -bh input.bam chr$i > chr$i.bam;done

Last edited by vivek_; 06-05-2013 at 01:16 PM.
vivek_ is offline   Reply With Quote
Old 06-05-2013, 01:29 PM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Quote:
Originally Posted by vivek_ View Post
A unix oneliner should work right?

Code:
for i in {1..22};do samtools view -bh input.bam chr$i > chr$i.bam;done
Depends, I've seen examples where the contigs weren't sequentially numbered (presumably due to some contigs becoming merged as latter data came in)

Also, for a file with a large number of contigs (have a look at some of the mouse lines from the Sanger Institute), looping over the whole file many many times will get super slow. You could probably write a script to process the whole thing in one go in a fraction of the time.
dpryan is offline   Reply With Quote
Old 06-05-2013, 01:37 PM   #5
vivek_
PhD Student
 
Location: Denmark

Join Date: Jul 2012
Posts: 164
Default

That's why you have the BAM index right, so you are not reading the entire file to export each coordinate?

for the sequentiality issue, you can extract the contig names from the bam header into a file and loop over them:

Code:
samtools view -H input.bam | awk '{print $2}' | awk '{gsub(/SN\:/,""); print}'  > contigs.txt
vivek_ is offline   Reply With Quote
Old 06-06-2013, 08:23 AM   #6
syfo
Just a member
 
Location: Southern EU

Join Date: Nov 2012
Posts: 103
Default

Quote:
Originally Posted by vivek_ View Post
That's why you have the BAM index right, so you are not reading the entire file to export each coordinate?

for the sequentiality issue, you can extract the contig names from the bam header into a file and loop over them:

Code:
samtools view -H input.bam | awk '{print $2}' | awk '{gsub(/SN\:/,""); print}'  > contigs.txt
Watch out, you don't want the first ("@HD ...") nor the last ("@PG ...") line of the header.

Try this instead:
Code:
samtools view -H all.bam | sed '1d;s/.*SN:\(.*\)\t.*/\1/;$d' > contigs.list
Or, if you prefer awk:
Code:
samtools view -H all.bam | awk '/^@SQ/{gsub(/SN\:/,"");print $2}' > contigs.list
or even (just for fun):
Code:
samtools idxstats all.bam | cut -f1 > contigs.list
All those should give you the same list of contigs.

Then,
Code:
for c in `cat contigs.list` ; do
echo processing $c
samtools view -bh all.bam $c > $c.bam
done
But I agree it might take a while...
syfo is offline   Reply With Quote
Old 06-06-2013, 10:59 AM   #7
jazz
Member
 
Location: us

Join Date: Nov 2008
Posts: 28
Default

Thanks everyone. I will give these suggestions a try and let you know how it went.
jazz is offline   Reply With Quote
Old 02-09-2015, 05:09 PM   #8
jjlaisnoopy
Junior Member
 
Location: Taiwan

Join Date: Sep 2011
Posts: 6
Default

This is good method to split bam file. But I got a question.
The following is idxstats of a bam file:
chrM 16571 2073252 32042
chr1 249250621 115733016 1937746
chr2 243199373 104133908 2244387
chr3 198022430 96577573 1501432
chr4 191154276 89582368 1825761
chr5 180915260 94818025 1486923
chr6 171115067 84533173 1273600
chr7 159138663 71186849 1531851
chr8 146364022 65630236 1315785
chr9 141213431 59368028 1184324
chr10 135534747 63503839 2018103
chr11 135006516 59963670 1373030
chr12 133851895 63898721 1180836
chr13 115169878 41939790 616704
chr14 107349540 43647215 758336
chr15 102531392 39227879 624791
chr16 90354753 42298502 991456
chr17 81195210 49043800 916042
chr18 78077248 75701725 1614444
chr19 59128983 26119207 755016
chr20 63025520 32668117 644231
chr21 48129895 19226969 547857
chr22 51304566 15797809 277926
chrX 155270560 74715396 1365662
chrY 59373566 2021162 642042
* 0 0 42640654

How could I get the "*" 42640654 reads, those were not mapped to any contigs ?
jjlaisnoopy is offline   Reply With Quote
Old 02-09-2015, 05:35 PM   #9
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,054
Default

Quote:
Originally Posted by jjlaisnoopy View Post
This is good method to split bam file. But I got a question.

* 0 0 42640654

How could I get the "*" 42640654 reads, those were not mapped to any contigs ?
https://www.biostars.org/p/56246/
GenoMax is offline   Reply With Quote
Old 02-10-2015, 04:50 PM   #10
jjlaisnoopy
Junior Member
 
Location: Taiwan

Join Date: Sep 2011
Posts: 6
Default

Quote:
Originally Posted by GenoMax View Post
I tried the parameter: -f 4
And then index the result bam file
Here is the idxstats from it:

chrM 16571 0 32042
chr1 249250621 0 1937746
chr2 243199373 0 2244387
chr3 198022430 0 1501432
chr4 191154276 0 1825761
chr5 180915260 0 1486923
chr6 171115067 0 1273600
chr7 159138663 0 1531851
chr8 146364022 0 1315785
chr9 141213431 0 1184324
chr10 135534747 0 2018103
chr11 135006516 0 1373030
chr12 133851895 0 1180836
chr13 115169878 0 616704
chr14 107349540 0 758336
chr15 102531392 0 624791
chr16 90354753 0 991456
chr17 81195210 0 916042
chr18 78077248 0 1614444
chr19 59128983 0 755016
chr20 63025520 0 644231
chr21 48129895 0 547857
chr22 51304566 0 277926
chrX 155270560 0 1365662
chrY 59373566 0 642042
* 0 0 42640654

Any suggestions ?
jjlaisnoopy is offline   Reply With Quote
Old 02-10-2015, 05:55 PM   #11
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,054
Default

Are you asking about why the number is not adding up to 42640654? See possible explanation here: https://www.biostars.org/p/18949/

Last edited by GenoMax; 02-10-2015 at 06:05 PM.
GenoMax is offline   Reply With Quote
Old 02-10-2015, 06:54 PM   #12
jjlaisnoopy
Junior Member
 
Location: Taiwan

Join Date: Sep 2011
Posts: 6
Default

All I want is the Unmapped reads with no mate or an unmapped mate are assigned to chrom "*" , not include unmapped mate reads which assigned to chr1, chr2, ...
just reads assigned to
"*" 0 0 42640654
jjlaisnoopy is offline   Reply With Quote
Old 02-10-2015, 11:31 PM   #13
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Just write a quick little script in python with pysam to do this. There isn't always a premade program to do everything.
dpryan is offline   Reply With Quote
Old 02-11-2015, 05:34 AM   #14
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,054
Default

A command line solution. See if this works:

Code:
$ samtools view -h file.bam | awk -F'\t' '{OFS = "\n"; ORS = "\n";}{ if ($3 == "*" ) print "@"$1,$10,"+",$11}' > outfile.fastq

Last edited by GenoMax; 02-11-2015 at 05:38 AM. Reason: correction
GenoMax is offline   Reply With Quote
Old 02-11-2015, 05:37 AM   #15
sarvidsson
Senior Member
 
Location: Berlin, Germany

Join Date: Jan 2015
Posts: 137
Default

Quote:
Originally Posted by GenoMax View Post
A command line solution. See if this works:

Code:
$ samtools view -h file.bam | awk -F'\t' '{OFS = "\n"; ORS = "\n";}{ if ($3 == "*" ) print "@"$1,$10,"+",$11}' > outfile.bam
GenoMax, you probably mean to call the output file "outfile.fastq", right?
Code:
$ samtools view -h file.bam | awk -F'\t' '{OFS = "\n"; ORS = "\n";}{ if ($3 == "*" ) print "@"$1,$10,"+",$11}' > outfile.fastq
sarvidsson is offline   Reply With Quote
Old 02-11-2015, 05:39 AM   #16
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,054
Default

Ah yes. Corrected.
GenoMax is offline   Reply With Quote
Old 02-11-2015, 09:25 PM   #17
jjlaisnoopy
Junior Member
 
Location: Taiwan

Join Date: Sep 2011
Posts: 6
Default

Thank you all !
jjlaisnoopy 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 05:54 PM.


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