Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • How to get all contig boundaries from a sorted bam file

    Hi everyone,

    Is there a smart way to assemble reads in the deep sequencing data?
    After aligning reads to the genome using bowtie and samtools, I got the bam file. I want to extend the reads to get contigs if these reads have the same ref chromosome and are on the DNA strand.
    Can I get the start position, end position, ref strand and chromosome number of all these contigs?

    I wrote the following code but it ignore the strand (+/-) the reads come from. Reads from - strand might connected to a + strand one. That's not good.
    Is that possible to extract ref strand info from samtools mpileup?

    My code:

    samtools mpileup input.sorted.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); } }'>output.txt

    And I really want to achieve this in the shell. Any suggestion?

    Thanks a lot.

    Dadi

  • #2
    Ignore me: Thought you were working from the SAM/BAM file directly:
    For the strand you'll have to parse the FLAG field, and restricting yourself to shell script makes this much more difficult than it need be. Is using a scripting language really not an option for you?

    Comment


    • #3
      Originally posted by maubp View Post
      Ignore me: Thought you were working from the SAM/BAM file directly:
      For the strand you'll have to parse the FLAG field, and restricting yourself to shell script makes this much more difficult than it need be. Is using a scripting language really not an option for you?
      No No. Using a scripting language, especially Python, is fine for me. Any suggestion on this please?
      I just thought samtools might be better to deal with sequencing data coz is more "professional" than a script...

      Comment


      • #4
        You can use the samtools C API from a script, e.g. for Python there is PySam,

        Comment

        Latest Articles

        Collapse

        • seqadmin
          Strategies for Sequencing Challenging Samples
          by seqadmin


          Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
          03-22-2024, 06:39 AM
        • seqadmin
          Techniques and Challenges in Conservation Genomics
          by seqadmin



          The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

          Avian Conservation
          Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
          03-08-2024, 10:41 AM

        ad_right_rmr

        Collapse

        News

        Collapse

        Topics Statistics Last Post
        Started by seqadmin, Yesterday, 06:37 PM
        0 responses
        8 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, Yesterday, 06:07 PM
        0 responses
        8 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 03-22-2024, 10:03 AM
        0 responses
        49 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 03-21-2024, 07:32 AM
        0 responses
        67 views
        0 likes
        Last Post seqadmin  
        Working...
        X