Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Blasr/Quiver - Best Practice for Current Versions?

    So, first time attempting to use quiver to polish an assembly I created in Canu. All of the software is installed and functional via command line (quiver, pbalign, bax2bam, quiver, arrow, etc). Our systems people installed using pitchfork.

    I have the assembly in .fasta format and I also have 30ea bax.h5 files to utilize. The issue comes when I start looking at blasr which apparently no longer outputs cmp.h5 files that are preferred by pbalign. blasr now defaults to .bam files (nice for standardization I guess) and it isn't obvious that pbalign will utilize .bam as input. Found this little snippet which was somewhat informative:

    GitHub is where people build software. More than 100 million people use GitHub to discover, fork, and contribute to over 420 million projects.


    Ultimately, in short, I simply want to utilize the bax.h5 files with quiver (or arrow) to polish the assembly I did in Canu.

  • #2
    Basic workflow using the data you have and the latest versions of software installed using pitchfork:
    bax.h5 -> bax2bam -> reads.bam
    reads.bam + reference.fasta -> pbalign -> aligned_reads.bam
    aligned_reads.bam -> variant_caller (quiver or arrow) -> consensus.fasta

    Comment


    • #3
      Originally posted by rhall View Post
      Basic workflow using the data you have and the latest versions of software installed using pitchfork:
      bax.h5 -> bax2bam -> reads.bam
      reads.bam + reference.fasta -> pbalign -> aligned_reads.bam
      aligned_reads.bam -> variant_caller (quiver or arrow) -> consensus.fasta
      Thanks so much, rhall! Just getting my data set up for a trial now ;-)

      I'll update the thread when I'm done!

      Comment


      • #4
        Just let me know if you need any more details on a particular step. If it's a particularly large dataset http://pacificbiosciences.github.io/...coretools.html has info on how to do a lightweight chunking without using the full workflow engine.

        Comment


        • #5
          Originally posted by rhall View Post
          Just let me know if you need any more details on a particular step. If it's a particularly large dataset http://pacificbiosciences.github.io/...coretools.html has info on how to do a lightweight chunking without using the full workflow engine.
          Genome size is about 170Mb....and I have ~65GB of bax.h5 files. BUT...I also have a pretty robust linux server to run it on.

          Comment


          • #6
            Hmmm...apparently some misunderstanding about "movies" on my part.

            Using the following bax2bam command:
            bax2bam Q1133.bax.h5/*.bax.h5 -o movie --subread


            ERROR: multiple movies detected:
            ERROR: m160611_230902_00125_c101003722550000001823226809161614_s1_p0
            ERROR: m160612_032816_00125_c101003722550000001823226809161615_s1_p0
            ERROR: m160612_074931_00125_c101003722550000001823226809161616_s1_p0
            ERROR: m160614_131714_00125_c101010592550000001823230710211662_s1_p0
            ERROR: m160614_194151_00125_c101010592550000001823230710211663_s1_p0
            ERROR: m160615_020159_00125_c101010592550000001823230710211664_s1_p0
            ERROR: m160615_062522_00125_c101010592550000001823230710211665_s1_p0
            ERROR: m160615_104711_00125_c101010592550000001823230710211666_s1_p0
            ERROR: m160615_150655_00125_c101010592550000001823230710211667_s1_p0
            ERROR: m160615_213630_00125_c101009882550000001823230710211620_s1_p0

            Contents of my bax.h5 directory:

            [jpummil@tres-l1 Q1133.bax.h5]$ ls
            m160611_230902_00125_c101003722550000001823226809161614_s1_p0.1.bax.h5
            m160611_230902_00125_c101003722550000001823226809161614_s1_p0.2.bax.h5
            m160611_230902_00125_c101003722550000001823226809161614_s1_p0.3.bax.h5
            m160612_032816_00125_c101003722550000001823226809161615_s1_p0.1.bax.h5
            m160612_032816_00125_c101003722550000001823226809161615_s1_p0.2.bax.h5
            m160612_032816_00125_c101003722550000001823226809161615_s1_p0.3.bax.h5
            m160612_074931_00125_c101003722550000001823226809161616_s1_p0.1.bax.h5
            m160612_074931_00125_c101003722550000001823226809161616_s1_p0.2.bax.h5
            m160612_074931_00125_c101003722550000001823226809161616_s1_p0.3.bax.h5
            m160614_131714_00125_c101010592550000001823230710211662_s1_p0.1.bax.h5
            m160614_131714_00125_c101010592550000001823230710211662_s1_p0.2.bax.h5
            m160614_131714_00125_c101010592550000001823230710211662_s1_p0.3.bax.h5
            m160614_194151_00125_c101010592550000001823230710211663_s1_p0.1.bax.h5
            m160614_194151_00125_c101010592550000001823230710211663_s1_p0.2.bax.h5
            m160614_194151_00125_c101010592550000001823230710211663_s1_p0.3.bax.h5
            m160615_020159_00125_c101010592550000001823230710211664_s1_p0.1.bax.h5
            m160615_020159_00125_c101010592550000001823230710211664_s1_p0.2.bax.h5
            m160615_020159_00125_c101010592550000001823230710211664_s1_p0.3.bax.h5
            m160615_062522_00125_c101010592550000001823230710211665_s1_p0.1.bax.h5
            m160615_062522_00125_c101010592550000001823230710211665_s1_p0.2.bax.h5
            m160615_062522_00125_c101010592550000001823230710211665_s1_p0.3.bax.h5
            m160615_104711_00125_c101010592550000001823230710211666_s1_p0.1.bax.h5
            m160615_104711_00125_c101010592550000001823230710211666_s1_p0.2.bax.h5
            m160615_104711_00125_c101010592550000001823230710211666_s1_p0.3.bax.h5
            m160615_150655_00125_c101010592550000001823230710211667_s1_p0.1.bax.h5
            m160615_150655_00125_c101010592550000001823230710211667_s1_p0.2.bax.h5
            m160615_150655_00125_c101010592550000001823230710211667_s1_p0.3.bax.h5
            m160615_213630_00125_c101009882550000001823230710211620_s1_p0.1.bax.h5
            m160615_213630_00125_c101009882550000001823230710211620_s1_p0.2.bax.h5
            m160615_213630_00125_c101009882550000001823230710211620_s1_p0.3.bax.h5

            So, does one just process the three files per group into a single .bam file, do the next, etc...then repeat the PBalign step per each .bam created?

            Comment


            • #7
              Yes the bax2bam should be ran per movie (3 files). PBalign can then be ran using all the data.

              Comment


              • #8
                Originally posted by rhall View Post
                Yes the bax2bam should be ran per movie (3 files). PBalign can then be ran using all the data.
                Thanks rhall! Was just about to post an update when you posted this ;-)

                -rw-r--r-- 1 jpummil jpummil 2.3G Dec 16 12:06 m160611_230902_00125_c101003722550000001823226809161614_s1_p0.1.bax.h5

                -rw-r--r-- 1 jpummil jpummil 2.2G Dec 16 12:06 m160611_230902_00125_c101003722550000001823226809161614_s1_p0.2.bax.h5

                -rw-r--r-- 1 jpummil jpummil 2.2G Dec 16 12:06 m160611_230902_00125_c101003722550000001823226809161614_s1_p0.3.bax.h5

                BECOMES (using bax2bam):

                -rw-rw-r-- 1 jpummil jpummil 1.3G Dec 16 12:18 movie.scraps.bam

                -rw-rw-r-- 1 jpummil jpummil 2.3G Dec 16 12:18 movie.subreads.bam


                Now running PBalign (blasr aligner) of the one lane of converted data to create "Q1133.aligned.bam".

                I'll then see about a Quiver run of that file just to validate the process. Afterwards, I'll go back and tidy things up to process all ten lanes of data, put thru PBalign, then Quiver.

                Comment


                • #9
                  FYI, the scraps file includes all the sequence not in the HQ (high quality) region, and the adapter sequences. The old bax format saved everything with indexes into the usable data. The subreads.bam contains all the usable data.

                  Comment


                  • #10
                    Used bax2bam to convert all of the bax.h5 files to .bam format...no problem.

                    When I try to use pbalign on the 10ea .bam files, it complains after it appears to read and accept the first three files:

                    pbalign: error: unrecognized arguments: movie_04.subreads.bam movie_05.subreads.bam movie_06.subreads.bam movie_07.subre
                    ads.bam movie_08.subreads.bam movie_09.subreads.bam movie_10.subreads.bam Q1133.092116.GS170M.contigs.fasta Q1133.aligne
                    d.bam

                    PBalign comand as follows:
                    pbalign --nproc 32 *.subreads.bam Q1133.092116.GS170M.contigs.fasta Q1133.aligned.bam

                    *.subreads.bam corresponds to the following filenames:
                    movie_01.subreads.bam
                    thru
                    movie_10.subreads.bam

                    I don't see anything obvious in the pbalign --help details about running multiple .bam files for alignment.

                    My initial tests using just a single lane worked as expected. It's when I get to the task of running 10ea lanes thru pbalign that things go south.

                    I'll keep looking for docs, etc...but any tips are appreciated. It just seems as if the PacBio help documents are a bit lacking on details...

                    Comment


                    • #11
                      Originally posted by jpummil View Post
                      Used bax2bam to convert all of the bax.h5 files to .bam format...no problem.

                      When I try to use pbalign on the 10ea .bam files, it complains after it appears to read and accept the first three files:

                      pbalign: error: unrecognized arguments: movie_04.subreads.bam movie_05.subreads.bam movie_06.subreads.bam movie_07.subre
                      ads.bam movie_08.subreads.bam movie_09.subreads.bam movie_10.subreads.bam Q1133.092116.GS170M.contigs.fasta Q1133.aligne
                      d.bam

                      PBalign comand as follows:
                      pbalign --nproc 32 *.subreads.bam Q1133.092116.GS170M.contigs.fasta Q1133.aligned.bam

                      *.subreads.bam corresponds to the following filenames:
                      movie_01.subreads.bam
                      thru
                      movie_10.subreads.bam

                      I don't see anything obvious in the pbalign --help details about running multiple .bam files for alignment.

                      My initial tests using just a single lane worked as expected. It's when I get to the task of running 10ea lanes thru pbalign that things go south.

                      I'll keep looking for docs, etc...but any tips are appreciated. It just seems as if the PacBio help documents are a bit lacking on details...

                      $ ls *.subreads.bam >input.fofn
                      $ pbalign --nproc 32 input.fofn Q1133.092116.GS170M.contigs.fasta Q1133.aligned.bam


                      pbalign takes individual sequence files as input, or a File of File Names (*.fofn)

                      From pbalign --help:
                      "Input can be a fasta, pls.h5, bas.h5 or ccs.h5 file or a fofn (file of file names). Output can be in SAM or BAM format."
                      Last edited by gconcepcion; 12-19-2016, 04:34 PM.

                      Comment


                      • #12
                        Originally posted by gconcepcion View Post
                        $ ls *.subreads.bam >input.fofn
                        $ pbalign --nproc 32 input.fofn Q1133.092116.GS170M.contigs.fasta Q1133.aligned.bam


                        pbalign takes individual sequence files as input, or a File of File Names (*.fofn)

                        From pbalign --help:
                        "Input can be a fasta, pls.h5, bas.h5 or ccs.h5 file or a fofn (file of file names). Output can be in SAM or BAM format."
                        Thanks gconcepcion! That definitely did the trick. Ran into a /tmp space issue next. THINK I have it solved by supplying the --tmpDir option in the command line and setting it to the /local_scratch on the compute node as well...almost a terabyte available there. Will know very soon ;-)

                        Onwards!

                        Comment


                        • #13
                          SUCCESS!

                          Not a huge improvement on the original assembly (using Canu)...but improved gene's found using Quast (comparison attached).

                          I used default parameters in Quiver, not sure tweaking would've changed a whole lot.

                          I'm "officially" on vacation 'til Jan 9th, but will try and toss together a short "How-To" based on the steps I followed (and those tips graciously supplied by rhall and gconcepcion!) so that others just starting to work with PacBio data and tools can do this as well.
                          Attached Files

                          Comment


                          • #14
                            So, I still need to write up a brief process workflow to describe what I've doe thus far, but thought I'd update anyone following this thread on another step I heard about thru the grapevine (people at PAG this week sending me email).

                            APPARENTLY....many are performing an additional step after the Quiver polishing and running the updated assembly thru BROAD's Pilon tool. Of course, you have to create a new alignment file with the original bax.h5 files (converted to .bams with bax2bam utility) and the Quiver polished assembly....

                            It's running now, so I'll report back (bad or good) when it completes and I can run a quick comparison of the three with Quast.

                            It should be said that I have been just using default settings with most of the tools, so probably some incremental gains to be had by tweaking settings IF you know your data well. That said, it seems developers do a pretty nice job of setting defaults these days!

                            Comment


                            • #15
                              Results of using Pilon AFTER having used Quiver....

                              So, Pilon just completed after 40+ hours of computation....another small increase in file size from what Quiver provided...

                              -rw-rw-r-- 1 jpummil jpummil 121958163 Jan 20 08:50 Q1133.012017.GS170M.Canu-Base.fasta

                              -rw-rw-r-- 1 jpummil jpummil 122467877 Jan 20 02:16 Q1133.122116.GS170M.Quiver.fasta (+ 509,714 bytes)

                              -rw-rw-r-- 1 jpummil jpummil 122649681 Jan 20 02:16 Q1133.012017.GS170M.Pilon.fasta (+ 181,804 bytes)

                              Still need to do a Quast assessment, maybe sometime today...

                              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
                              10 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, Yesterday, 06:07 PM
                              0 responses
                              10 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 03-22-2024, 10:03 AM
                              0 responses
                              51 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