Request for a protocol on annotation of DNA mismatches with REDItools, and help when DNA and RNA annotation don't match

[jonng1000]

Thanks so much for all your hard work replying everyone's questions here, I'm sure we all appreciate it =). I've successfully ran REDItools 3 based on the info provided in this repo python3 -m reditools analyze "$BAM" -o "$OUTPUT_DIR/${SAMPLE}_reditools.txt" -q 255 -me 5 -s 2 -C -t 24.

For some of my samples, I do have WGS/WES data. Hence I want to annotate REDItools3 output from the RNA bam files, with the DNA mismatches, and use this to help me filter away edit sites which could be DNA mutations. May I please know exactly how do I do it? Do I follow the standard DNA sequence analysis protocol which is the below?

1. Read QC and trimming
2. bwa-mem2 alignment
3. Read duplicate marking
4. BQSR
5. Alignment QC

If so, may I please confirm that typically used parameters for the above steps would be sufficient?

If my understanding is correct, are the next steps to use analyze and annotate in REDItools? If so, may I please know what exact commands to use?

Thank you so much. Please pardon my ignorance. I think I understand how to use REDItools for RNA-seq only, but am now learning how to use it for DNA bam files as well

[jonng1000]

Also, I've another issue which I would like your help on. I have ran reditools annotate after analyze. But even overlapping DNA and RNA positions are not annotated. May I get some help help here please?

((bioinfo) ) jonathan@comp1:~/work/rna_edit/pancreatic_lines/elly_dataset/jon_workflow/
dna_sites_wes$ grep 15820 annotated_output/SRR8615241_SRR8619007_annotated.txt | head
chr1    15820   G       *       106     37.48   [0, 0, 64, 42]  GT      0.40    -       - -    --

This is my annotated RNA (SRR8615241) and DNA (SRR8619007) file. This shows that chr1 15820 is my RNA site

((bioinfo) ) jonathan@comp1:~/work/rna_edit/pancreatic_lines/elly_dataset/jon_workflow/
dna_sites_wes$ grep 15820 reditools_output/SRR8619007_reditools.txt | head
chr1    15820   G       *       4       37.75   [0, 0, 2, 2]    GT      0.50    -       -       - -- 

And this shows that my DNA file has the same site. But in the above combined annotation file, the DNA info is not present.

From the RNA and DNA reditools output file, I see that the order of the chromosomes are not the same. Does this affect it? If so, lets say my RNA has the chr1-22, X,Y and M chromosomes. But it also has other contigs like chr1_ABCDEF. In my DNA file, can I remove all of these contigs, and just order the chr1-22, X,Y and M chromosomes in the same order as the RNA file? Will this work? If my RNA file has contigs which doesn't appear in the DNA file, will this affect things?

My RNA and DNA bam files were aligned against the same reference files. But the RNA bam files have the contigs renamed to the UCSC names, but the DNA bam files werent, as the UCSC renaming tool has an error here. But for both of them, the main chromosomes like chr1, chr2 etc have the same names.

[ahanden]

For some of my samples, I do have WGS/WES data. Hence I want to annotate REDItools3 output from the RNA bam files, with the DNA mismatches, and use this to help me filter away edit sites which could be DNA mutations. May I please know exactly how do I do it? Do I follow the standard DNA sequence analysis protocol which is the below?

1. Read QC and trimming
2. bwa-mem2 alignment
3. Read duplicate marking
4. BQSR
5. Alignment QC

This setup looks fine. In general, you should use the same preprocessing steps as you would for variant calling.

Marking duplicates is particularly important for DNA sequencing data as REDItools will discard any read marked as a duplicate during processing.

Personally, I do not recommend BQSR. I have read too many reports of reduced accuracy of variant calling from the tool, and most modern systems seem to have overcome the PHRED score issue that BQSR was originally designed to address. You may want to try running REDItools with and without BQSR adjusted BAM files to see if there is a difference.

If my understanding is correct, are the next steps to use analyze and annotate in REDItools? If so, may I please know what exact commands to use?

For your DNA sample, based on your RNA run, I would recommend the following:

python3 -m reditools analyze \
    "$BAM" \
    --output-file "$OUTPUT_DIR/${SAMPLE}_dna_reditools.txt" \
    --min-read-quality 255 \
    --min-edits 0 \
    --min-read-depth 5 \
    --threads 24 \
    --dna

I used long-hand for the arguments for clarity here. The two major changes here are --min-edits 0 and --min-read-depth 5. Ideally, you want to find elements of your RNAseq that differ from your DNA, so you want to include unedited sites from the DNA data. The --min-read-depth 5 should give you reasonable coverage, though I have seen some people use values as high as 15.

Note: There is a known bug in v3.6 when setting --min-edits 0. This bug has been patched, but you need to download the development version of the code or v3.7+ (when available).

After you have both the DNA and RNA results, you can use the annotate tool. Again, using what code you've provided.

python3 -m reditools annotate \
    "$OUTPUT_DIR/${SAMPLE}_reditools.txt" \
    "$OUTPUT_DIR/${SAMPLE}_dna_reditools.txt" \
    --bam "$BAM"

The --bam "$BAM" flag is optional, but will speed up your results.

Note: There is a known bug in v3.6 with the annotate tool caused by chromosomes/contigs missing from the DNA/RNA file. This bug has been patched, but you need to download the development version of the code or v3.7+ (when available).

[jonng1000]

Ok thank you so much, I will try to git clone from the REDItools repo and see if it works =)

[jonng1000]

Hi, I've tried to run REDItools on my DNA BAM files, using GATK best practices workflow, which uses BQSR. However, the first error I got is this:

jonathan@comp1:~/work/rna_edit/pancreatic_lines/elly_dataset/jon_workflow/dna_sites_wes/gatk_workflow$ python3 -m reditools analyze ~/work/rna_edit/pancreatic_lines/elly_dataset/jon_workflow/wes_align/bqsr_output/SRR8619006.bqsr.bam -N -o test_output.txt -q 30 -l 10 -me 0 -t 24 

[ERROR] (<class 'IndexError'>) list index out of range 
[ERROR] (<class 'IndexError'>) list index out of range 
[ERROR] (<class 'IndexError'>) list index out of range 
[ERROR] Killing job 

From what I understand, the GATK ref genome has HLA contigs, which has names like HLA-A*01:01. The GENCODE ref genome does not have such contig names and only has chr names like chr1, chr2 etc. It appears REDItools has been hardcoded to only accept chr names from GENCODE and Ensembl ref genomes. Would it be possible for REDItools to be accept the GATK ref genome as well?

I solved the first error by removing all HLA contigs from my DNA bam fille header, as well as all reads mapping to the HLA contigs. However, I'm not sure how robust this solution is.

Secondly, I hit another error after solving the first one. This is similar to the first error as I had a KeyError when running REDItools, as it doesn't recognize 'R'. Apparently REDItools only recognises the 4 bases in the ref genome, but the GATK ref genome uses more letters besides the 4 main ones, as it has ambiguous bases to represent variation. Thus, for now, I decided not to do BQSR as that uses the GATK ref genome. Hence, would it also be possible for me to request for REDItools to accept letters which arent the 4 main ones in the ref genome?

If any of my understanding/logic here is wrong, please correct me. Thank you =).

[ahanden]

From what I understand, the GATK ref genome has HLA contigs, which has names like HLA-A*01:01. The GENCODE ref genome does not have such contig names and only has chr names like chr1, chr2 etc. It appears REDItools has been hardcoded to only accept chr names from GENCODE and Ensembl ref genomes. Would it be possible for REDItools to be accept the GATK ref genome as well?

To be clear, REDItools is not harcoded for chr names. What you've noticed is that if a chromosome name from the BAM file is not in the reference genome, then REDItools tries adding "chr" to the chromosome name. This saves some headaches when a BAM file labels a chromosome as "1" while the reference FASTA uses "chr1".

The problem you have is that the reference FASTA file you are giving REDItools does not have the HLA-A*01:01 chromosomes. Double check that you are using the correct filenames.

Secondly, I hit another error after solving the first one. This is similar to the first error as I had a KeyError when running REDItools, as it doesn't recognize 'R'. Apparently REDItools only recognises the 4 bases in the ref genome, but the GATK ref genome uses more letters besides the 4 main ones, as it has ambiguous bases to represent variation. Thus, for now, I decided not to do BQSR as that uses the GATK ref genome. Hence, would it also be possible for me to request for REDItools to accept letters which arent the 4 main ones in the ref genome?

There are currently no plans for including IUPAC ambiguity codes in REDItools analysis. If enough users show interest, I may add it to the feature todo list. But there are a number of complications that would come from trying to incorporate ambiguity codes, especially with stranded analysis.

If you still want to use BQSR, I think you can solve both of your problems here by using a reference genome that does not include ambiguity codes (e.g. GENCODE).

[jonng1000]

Ok thanks. By the way, may I please clarify, "Marking duplicates is particularly important for DNA sequencing data"? I've run my reditools annotation workflow on 2 different DNA alignment workflows. One using mark duplicates and one without. I thought that using mark duplicates will give me better results but it seems the opposite has occured. The number of edit sites I call after passing the DNA filter has dropped to very low levels. Eg for SRR8616212, I have 178 sites after DNA filtering (no mark duplicates). However, I only have 68 sites after DNA filtering (with mark duplicates). I am now worried about sensitivity and hence am thinking about dropping mark duplicates.

I note that the reditools 2020 protoccol does not use mark duplicates and your advice in #63 says mark duplicates doesn't really help. Hence I'm confused. May I please ask, for the purpose of using reditools on DNA samples for filtering RNA edit sites, is mark duplicates essential or optional? Thank you =)

[ahanden]

may I please clarify, "Marking duplicates is particularly important for DNA sequencing data"?

Thank you for bringing up #63. I had forgotten about that paper (https://pmc.ncbi.nlm.nih.gov/articles/PMC4965708/ for future readers). The authors show that removing PCR duplicates has a minimal effect on variant calling (Despite this finding, GATK best practices still recommends marking duplicates). So why does it have an effect on your analysis?

My best guess would be the --min-read-depth 5 parameter is the root issue. By flagging duplicates in your BAM file, some regions will have lower coverage. If there are fewer than 5 reads, in this case, edits will not be reported.

You can test this by looking at the edit sites missing from your duplicate-marked analysis. Let's say you have a line like this in the analysis output.

Region  Position  Reference  Strand  Coverage  MeanQ  BaseCount[A,C,G,T]  AllSubs  Frequency
chr1    15975     A          *       5         36.00  [4, 0, 1, 0]        AG       0.20

If you then use samtools to count the read coverage for that base position, you should get at least 5 reads.

samtools view --count "$BAM" chr1:15975-15976

Now if you filter out your duplicates, my guess is you will get fewer than 5 reads.

samtools view --count --exclude-flags 1024 "$BAM" chr1:15975-15976

May I please ask, for the purpose of using reditools on DNA samples for filtering RNA edit sites, is mark duplicates essential or optional?

Based on the conclusions of Ebbert et al, my opinion is that marking duplicates is optional. You can remove duplicates and justify your actions from GATK, or keep potential duplicates and reference Ebbert. The choice is yours.

[jonng1000]

I see ok, thanks =). I think I will drop mark duplicates as it causes me to lose a lot of edit sites. I've also decided not to use BQSR, as it uses variant ref files from the GATK reference set. I could in theory mix the GENCODE ref genome with the GATK variant ref files, but I'm concerned if mixing files from different reference sources is wise, so I decided to drop BQSR.