Thanks Matt but we did actually request unaligned as well as the standard set of Casava processed stuff they issue, as we'll be running this through our own pipeline. I'm fairly certain they are supposed to be unaligned bams as Illumina's documentation provides details on how to extract the /1 and /2 fastq files from it. I wonder if it's something as simple as I've accidentally started on the wrong bam file.
Extracting Exome And Other Data From Whole Genome Unaligned Bam
Hi there,
Firstly, apologies for the overly long description, I'm rather tired and trying not to ramble too much!
To give you some background, I'm currently working with some human whole genome sequences which are in the form of paired-end, unaligned bam files that we received from Illumina. Ultimately we will be doing whole genome analysis, but initially we want to extract part of the sequences, namely we want to look at the exome and at a predefined set of loci of interest.
In order to achieve this I have created appropriate bed files with the intention of using intersectBed / pairedToBed from the Bed Tools package to extract the appropriate parts of the genome. However, after thoroughly confusing myself I am a little unsure as to the best way to proceed.
At first I had assumed that the reads would need to be aligned prior to hacking out the areas of interest using intersectBed. So I split the unaligned bam by chromosome and aligned it to hg37 using bwa (default values), it was split into the chromosomes to allow the sampe step to be run in parallel to reduce time.
After completing the alignment, converting to bam and sorting by queryname using Picard (apparently required in order to use pairToBed) I frustratingly ended up with errors when attempting to use pairToBed. From searching around I think the error was caused by singletons being aligned, something that pairToBed cannot currently handle. So instead I tried using intersectBed but I ended up getting an empty bam file as output, but no errors.
As an experiment I tried running intersectBed on the unaligned chromosome 1 bam and it gave a bam file that was appropriately sized for exome output!
What probably very obvious thing am I missing?
• 392 views
•
link
updated
by
Matt Shirley
•
written
by
Jamie •
•
written
by
Jamie •
1 answer
I highly doubt Illumina sent you an unaligned bam file, since your bam file has chromosome annotations (which come from mapping to a reference genome). Probably, your bam file was aligned to hg19 using Eland, and you might also have a file with SNV and indels. Unless you want to look at the aligned reads, you might consider starting with Illumina's variant calls, since you would need to call variants on the region of reads you extract. But, if you really want to extract region corresponding to your bed file from your aligned reads, you should consider using samtools:
samtools view -L file.bed file.bam > set.bam
• 0 views
•
link
written
by
Matt Shirley
Log in to answer this question.
More posts like this
-
Pseudo bulking of scATAC-seq
written by Manhezz •
-
DNA sequencing workflow involving UMIs
written by ww22runner •
-
how to split a Bam by chromosome and unplaced scaffolds.
written by Maxine •
-
PacBio merge runs
written by firestar
-
Union of unaligned fastq reads
written by Jeffin Rockey
-
Counting reads on genome
written by EVR
-
Split bam by multiple chromosome to a single bam file
written by l0o0
-
Genome data on hadoop chromosome level
written by shalini.ravishankar •
-
Remove Intronic Regions in .BAM
written by K.Nijbroek •
-
Advice On Adding Readgroups
written by Jamie •
how do you split your BAM per chromosome if the reads are unaligned? Are you sure your BAM contains only unaligned reads?
That could be it, I need to get more sleep... I'll check the data and confirm tomorrow, but that's a completely obvious thing I may have been missing.
I used samtools to split and also tried bam utils's split function (on separate occassions), which I guess should not have worked if the bam was unaligned? I'm still fairly sure these should be unaligned bams though... Curious, I'll look into further tomorrow.