[gconcepcion]

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

[ro2006]

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

[jpummil]

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

[jpummil]

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]

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]

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.

[jpummil]

$ 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!

[gconcepcion]

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

[jpummil]

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

[rhall]

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.

[jpummil]

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.

[rhall]

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

[jpummil]

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]

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.

[rhall]

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.