Answers.Body: <p>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 &quot;unmapped&quot; reads. Otherwise you will not have the complete information. You can get the unmapped reads by putting <code>unmapped</code> into the <code>-L</code> argument.</p>
  <p>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).</p>
  <p>You could split the BAM into shards but only the &quot;proper pairs&quot; ( filter after GATK with <code>samtools view</code>, 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.</p>
  <p>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.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-12-30T15:48:48.763Z
Answers.LastActivityDate: 2023-12-30T16:00:26.540Z
Answers.LastEditDate: 2023-12-30T16:00:26.540Z
Answers.UserDisplayName: Yossi Farjoun
Body: <p>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.</p>
  <pre><code>samtools: samtools view -bo ${pair_id}.\\\$chr.bam $bam \\\$chr
  bamtools:bamtools split -in $bam -reference 
  </code></pre>
  <p>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.</p>
  <p>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?</p>
  <p>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?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-16T17:56:12.570Z
Comments.Text: Please [edit] your question and show us the exact commands you are using and the exact error message. My guess is that you're giving the index for the entire genome, so the error makes sense, but we can't know without the commands. Can you explain why you are splitting the bam though? Why not get the stats for the entire bam, and then split those by chromosome?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-17T07:29:39.857Z
Comments.Text: Can you please explain what metrics you are interested in? If it's something equivalent to "number of reads mapped per chromosome", there's probably a better way to get those statistics.
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-17T18:50:51.713Z
Comments.Text: Thanks for the edit, but we still don't know what Picard commands you are giving. What index files are you using? Is it an index for the entire genome and you're feeding only a per-chr fasta file? Are you referring to the bam index? Have you indexed the individual bam files? Also, 70G isn't huge, that's a pretty normal file size for this kind of thing, so you might want to consider more powerful hardware if that is an option.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-12-18T02:12:56.070Z
Comments.Text: Thank you! So first I aligned with the whole genome, and then tried to split the bam into different chromsome, then to do the picard using the per chrom ref.  Thank you for your comment on the size(70G not huge), I have decided to just to use large memory rather than split them. I am OK if this post is to removed from this site, since it might not be helpful for anyone
Comments.UserDisplayName: cautree
Tags: <phylogenetics><picard>
Title: align to the whole hg38 genome, then split bam and collect metrics on each bam issue
ThreadText: # align to the whole hg38 genome, then split bam and collect metrics on each bam issue
  
  (Asked by: username_0 on 2023-12-16)
  
  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?
  
  **Tags:** phylogenetics, picard
  
  > **Comments:**
  > **username_1**: Please [edit] your question and show us the exact commands you are using and the exact error message. My guess is that you're giving the index for the entire genome, so the error makes sense, but we can't know without the commands. Can you explain why you are splitting the bam though? Why not get the stats for the entire bam, and then split those by chromosome?  
  > **username_2**: Can you please explain what metrics you are interested in? If it's something equivalent to "number of reads mapped per chromosome", there's probably a better way to get those statistics.  
  > **username_1**: Thanks for the edit, but we still don't know what Picard commands you are giving. What index files are you using? Is it an index for the entire genome and you're feeding only a per-chr fasta file? Are you referring to the bam index? Have you indexed the individual bam files? Also, 70G isn't huge, that's a pretty normal file size for this kind of thing, so you might want to consider more powerful hardware if that is an option.  
  > **username_0**: Thank you! So first I aligned with the whole genome, and then tried to split the bam into different chromsome, then to do the picard using the per chrom ref.  Thank you for your comment on the size(70G not huge), I have decided to just to use large memory rather than split them. I am OK if this post is to removed from this site, since it might not be helpful for anyone  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_3 on 2023-12-30. Score: 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.
  
  
  ---
  
  