A new tool to transfer RX tag from an unaligned bam. #7739

takutosato · 2022-03-28T20:05:02Z

Some read tags get lost when we convert SAM to fastq. This tool allows us to get those tags back once we are done processing the fastqs (some tools e.g. adapter clippers cannot take SAMs as input so the conversion is unavoidable.) So this tool works like Picard MergeBamAlignment, except that we are putting the tags from the unaligned bam to the aligned bam, rather than adding alignment info to the unaligned bam. We will use this in our new TCap RNA pipeline.

takutosato · 2022-03-28T20:10:29Z

@droazen, @lbergelson I have the following argument to the tool:

@Argument(fullName = "read-tags", doc = "read tag names to recover")
public List<String> readTags = DEFAULT_READ_TAGS;

On the command line I want to say

java -jar gatk.jar …
...
--read-tags RX

and want readTags to be a singleton list. But in my test it's not parsing the arguments correctly---what am I doing wrong here?

cmnbroad · 2022-03-28T20:35:23Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/RecoverReadTags.java

+    SAMFileGATKReadWriter writer;
+
+    final static String UMI_TAG = "RX";
+    final static List<String> DEFAULT_READ_TAGS = Arrays.asList(UMI_TAG);


Suggested change

final static List<String> DEFAULT_READ_TAGS = Arrays.asList(UMI_TAG);

final static List<String> DEFAULT_READ_TAGS = new ArrayList() {{

add(UMI_TAG);

}};

@takutosato The list returned by Arrays.asList() is immutable, so the parser can't modify the DEFAULT_READ_TAGS list. See the suggested alternative above.

I was getting errors (warnings, really) related to serializable UID stuff, which I don't understand well. The simpler approach seems to be to just have an empty list be the default, so that's what I did in the latest commit.

cmnbroad · 2022-03-28T20:39:16Z

...test/java/org/broadinstitute/hellbender/tools/walkers/qc/RecoverReadTagsIntegrationTest.java

+
+    @Test
+    public void test(){
+        final String home = "/Volumes/dsde_working/tsato/hydro.gen/Analysis/874_twist_RNA/add_umi/";


While I'm here, I noticed these local file refs...

gatk-bot · 2022-03-28T20:56:05Z

Travis reported job failures from build 38286
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
integration	openjdk11	38286.12	logs
integration	openjdk8	38286.2	logs

takutosato · 2022-03-28T21:17:46Z

@cmnbroad thanks for your review. I changed the default behavior to work around the issue and now the test is passing.

takutosato · 2022-03-31T18:40:41Z

@cmnbroad is my workaround satisfactory or should I try to figure out the serial version uid issue?

jamesemery

A few changes but for a new tool that is marked as beta/experimental it should be alright to merge. there are some tests and edge cases you might want to consider with your traversal though. You seem to have put a bunch of thought into iterating over sample1 assuming that sample2 might have a bunch of extra reads in it (not vice versa though) and my question is why?

jamesemery · 2022-04-05T14:38:27Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+@CommandLineProgramProperties(
+        summary = "Incorporate read tags in a SAM file to that of a matching SAM file",
+        oneLineSummary = "Incorporate read tags in a SAM file to that of a matching SAM file",
+        programGroup = QCProgramGroup.class // sato: change


jamesemery · 2022-04-05T14:38:42Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+        oneLineSummary = "Incorporate read tags in a SAM file to that of a matching SAM file",
+        programGroup = QCProgramGroup.class // sato: change
+)
+public class TransferReadTags extends GATKTool {


given that this is a new tool you might consider marking this as beta or experimental.

documentedfeature

Added Experimental and DocumentedFeature.

jamesemery · 2022-04-05T14:39:25Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+import java.util.List;
+
+/**
+ * This tool is designed for a pair of SAM files sharing the same read names (e.g.


designed for what on that pair?

elaborated.

jamesemery · 2022-04-05T14:40:25Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+ * due to the fact that some tools, such as the adapter clipping tools, are written for fastq only,
+ * while others are written specifically for SAM files.
+ *
+ * This tools behaves similarly to Picard MergeBamAlignment (MBA). The difference is that whereas


this can be a doc link with {@link ....MergeBamAlignment}

jamesemery · 2022-04-05T14:44:54Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+    public List<String> readTags = new ArrayList<>();
+
+
+    PeekableIterator<GATKRead> alignedSamIterator;


all of these should be private

jamesemery · 2022-04-05T14:48:41Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+    public File outSamFile;
+    SAMFileGATKReadWriter writer;
+
+    @Argument(fullName = "read-tags", doc = "read tag names to transfer")


add a note here that this will explode if any read lacks the requested tag (or change that behavior)

I think exiting (with a user exception) when the tag is not found is OK. Updated the doc.

jamesemery · 2022-04-05T14:53:55Z

...est/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTagsIntegrationTest.java

+
+    @Test
+    public void test() {
+        final File alignedSAMFile = createTempFile("sam1", ".bam");


Optional but I would consider pulling out all of this sam construction code DataProvider that just returns paths to the two bams and a list of tags that were artificially covered. This would make it easy to test other edge cases

Speaking of edge casees. I would recommend adding a test that asserts what the behavior is if both samples have the tag present. Should it overwrite the results? I believe that is what currently happens.

Done dataProvider

Done test with the case where the tag is already present.

jamesemery · 2022-04-05T15:00:10Z

src/main/java/org/broadinstitute/hellbender/utils/help/HelpConstants.java

@@ -56,6 +56,9 @@ public static String forumPost(String post) {
    public final static String DOC_CAT_METHYLATION_DISCOVERY = "Methylation-Specific Tools";
    public final static String DOC_CAT_METHYLATION_DISCOVERY_SUMMARY = "Tools that perform methylation calling, processing bisulfite sequenced, methylation-aware aligned BAM";

+    public final static String DOC_CAT_QC = "QC Tools";
+    public final static String DOC_CAT_QC_SUMMARY = "Tools that perform QC-related or other data processing tasks";


I actually think this tool fits under the banner of ReadDataManipulationProgramGroup or under Diagnostics and Quality Control which already exist as groups rather than making a new bespoke command line grouping for this one tool.

Deleted, using ~~Diagnostics and QC~~ ReadDataManipulation.

jamesemery · 2022-04-05T15:00:50Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+        programGroup = QCProgramGroup.class // sato: change
+)
+public class TransferReadTags extends GATKTool {
+    final static String UMI_TAG_NAME = "RX";


jamesemery · 2022-04-05T15:21:51Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+        // will need to use this later
+        while (alignedSamIterator.hasNext()){
+            currentRead1 = alignedSamIterator.next();
+            int diff = queryNameComparator.compareReadNames(currentRead1, currentRead2);


I'm a little worried about what happens when this QueryNameComparitor runs up against samtools sorted reads (which sorts the reads not lexicographicallly but lexicographically with the tile/x/y coordinates numerically sorted which we don't do by default) it seems like your traversal of pull on the target bam, then pull on the unmapped bam until you reach might run into trouble if your bam uses the samtools sorting for example on these reads: XXXXXXX:2:130:1000 followed by XXXXXXX:2:1200:1000 (which is still legal by the other style of queryname sorting). After processing the first read (at 2:130) and then loading currentRead1 with XXXXXXX:2:1200:1000 and currentRead2 will still be loaded with the last position XXXXXXX:2:130:1000. When you compute the diff with QueryNameComparitor itis going to consider XXXXXXX:2:130:1000 > XXXXXXX:2:1200:1000 and which i think will fall into the first IllegalStateException in this method.

Let me ask you. Is there a reason you cant just pull on both iterators in tandem and blow up if there is ever a mismatch between the reads and samflags? Do you expect the unmapped bam to have more reads in it for some reason? I might consider reading AbstractAlignmentMerger.mergeAlignment() in picard to see how they handle iterating the two lists in tandem.

Excellent point. STAR aligner drops some reads when it creates the transcriptome-aligned bam, so I cannot assume that there are an equal number of reads in both aligned and unaligned. I need to think more carefully about the case you brought up though.

OK phew it look like samtools query-name sorts the way you described, and Picard SortSam (which is what we are using in the RNA pipeline) does the string compare (see SAMRecordQueryNameComparator). So the code as is should work; how about we add the samtools version of comparator after this week's release (since we don't need it for this pipeline if not ever)

This might or might not be an issue? we don't support writing bams labeled as queryname ordered in the samtools sense but we do support reading them? That will mean you get cryptic errors here if you aren't careful. Its also somewhat of an edge cases that likely won't show up in your pipelines

Won't be an issue in my pipeline. But I updated the error message to make it less cryptic.

jamesemery · 2022-04-05T15:40:47Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+    PeekableIterator<GATKRead> alignedSamIterator;
+    PeekableIterator<GATKRead> unmappedSamIterator;
+
+    GATKRead currentRead1;


rename these to something more informative like "targetRead" and/or "annotationRead"

Good point. Done.

lbergelson · 2022-04-05T15:49:31Z

src/main/java/org/broadinstitute/hellbender/cmdline/programgroups/QCProgramGroup.java

+import org.broadinstitute.barclay.argparser.CommandLineProgramGroup;
+import org.broadinstitute.hellbender.utils.help.HelpConstants;
+
+public class QCProgramGroup implements CommandLineProgramGroup {


QC seems like a weird name for these. I would either put it either in ReadDataManipulation or DiagnosticsandQC. Alternatively, this seems like maybe a very specific tool for a certain pipeline, so maybe something like RNAReadData or something like that.

Good point, I think ReadDataManipulation makes sense the most.

takutosato · 2022-04-06T15:14:53Z

@jamesemery thanks for reviewing. back to you

gatk-bot · 2022-04-06T16:38:54Z

Travis reported job failures from build 38522
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
integration	openjdk11	38522.12	logs

jamesemery

Two very minor documentation comments and then I would say this is good to merge.

jamesemery · 2022-04-07T18:45:01Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+ * base and adds the read tags from the unaligned bam to the aligned bam.
+ *
+ * Currently, the tool is implemented for the specific case of transfering UMI read tags (RX)
+ * from an unaligned bam.


Add an example usage to this.

jamesemery · 2022-04-07T18:47:00Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/qc/TransferReadTags.java

+                        writer.addRead(updatedRead);
+                        break;
+                    } else {
+                        throw new IllegalStateException("Aligned read is lexicographically smaller than the unmapped read. This tool assumes " +


I would say the regular case here is that a read is missing from the unmapped bam file and I wold probably make that the primary takeaway (with a note that thsi could also be caused by sorting woes).

takutosato added 4 commits March 28, 2022 11:34

first commit, not add test

f84ceb4

test draft 1

bbdc101

test written

dae8cd6

take a list of tags

5e56f2c

cmnbroad reviewed Mar 28, 2022

View reviewed changes

takutosato added 2 commits March 28, 2022 17:03

chris edit

39f744a

more chris edit

a9077df

takutosato added 2 commits March 29, 2022 14:07

name change

be7a086

continue to rename

964ee51

jamesemery requested changes Apr 5, 2022

View reviewed changes

jamesemery reviewed Apr 5, 2022

View reviewed changes

lbergelson reviewed Apr 5, 2022

View reviewed changes

james comment

b24c4ae

slightly more helpful error message

ed3ef20

jamesemery approved these changes Apr 7, 2022

View reviewed changes

doc modification

94ce12b

takutosato merged commit 2b0a558 into master Apr 7, 2022

takutosato deleted the ts_add_tag branch April 7, 2022 21:53

-    final static List<String> DEFAULT_READ_TAGS = Arrays.asList(UMI_TAG);
+        final static List<String> DEFAULT_READ_TAGS = new ArrayList() {{
+        add(UMI_TAG);
+    }};

		public List<String> readTags = new ArrayList<>();


		PeekableIterator<GATKRead> alignedSamIterator;

A new tool to transfer RX tag from an unaligned bam. #7739

A new tool to transfer RX tag from an unaligned bam. #7739

Conversation

takutosato commented Mar 28, 2022

takutosato commented Mar 28, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gatk-bot commented Mar 28, 2022 • edited Loading

takutosato commented Mar 28, 2022

takutosato commented Mar 31, 2022

jamesemery left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

takutosato Apr 5, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

takutosato Apr 5, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

takutosato commented Apr 6, 2022

gatk-bot commented Apr 6, 2022

jamesemery left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gatk-bot commented Mar 28, 2022 •

edited

Loading

takutosato Apr 5, 2022 •

edited

Loading

takutosato Apr 5, 2022 •

edited

Loading