Request for a detailed tutorial on REDItools 3

[jonng1000]

Hi thank you for creating REDItools 3, it certainly sounds very useful. I'm new to RNA editing and am very keen to use your software. However I'm at a lost on how to run it. I did the protocol paper at https://pubmed.ncbi.nlm.nih.gov/31996844/ but understand its old and not all of the code is transferrable. Hence could I please request for a detailed tutorial with the specific code to run, to be given for beginers to learn and run REDItools 3? I am interested to run REDInet after REDItools 3

I've some questions which I hope I could get help in, either here or through the requested tutorial. I understand the workflow to use REDItools 3 should be:

Q1) RNA preprocessing - i understand that standard tools like fastp or fastqc + trimmomatic can be used. May I know if their default parameters can be used? Asking because I'm not sure if REDItools 3 requires special handling of data.

Q2) Alignment of reads - i understand STAR can be used. Is it's default parameters acceptable?

Q3) Do I need to identify the strand of the BAM files?

Q4) I understand the core command is analyze. Please correct me if I am wrong. But before it, do I need to run find-repeats? If so, may I please know what is exact comamnd here? What is find-repeat's output and how do I send it into analyze?

Also, I understand that find-repeats only identifies a homopolymeric regions, and that I need to mask other types of repeats. May I please know exactly how to do this masking?

Q5) I understand the next command is to use analyze. What is the exact command for this?

Q6) Once analyze has completed, can I use its output for REDInet? As in, can I follow the tutorial from, "Compress and Tabix indexing the REDItools output tables" (https://github.com/BioinfoUNIBA/REDInet) onwards?

[ahanden]

A tutorial is definitely needed - and I'm actively writing one now. With perseverance, it should be done by the end of March.

On to your questions.

Q1) RNA preprocessing - i understand that standard tools like fastp or fastqc + trimmomatic can be used. May I know if their default parameters can be used? Asking because I'm not sure if REDItools 3 requires special handling of data.

You can and should use adapter trimming as you would for a typical NGS workflow. No special considerations need to be taken for REDItools.

Q2) Alignment of reads - i understand STAR can be used. Is it's default parameters acceptable?

With any aligner, I strongly encourage using whatever parameters are needed to include MD tags in the output. A BAM file with MD tags will be processed many times faster than one without. For STAR specifically, you'll need this: --outSAMattributes NH HI AS nM MD

Q3) Do I need to identify the strand of the BAM files?

REDItools can be run in unstranded mode (default) or stranded and can determine strand from the BAM file.

Q4) I understand the core command is analyze. Please correct me if I am wrong. But before it, do I need to run find-repeats? If so, may I please know what is exact comamnd here? What is find-repeat's output and how do I send it into analyze?

The find-repeats tool is legacy code. You do not need to run it. If you want to focus your analysis on repetitive elements (or exclude them), I suggest using Repeat Masker data from UCSC.

Q5) I understand the next command is to use analyze. What is the exact command for this?

This is my best translation of the REDI1 parameters from the REDInet tutorial:

python3 -m reditools analyze -l 1 -m 255 -me 1 -bq 30 -Men 1 -s 2 -C.

You may want to use --threads and --window for parallel processing. And if you want to exclude repetitive regions, the -k option can be used with a BED file.

Q6) Once analyze has completed, can I use its output for REDInet? As in, can I follow the tutorial from, "Compress and Tabix indexing the REDItools output tables" (https://github.com/BioinfoUNIBA/REDInet) onwards?

As far as I'm aware, yes, you can go by the REDInet tutorial. My role in the Briefings in bioinformatics paper was exclusively on REDItools development, so any questions regarding REDInet should be put in that repo.

[jonng1000]

Thank so you much for your prompt reply! I appreciate it a lot, especially since I have a tight deadline on running reditools 3 successfully =). If you are interested, I'm happy to read a draft of the tutorial and give feedback. I'm also happy to contribute to writing it. I come from a bio background with some experience in programming, but am new to RNA editing. Hence I could give my perspective of how a new user would read and understand REDItools 3.

From your answers, I infer that since REDItools 3 has been updated, I no longer need to identify strandedness and repetitive elements before running it, as it is able to handle such things on its own. Hence I will not do these steps. Please correct me if I'm wrong. I will explore Repeat Masker only if I have specific reasons to look at/exclude repetitive elements.

Thank you for giving me sample code. Would the following code be correct given a set of input and output files, for a general workflow?
python3 -m reditools analyze <BAM file> -r <hg38_reference.fasta.gz> -o <reditools_output_table.txt> -q 255 -s 2 -C True

I made the following modifications from the sample code:

• Removed arguments where the value matches the default in the REDItools 3 help output in the command line
• Sample code has -m 255 -> realised this maps to -q MIN_READ_QUALITY in REDItools 3
• Sample code has -Men 1 -> removed this to use REDItools 3 default, which is 4
• Sample code has -s 2 -> remove this to use REDItools 3 default, which is 0. This means unstranded, which I understand is the default setting
• Sample code has -C, realised I need to put True here, based on the REDItools 3 help output

For my STAR aligned bam files, if I happen to select only uniquely mapped reads into a new bam file, would this be similar to using Repeat Masker? I know these two would not be exactly the same, but uniquely mapped reads should not map to repetitive elements in the genome, hence I would be able to get rid of such elements too right? I am considering doing this before sending the bam file to REDItools 3.

Thank you for explaining more about REDInet. I will contact them if I need help with it

[ahanden]

I'm happy to read a draft of the tutorial and give feedback

That would be great! You can send me an email and we can correspond.

I no longer need to identify strandedness and repetitive elements before running it, as it is able to handle such things on its own.

Regarding strandedness - I apologize. I didn't write my response correctly. You do still need to determine the strandedness of your BAM file and supply that using the --strand option. Your understanding of repetitive elements is correct.

Thank you for giving me sample code. Would the following code be correct given a set of input and output files, for a general workflow?
python3 -m reditools analyze <BAM file> -r <hg38_reference.fasta.gz> -o <reditools_output_table.txt> -q 255 -s 2 -C True

If your BAM file has MD tags, do not provide a reference with -r. The reference option will override the presence of MD tags and will slow REDItools down significantly. You also do not need to add "True" after the -C option. Simply adding -C should turn on strand correction.

For my STAR aligned bam files, if I happen to select only uniquely mapped reads into a new bam file, would this be similar to using Repeat Masker?

Since it looks like you're using human data, I would not recommend this. Simply get Repeat Masker data directly from UCSC: https://genome.ucsc.edu/cgi-bin/hgTables

[jonng1000]

Thanks =) have sent you an email.

I see, I understand that strandedness needs to be determined and will figure out how to do it with RSeQC. Thanks for the additonal advice. Have updated my code as follows, and my BAM files do have the MD tags.

python3 -m reditools analyze <BAM file> -o <reditools_output_table.txt> -q 255 -s <0/1/2 here depending on strand> -C -t <# of threads>

I've added -t since I have multiple cores, and I think its good to show this in the example code in the tutorial, since many bioinfo users can be assumed to have access to a multicore computer. Not setting the -w argument as I will let it use the default value of 0. Hope my run works =)

[jonng1000]

From RSeQC, I understand its strandedness output is 1++, 1--, 2+-, 2-+ for the forward strand and 1+-, 1-+, 2++, 2-- for the reverse strand. How does this map to reditool's strandedness value for the -s parameter? Does the forward strand map to, "1 (secondstrand oriented)" and the reverse strand map to, "2 (firststrand oriented)"?

[jonng1000]

Another question, do my BAM files need to have duplicates marked by using MarkDuplicates from GATK Picard?

[ahanden]

From RSeQC, I understand its strandedness output is 1++, 1--, 2+-, 2-+ for the forward strand and 1+-, 1-+, 2++, 2-- for the reverse strand. How does this map to reditool's strandedness value for the -s parameter? Does the forward strand map to, "1 (secondstrand oriented)" and the reverse strand map to, "2 (firststrand oriented)"?

Short answer: here is how RSeQC corresponds to REDItools

• 1++, 1--, 2+-, 2-+ => --strand 1
• 1+-, 1-+, 2++, 2-- => --strand 2

Some clarification on how the strand option works under the hood.

Unstranded (0) ignores any strand data from the BAM file, essentially assuming all reads are on the forward strand. The downside of this method is that if you do have a reverse reads, the detected edit would be inverse (i.e. An A-to-I would show up as T-to-G).

The other two modes make use of three SAM flags to determine strandedness: 16 (read reverse strand), 64 (first in pair), and 128 (second in pair).

In secondstrand oriented (1) forward strand is defined as first in pair and not reversed OR second in pair and reversed.

In firststrand oriented (2) forward strand is defined as first in pair and reversed OR second in pair and not reversed.

In either case, reverse reads have edits reported as their complement.

I think our wording choice for "firststrand" and "secondstrand" is really confusing. Admittedly, they've been copy pasted from REDI2, which copied from REDI1.

Another question, do my BAM files need to have duplicates marked by using MarkDuplicates from GATK Picard?

No, you do not need to mark duplicates and you should not filter duplicates. A lot of papers do filter duplicates in RNAseq, but the reality is that in RNAseq, there is no way to differentiate between a PCR duplicate and actual transcript. The bioinformatics community has come to a consensus at this point that you should not mark or remove duplicates from RNAseq.

[jonng1000]

Ok thank you, I understand now =)

[jonng1000]

Sorry, I've another question. May I please know what sort of post processing is required after edit site calling? I understand that post processing of reditools' output is done by people in the field. I understand this includes filtering for a certain number of coverage, edited reads, homopolymeric regions, and also for polymorphic regions. If there's stuff which I missed out, please let me know.

Perhaps I could get these instructions from you, and we can write them in the tutorial for others to follow =)

[ahanden]

I haven't done a real review of how people are post-processing REDItools output. But since the results are simple TSVs, you can do a lot with awk. I do recommend skipping "problematic" sections of the genome (I like the DAC exclusion list for this). And you can reduce your search area depending on the mechanism you're interested in. For example, ADAR1 likes to target Alu repeats. You can download the coordinates of those from UCSC's Repeat Masker table.

[jonng1000]

I see ok thank you. Will explore this. Perhaps the tutorial could have 3 sections. First, how to preprocess RNA sequence data for reditools (I can help with this). Second, how to use REDItools on the preprocessed data (I can give feedback from a perspective of an end user). Third, an example application, such as how to use REDItools to identify A-to-I edit sites and filter them computationally (I'm working in this area, so I can help with this, and read up on the strategies you mentioned here)

[jonng1000]

I've another question on the proper interpretation of REDItools output. Lets say I want to identify A-to-I edit sites and see the below output, when REDItools was ran with the -s 2 parameter.

Region  Position  Reference  Strand  Coverage  MeanQ  BaseCount[A,C,G,T]  AllSubs  Frequency  gCoverage ...
1       14522     C          -       9         41.00  [0,2,0,7]           CT      0.78       -
1       14523     A          -       9         41.00  [2,0,7,0]           AG      0.78       -

RNA editing is on mRNA, which is complementary to DNA. However, I understand REDItools reports all the bases in the ouput above, in terms of DNA bases. Thus would it be correct for me to say the following?

RNA edits (mRNA strand) were detected on the DNA - strand. This means that for position 14523, the DNA - strand would have the equivalent of a TC change. Thus, since Region and Position values are given with respect to the + strand, I am looking for the complent of TC on the + strand, which would be AG on the + strand. Since I is interpreted as G by NGS technology, hence chr1:14523 would contain my edit site.

If the above is true, then what happens when I see edit sites given on the + strand?

Region  Position  Reference  Strand  Coverage  MeanQ  BaseCount[A,C,G,T]  AllSubs  Frequency  gCoverage ...
1    829009    T    +    41    39.85    [0, 1, 0, 40]    TC    0.02    -

Would the following be true?

T-to-C RNA edit (mRNA strand) was detected on the DNA + strand. Thus, I am looking for the complent of TC on the mRNA strand, which would be AG. Hence, chr1:829009 would also contain my A-to-I edit site.

If all the above is true, does that mean that to filter for A-to-I edit sites only, I need to select AG values for AllSubs when the Strand is -, and TC values for AllSubs when the Strand is +?

Thank you so much for your patience, and sorry for the trouble.

[ahanden]

T-to-C RNA edit (mRNA strand) was detected on the DNA + strand. Thus, I am looking for the complent of TC on the mRNA strand, which would be AG. Hence, chr1:829009 would also contain my A-to-I edit site.

One of the specific advantages of using REDItools is that you don't have to do this kind of logic. When run in stranded mode, REDItools reports the edits relative to the RNA, not the genome, so you only need to filter for AG. The only exception is when the strand is ambiguous (*). In that case, if you see an AG variant, you can't really know whether it is A-to-I or T-to-C.

An aside - I just release v3.6 today.

[jonng1000]

Thanks, will run pip install to update my REDItools 3. Re my main question, ok I understand that AllSubs is relative to the RNA. But how do I interpret the "Position" and "Reference" columns then? These are relative to DNA right?

So for:

Region  Position  Reference  Strand  Coverage  MeanQ  BaseCount[A,C,G,T]  AllSubs  Frequency  gCoverage ...
1       14523     A          -       9         41.00  [2,0,7,0]           AG      0.78       -

The reference base at the genome is A? And the mRNA strand is detected at the - strand of the DNA, and has a AG subsitution?

This is because the REDItools 3 readme says that the Reference column is the base at the reference sequence, which I assume would be DNA.