How To Get All Contig Boundaries In A Bam File.

Question

[dustar1986]

Hi everyone,

Is there a smart way to get all contig boundaries after assembly? I know samtools can get number of sequences covering a defined region. Then I could write a pipeline to find each contig's boundary by looking for a region flank by coverage of 0. It seems this costs a lot of time. Is there an alternative and quicker way to know the start and end position of each contig?

Thanks a lot.

Dadi

[Aaron Statham]

I'm not sure I understand your question - are you aligning reads to already assembled contigs or assembling contigs from reads?

[dustar1986]

I am talking about assembling contigs from reads. Sorry to confuse you.

Answer 1

[Pierre Lindenbaum]

You could use samtools mpileup and awk to get the contigs:

samtools  mpileup file.bam |\
cut -d '        ' -f 1,2 |\
awk -F '        ' 'BEGIN {chr="";start=-1;end=-1} {if(chr!=$1 || int($2)!=end+1) { if(chr!="") {printf("%s:%d-%d\n",chr,start,end);} chr=$1;start=int($2);end=int($2);} else { end=end+1;}} END {if(chr!="") {printf("%s:%d-%d\n",chr,start,end); } }'
chr1:10044-10300
chr1:10330-10468
chr1:10504-10557
chr1:10585-10633
chr1:11302-11355
chr1:11594-11647
chr1:11696-11823
chr1:11827-13107
chr1:13122-13840
chr1:13874-13973
(...)

[dustar1986]

Thank you very much. Extremely helpful. That's exact what I wanna get.

[Pierre Lindenbaum]

update 10 years later: just pipe into `bedtools merge|

Answer 2

[Pablo]

I don't know if I got your question right, but if you have a BAM file, doing this:

samtools view -H file.bam| grep "^@SQ"

you will get the length of each contig in the "LN" sub-field (the start is always 1).

Apologies in advanced if I misunderstood you...

[dustar1986]

I'm sorry maybe I didn't make it clear. If there are only 3 reads within the sequencing data. Sequence 1: chr1:100-150; Sequence 2: chr1:110-155; Sequence 3: chr1:130-175. Then it should be assembled into a contig:chr1:100-175. I am asking if I can get all the contigs position after reads assembly.

Answer 3

[Alex Reynolds]

You could use a sam2bed conversion step with a bedops --merge operation:

$ sam2bed < reads.sam | bedops --merge - > mergedReadCoordinates.bed

Note that this conversion from SAM to BED will change 1-indexed coordinates to 0-indexed coordinates.

Answer 4

[sureshbioinfo]

if you do in samtools idxstats with your file.bam

samtools idxstats file.bam > file.stats.txt

you will get output three columns with Contigname, contiglenth, no of reads

[zjamal.msbi19rcms]

how to obtain contig sequence?