align to the whole hg38 genome, then split bam and collect metrics on each bam issue

Question

I am doing some coverage analysis on deep sequencing human genome. I first align to the whole hg38 genome, then split bam to each chrom, and collect metrics on bam for each chom. I first split using samtools. and then tried bamtools.

samtools: samtools view -bo ${pair_id}.\\\$chr.bam $bam \\\$chr
bamtools:bamtools split -in $bam -reference 

When I collect metrics I am using the reference fasta from each chrom. The picard metrics complains the index dictionary size is different: Sequence dictionaries are not the same size (25, 1). The picard method I used is: MarkDuplicates, CollectAlignmentSummaryMetrics, CollectGcBiasMetrics, CollectInsertSizeMetrics, CollectWgsMetrics.

My data is really huge, 70G for a sample, that is why I would like to split when doing the picard metrics. Seems picard metrics collection use a lot of memory. So for the 70G human fastq(30X), how big the memory is needed for running picard wgs metrics?

How can I deal with this? Can I modify the split bam header? How to do it? Or not split it into chromosmes, rather, using an interval list when collecting the metrics? will that use less memory?

Answer 1

This can be done, but you have to split the BAMs while keeping the full header (with the dictionary) in place for all the shards. This can be done with gatk PrintReads. One thing to remember is that you all need to have one shard with all the "unmapped" reads. Otherwise you will not have the complete information. You can get the unmapped reads by putting unmapped into the -L argument.

That said, MarkDuplicates cannot be run on shards since it needs both reads from every pair and some of the read pairs might be split. So you'll need to mark-duplicates on the entire 70GB bam (sorry).

You could split the BAM into shards but only the "proper pairs" ( filter after GATK with samtools view, but only on the sam flag so you don't mess with the dictionary) which will make sure that your pairs are on the same contig. Then put all the not proper pairs into their own shard.

You can also do the proper/improper filtering with GATK, but it involves read filters which are a bit more complex than samtools flag filters.