This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How To Get All Contig Boundaries In A Bam File.

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

next-gen sequencing contigs
• 8,984 views • link
updated by Alex Reynolds Alex Reynolds • written by dustar1986 dustar1986

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

• 0 views • link
written by Aaron Statham Aaron Statham

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

• 0 views • link
written by dustar1986 dustar1986

4 answers

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
(...)
• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

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

• 0 views • link
written by dustar1986 dustar1986

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

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

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...

• 1 views • link
updated by Ram Ram • written by Pablo Pablo

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.

• 0 views • link
written by dustar1986 dustar1986

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.

• 296 views • link
updated by Ram Ram • written by Alex Reynolds Alex Reynolds

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

• 0 views • link
updated by Ram Ram • written by sureshbioinfo • sureshbioinfo

Log in to answer this question.

More posts like this