SEQanswers

Go Back   SEQanswers > Sequencing Technologies/Companies > Pacific Biosciences



Similar Threads
Thread Thread Starter Forum Replies Last Post
blasr SAM output tags Aizon Pacific Biosciences 5 04-22-2015 12:27 AM
BLASR becoming very slow. What's the cause? myrs Pacific Biosciences 0 02-05-2015 09:27 AM
Get consensus sequence from blasr alignment guibar Pacific Biosciences 4 04-10-2014 09:55 AM
BLASR output - what does each column mean? ritzriya Pacific Biosciences 7 01-20-2014 10:52 PM

Reply
 
Thread Tools
Old 12-13-2016, 07:40 PM   #1
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default 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:

https://github.com/PacificBioscience...-is-deprecated

Ultimately, in short, I simply want to utilize the bax.h5 files with quiver (or arrow) to polish the assembly I did in Canu.
jpummil is offline   Reply With Quote
Old 12-14-2016, 10:38 AM   #2
rhall
Senior Member
 
Location: San Francisco

Join Date: Aug 2012
Posts: 318
Default

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
rhall is offline   Reply With Quote
Old 12-14-2016, 10:43 AM   #3
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

Quote:
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!
jpummil is offline   Reply With Quote
Old 12-14-2016, 10:47 AM   #4
rhall
Senior Member
 
Location: San Francisco

Join Date: Aug 2012
Posts: 318
Default

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.
rhall is offline   Reply With Quote
Old 12-14-2016, 10:49 AM   #5
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

Quote:
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.
jpummil is offline   Reply With Quote
Old 12-15-2016, 05:39 PM   #6
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

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?
jpummil is offline   Reply With Quote
Old 12-16-2016, 10:52 AM   #7
rhall
Senior Member
 
Location: San Francisco

Join Date: Aug 2012
Posts: 318
Default

Yes the bax2bam should be ran per movie (3 files). PBalign can then be ran using all the data.
rhall is offline   Reply With Quote
Old 12-16-2016, 10:57 AM   #8
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

Quote:
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.
jpummil is offline   Reply With Quote
Old 12-16-2016, 11:19 AM   #9
rhall
Senior Member
 
Location: San Francisco

Join Date: Aug 2012
Posts: 318
Default

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.
rhall is offline   Reply With Quote
Old 12-19-2016, 04:09 PM   #10
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

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...
jpummil is offline   Reply With Quote
Old 12-19-2016, 04:31 PM   #11
gconcepcion
Member
 
Location: Menlo Park

Join Date: Dec 2010
Posts: 68
Default

Quote:
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 at 04:34 PM.
gconcepcion is offline   Reply With Quote
Old 12-20-2016, 11:43 AM   #12
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

Quote:
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!
jpummil is offline   Reply With Quote
Old 12-21-2016, 09:15 AM   #13
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

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
File Type: txt report.txt (2.0 KB, 13 views)
jpummil is offline   Reply With Quote
Old 01-18-2017, 07:08 AM   #14
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

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!
jpummil is offline   Reply With Quote
Old 01-20-2017, 08:02 AM   #15
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

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...
jpummil is offline   Reply With Quote
Old 01-21-2017, 06:57 PM   #16
jpummil
Member
 
Location: Fayetteville, AR

Join Date: Apr 2014
Posts: 82
Default

Apologies, attached is Quast report comparing base Canu assembly, Base + Quiver polishing, and Quiver + Pilon....

Note, still using --glimmer for gene finding in Quast. Need to switch over to GeneMark-ES (soon).
Attached Files
File Type: txt quast-report.txt (3.2 KB, 18 views)
jpummil is offline   Reply With Quote
Old 11-20-2018, 06:01 AM   #17
ro2006
Junior Member
 
Location: NL

Join Date: Nov 2018
Posts: 1
Default

Quote:
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
Hi everybody,
I hope my question is not too naive, but I'm quite new to the De Novo assembly. I did a De Novo assembly starting from PacBio reads using Canu.
Now I should polish the assembly and perform variant calling using arrow. My question is: in order to run pbalign as mentioned above, what's your reference.fasta file from the Canu output? Would this be the Contigs.fasta file?
Thanks in advance for any help
ro2006 is offline   Reply With Quote
Old 11-29-2018, 02:27 PM   #18
gconcepcion
Member
 
Location: Menlo Park

Join Date: Dec 2010
Posts: 68
Default

Quote:
Originally Posted by ro2006 View Post
Hi everybody,
I hope my question is not too naive, but I'm quite new to the De Novo assembly. I did a De Novo assembly starting from PacBio reads using Canu.
Now I should polish the assembly and perform variant calling using arrow. My question is: in order to run pbalign as mentioned above, what's your reference.fasta file from the Canu output? Would this be the Contigs.fasta file?
Thanks in advance for any help
Yes, the contigs.fasta that is output from Canu is what you want to `polish` so you will select that as your `reference.fasta` to which you map your raw reads using pbalign. Then variantCaller can subsequently call consensus for each base in your reference.fasta based on the raw read coverage
gconcepcion 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 09:00 PM.


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